dist/Math-BigInt/lib/Math/BigFloat.pm

   1 package Math::BigFloat;
   2
   3 #
   4 # Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
   5 #
   6
   7 # The following hash values are internally used:
   8 #   _e  : exponent (ref to $CALC object)
   9 #   _m  : mantissa (ref to $CALC object)
  10 #   _es : sign of _e
  11 # sign  : +,-,+inf,-inf, or "NaN" if not a number
  12 #   _a  : accuracy
  13 #   _p  : precision
  14
  15 $VERSION = '1.99_04';
  16 require 5.006002;
  17
  18 require Exporter;
  19 @ISA            = qw/Math::BigInt/;
  20 @EXPORT_OK      = qw/bpi/;
  21
  22 use strict;
  23 # $_trap_inf/$_trap_nan are internal and should never be accessed from outside
  24 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode
  25             $upgrade $downgrade $_trap_nan $_trap_inf/;
  26 my $class = "Math::BigFloat";
  27
  28 use overload
  29 '<=>'   =>      sub { my $rc = $_[2] ?
  30                       ref($_[0])->bcmp($_[1],$_[0]) :
  31                       ref($_[0])->bcmp($_[0],$_[1]);
  32                       $rc = 1 unless defined $rc;
  33                       $rc <=> 0;
  34                 },
  35 # we need '>=' to get things like "1 >= NaN" right:
  36 '>='    =>      sub { my $rc = $_[2] ?
  37                       ref($_[0])->bcmp($_[1],$_[0]) :
  38                       ref($_[0])->bcmp($_[0],$_[1]);
  39                       # if there was a NaN involved, return false
  40                       return '' unless defined $rc;
  41                       $rc >= 0;
  42                 },
  43 'int'   =>      sub { $_[0]->as_number() },             # 'trunc' to bigint
  44 ;
  45
  46 ##############################################################################
  47 # global constants, flags and assorted stuff
  48
  49 # the following are public, but their usage is not recommended. Use the
  50 # accessor methods instead.
  51
  52 # class constants, use Class->constant_name() to access
  53 # one of 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'
  54 $round_mode = 'even';
  55 $accuracy   = undef;
  56 $precision  = undef;
  57 $div_scale  = 40;
  58
  59 $upgrade = undef;
  60 $downgrade = undef;
  61 # the package we are using for our private parts, defaults to:
  62 # Math::BigInt->config()->{lib}
  63 my $MBI = 'Math::BigInt::FastCalc';
  64
  65 # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
  66 $_trap_nan = 0;
  67 # the same for infinity
  68 $_trap_inf = 0;
  69
  70 # constant for easier life
  71 my $nan = 'NaN';
  72
  73 my $IMPORT = 0; # was import() called yet? used to make require work
  74
  75 # some digits of accuracy for blog(undef,10); which we use in blog() for speed
  76 my $LOG_10 =
  77  '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
  78 my $LOG_10_A = length($LOG_10)-1;
  79 # ditto for log(2)
  80 my $LOG_2 =
  81  '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
  82 my $LOG_2_A = length($LOG_2)-1;
  83 my $HALF = '0.5';                       # made into an object if nec.
  84
  85 ##############################################################################
  86 # the old code had $rnd_mode, so we need to support it, too
  87
  88 sub TIESCALAR   { my ($class) = @_; bless \$round_mode, $class; }
  89 sub FETCH       { return $round_mode; }
  90 sub STORE       { $rnd_mode = $_[0]->round_mode($_[1]); }
  91
  92 BEGIN
  93   {
  94   # when someone sets $rnd_mode, we catch this and check the value to see
  95   # whether it is valid or not.
  96   $rnd_mode   = 'even'; tie $rnd_mode, 'Math::BigFloat';
  97
  98   # we need both of them in this package:
  99   *as_int = \&as_number;
 100   }
 101
 102 ##############################################################################
 103
 104 {
 105   # valid method aliases for AUTOLOAD
 106   my %methods = map { $_ => 1 }
 107    qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
 108         fint facmp fcmp fzero fnan finf finc fdec ffac fneg
 109         fceil ffloor frsft flsft fone flog froot fexp
 110       /;
 111   # valid methods that can be handed up (for AUTOLOAD)
 112   my %hand_ups = map { $_ => 1 }
 113    qw / is_nan is_inf is_negative is_positive is_pos is_neg
 114         accuracy precision div_scale round_mode fabs fnot
 115         objectify upgrade downgrade
 116         bone binf bnan bzero
 117         bsub
 118       /;
 119
 120   sub _method_alias { exists $methods{$_[0]||''}; }
 121   sub _method_hand_up { exists $hand_ups{$_[0]||''}; }
 122 }
 123
 124 ##############################################################################
 125 # constructors
 126
 127 sub new
 128   {
 129   # create a new BigFloat object from a string or another bigfloat object.
 130   # _e: exponent
 131   # _m: mantissa
 132   # sign  => sign (+/-), or "NaN"
 133
 134   my ($class,$wanted,@r) = @_;
 135
 136   # avoid numify-calls by not using || on $wanted!
 137   return $class->bzero() if !defined $wanted;   # default to 0
 138   return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
 139
 140   $class->import() if $IMPORT == 0;             # make require work
 141
 142   my $self = {}; bless $self, $class;
 143   # shortcut for bigints and its subclasses
 144   if ((ref($wanted)) && UNIVERSAL::can( $wanted, "as_number"))
 145     {
 146     $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy
 147     $self->{_e} = $MBI->_zero();
 148     $self->{_es} = '+';
 149     $self->{sign} = $wanted->sign();
 150     return $self->bnorm();
 151     }
 152   # else: got a string or something masquerading as number (with overload)
 153
 154   # handle '+inf', '-inf' first
 155   if ($wanted =~ /^[+-]?inf\z/)
 156     {
 157     return $downgrade->new($wanted) if $downgrade;
 158
 159     $self->{sign} = $wanted;            # set a default sign for bstr()
 160     return $self->binf($wanted);
 161     }
 162
 163   # shortcut for simple forms like '12' that neither have trailing nor leading
 164   # zeros
 165   if ($wanted =~ /^([+-]?)([1-9][0-9]*[1-9])$/)
 166     {
 167     $self->{_e} = $MBI->_zero();
 168     $self->{_es} = '+';
 169     $self->{sign} = $1 || '+';
 170     $self->{_m} = $MBI->_new($2);
 171     return $self->round(@r) if !$downgrade;
 172     }
 173
 174   my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
 175   if (!ref $mis)
 176     {
 177     if ($_trap_nan)
 178       {
 179       require Carp;
 180       Carp::croak ("$wanted is not a number initialized to $class");
 181       }
 182
 183     return $downgrade->bnan() if $downgrade;
 184
 185     $self->{_e} = $MBI->_zero();
 186     $self->{_es} = '+';
 187     $self->{_m} = $MBI->_zero();
 188     $self->{sign} = $nan;
 189     }
 190   else
 191     {
 192     # make integer from mantissa by adjusting exp, then convert to int
 193     $self->{_e} = $MBI->_new($$ev);             # exponent
 194     $self->{_es} = $$es || '+';
 195     my $mantissa = "$$miv$$mfv";                # create mant.
 196     $mantissa =~ s/^0+(\d)/$1/;                 # strip leading zeros
 197     $self->{_m} = $MBI->_new($mantissa);        # create mant.
 198
 199     # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
 200     if (CORE::length($$mfv) != 0)
 201       {
 202       my $len = $MBI->_new( CORE::length($$mfv));
 203       ($self->{_e}, $self->{_es}) =
 204         _e_sub ($self->{_e}, $len, $self->{_es}, '+');
 205       }
 206     # we can only have trailing zeros on the mantissa if $$mfv eq ''
 207     else
 208       {
 209       # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
 210       # because that is faster, especially when _m is not stored in base 10.
 211       my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/;
 212       if ($zeros != 0)
 213         {
 214         my $z = $MBI->_new($zeros);
 215         # turn '120e2' into '12e3'
 216         $MBI->_rsft ( $self->{_m}, $z, 10);
 217         ($self->{_e}, $self->{_es}) =
 218           _e_add ( $self->{_e}, $z, $self->{_es}, '+');
 219         }
 220       }
 221     $self->{sign} = $$mis;
 222
 223     # for something like 0Ey, set y to 1, and -0 => +0
 224     # Check $$miv for being '0' and $$mfv eq '', because otherwise _m could not
 225     # have become 0. That's faster than to call $MBI->_is_zero().
 226     $self->{sign} = '+', $self->{_e} = $MBI->_one()
 227      if $$miv eq '0' and $$mfv eq '';
 228
 229     return $self->round(@r) if !$downgrade;
 230     }
 231   # if downgrade, inf, NaN or integers go down
 232
 233   if ($downgrade && $self->{_es} eq '+')
 234     {
 235     if ($MBI->_is_zero( $self->{_e} ))
 236       {
 237       return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
 238       }
 239     return $downgrade->new($self->bsstr());
 240     }
 241   $self->bnorm()->round(@r);                    # first normalize, then round
 242   }
 243
 244 sub copy
 245   {
 246   # if two arguments, the first one is the class to "swallow" subclasses
 247   if (@_ > 1)
 248     {
 249     my  $self = bless {
 250         sign => $_[1]->{sign},
 251         _es => $_[1]->{_es},
 252         _m => $MBI->_copy($_[1]->{_m}),
 253         _e => $MBI->_copy($_[1]->{_e}),
 254     }, $_[0] if @_ > 1;
 255
 256     $self->{_a} = $_[1]->{_a} if defined $_[1]->{_a};
 257     $self->{_p} = $_[1]->{_p} if defined $_[1]->{_p};
 258     return $self;
 259     }
 260
 261   my $self = bless {
 262         sign => $_[0]->{sign},
 263         _es => $_[0]->{_es},
 264         _m => $MBI->_copy($_[0]->{_m}),
 265         _e => $MBI->_copy($_[0]->{_e}),
 266         }, ref($_[0]);
 267
 268   $self->{_a} = $_[0]->{_a} if defined $_[0]->{_a};
 269   $self->{_p} = $_[0]->{_p} if defined $_[0]->{_p};
 270   $self;
 271   }
 272
 273 sub _bnan
 274   {
 275   # used by parent class bone() to initialize number to NaN
 276   my $self = shift;
 277
 278   if ($_trap_nan)
 279     {
 280     require Carp;
 281     my $class = ref($self);
 282     Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
 283     }
 284
 285   $IMPORT=1;                                    # call our import only once
 286   $self->{_m} = $MBI->_zero();
 287   $self->{_e} = $MBI->_zero();
 288   $self->{_es} = '+';
 289   }
 290
 291 sub _binf
 292   {
 293   # used by parent class bone() to initialize number to +-inf
 294   my $self = shift;
 295
 296   if ($_trap_inf)
 297     {
 298     require Carp;
 299     my $class = ref($self);
 300     Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
 301     }
 302
 303   $IMPORT=1;                                    # call our import only once
 304   $self->{_m} = $MBI->_zero();
 305   $self->{_e} = $MBI->_zero();
 306   $self->{_es} = '+';
 307   }
 308
 309 sub _bone
 310   {
 311   # used by parent class bone() to initialize number to 1
 312   my $self = shift;
 313   $IMPORT=1;                                    # call our import only once
 314   $self->{_m} = $MBI->_one();
 315   $self->{_e} = $MBI->_zero();
 316   $self->{_es} = '+';
 317   }
 318
 319 sub _bzero
 320   {
 321   # used by parent class bone() to initialize number to 0
 322   my $self = shift;
 323   $IMPORT=1;                                    # call our import only once
 324   $self->{_m} = $MBI->_zero();
 325   $self->{_e} = $MBI->_one();
 326   $self->{_es} = '+';
 327   }
 328
 329 sub isa
 330   {
 331   my ($self,$class) = @_;
 332   return if $class =~ /^Math::BigInt/;          # we aren't one of these
 333   UNIVERSAL::isa($self,$class);
 334   }
 335
 336 sub config
 337   {
 338   # return (later set?) configuration data as hash ref
 339   my $class = shift || 'Math::BigFloat';
 340
 341   if (@_ == 1 && ref($_[0]) ne 'HASH')
 342     {
 343     my $cfg = $class->SUPER::config();
 344     return $cfg->{$_[0]};
 345     }
 346
 347   my $cfg = $class->SUPER::config(@_);
 348
 349   # now we need only to override the ones that are different from our parent
 350   $cfg->{class} = $class;
 351   $cfg->{with} = $MBI;
 352   $cfg;
 353   }
 354
 355 ##############################################################################
 356 # string conversion
 357
 358 sub bstr
 359   {
 360   # (ref to BFLOAT or num_str ) return num_str
 361   # Convert number from internal format to (non-scientific) string format.
 362   # internal format is always normalized (no leading zeros, "-0" => "+0")
 363   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
 364
 365   if ($x->{sign} !~ /^[+-]$/)
 366     {
 367     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
 368     return 'inf';                                       # +inf
 369     }
 370
 371   my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
 372
 373   # $x is zero?
 374   my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
 375   if ($not_zero)
 376     {
 377     $es = $MBI->_str($x->{_m});
 378     $len = CORE::length($es);
 379     my $e = $MBI->_num($x->{_e});
 380     $e = -$e if $x->{_es} eq '-';
 381     if ($e < 0)
 382       {
 383       $dot = '';
 384       # if _e is bigger than a scalar, the following will blow your memory
 385       if ($e <= -$len)
 386         {
 387         my $r = abs($e) - $len;
 388         $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
 389         }
 390       else
 391         {
 392         substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
 393         $cad = -$cad if $x->{_es} eq '-';
 394         }
 395       }
 396     elsif ($e > 0)
 397       {
 398       # expand with zeros
 399       $es .= '0' x $e; $len += $e; $cad = 0;
 400       }
 401     } # if not zero
 402
 403   $es = '-'.$es if $x->{sign} eq '-';
 404   # if set accuracy or precision, pad with zeros on the right side
 405   if ((defined $x->{_a}) && ($not_zero))
 406     {
 407     # 123400 => 6, 0.1234 => 4, 0.001234 => 4
 408     my $zeros = $x->{_a} - $cad;                # cad == 0 => 12340
 409     $zeros = $x->{_a} - $len if $cad != $len;
 410     $es .= $dot.'0' x $zeros if $zeros > 0;
 411     }
 412   elsif ((($x->{_p} || 0) < 0))
 413     {
 414     # 123400 => 6, 0.1234 => 4, 0.001234 => 6
 415     my $zeros = -$x->{_p} + $cad;
 416     $es .= $dot.'0' x $zeros if $zeros > 0;
 417     }
 418   $es;
 419   }
 420
 421 sub bsstr
 422   {
 423   # (ref to BFLOAT or num_str ) return num_str
 424   # Convert number from internal format to scientific string format.
 425   # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
 426   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
 427
 428   if ($x->{sign} !~ /^[+-]$/)
 429     {
 430     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
 431     return 'inf';                                       # +inf
 432     }
 433   my $sep = 'e'.$x->{_es};
 434   my $sign = $x->{sign}; $sign = '' if $sign eq '+';
 435   $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
 436   }
 437
 438 sub numify
 439   {
 440   # Make a number from a BigFloat object
 441   # simple return a string and let Perl's atoi()/atof() handle the rest
 442   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
 443   $x->bsstr();
 444   }
 445
 446 ##############################################################################
 447 # public stuff (usually prefixed with "b")
 448
 449 sub bneg
 450   {
 451   # (BINT or num_str) return BINT
 452   # negate number or make a negated number from string
 453   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
 454
 455   return $x if $x->modify('bneg');
 456
 457   # for +0 dont negate (to have always normalized +0). Does nothing for 'NaN'
 458   $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
 459   $x;
 460   }
 461
 462 # tels 2001-08-04
 463 # XXX TODO this must be overwritten and return NaN for non-integer values
 464 # band(), bior(), bxor(), too
 465 #sub bnot
 466 #  {
 467 #  $class->SUPER::bnot($class,@_);
 468 #  }
 469
 470 sub bcmp
 471   {
 472   # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
 473
 474   # set up parameters
 475   my ($self,$x,$y) = (ref($_[0]),@_);
 476
 477   # objectify is costly, so avoid it
 478   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 479     {
 480     ($self,$x,$y) = objectify(2,@_);
 481     }
 482
 483   return $upgrade->bcmp($x,$y) if defined $upgrade &&
 484     ((!$x->isa($self)) || (!$y->isa($self)));
 485
 486   # Handle all 'nan' cases.
 487
 488   return undef if ($x->{sign} eq $nan) || ($y->{sign} eq $nan);
 489
 490   # Handle all '+inf' and '-inf' cases.
 491
 492   return  0 if ($x->{sign} eq '+inf' && $y->{sign} eq '+inf' ||
 493                 $x->{sign} eq '-inf' && $y->{sign} eq '-inf');
 494   return +1 if $x->{sign} eq '+inf';    # x = +inf and y < +inf
 495   return -1 if $x->{sign} eq '-inf';    # x = -inf and y > -inf
 496   return -1 if $y->{sign} eq '+inf';    # x < +inf and y = +inf
 497   return +1 if $y->{sign} eq '-inf';    # x > -inf and y = -inf
 498
 499   # Handle all cases with opposite signs.
 500
 501   return +1 if $x->{sign} eq '+' && $y->{sign} eq '-';  # also does 0 <=> -y
 502   return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';  # also does -x <=> 0
 503
 504   # Handle all remaining zero cases.
 505
 506   my $xz = $x->is_zero();
 507   my $yz = $y->is_zero();
 508   return  0 if $xz && $yz;                              # 0 <=> 0
 509   return -1 if $xz && $y->{sign} eq '+';                # 0 <=> +y
 510   return +1 if $yz && $x->{sign} eq '+';                # +x <=> 0
 511
 512   # Both arguments are now finite, non-zero numbers with the same sign.
 513
 514   my $cmp;
 515
 516   # The next step is to compare the exponents, but since each mantissa is an
 517   # integer of arbitrary value, the exponents must be normalized by the length
 518   # of the mantissas before we can compare them.
 519
 520   my $mxl = $MBI->_len($x->{_m});
 521   my $myl = $MBI->_len($y->{_m});
 522
 523   # If the mantissas have the same length, there is no point in normalizing the
 524   # exponents by the length of the mantissas, so treat that as a special case.
 525
 526   if ($mxl == $myl) {
 527
 528       # First handle the two cases where the exponents have different signs.
 529
 530       if ($x->{_es} eq '+' && $y->{_es} eq '-') {
 531           $cmp = +1;
 532       }
 533
 534       elsif ($x->{_es} eq '-' && $y->{_es} eq '+') {
 535           $cmp = -1;
 536       }
 537
 538       # Then handle the case where the exponents have the same sign.
 539
 540       else {
 541           $cmp = $MBI->_acmp($x->{_e}, $y->{_e});
 542           $cmp = -$cmp if $x->{_es} eq '-';
 543       }
 544
 545       # Adjust for the sign, which is the same for x and y, and bail out if
 546       # we're done.
 547
 548       $cmp = -$cmp if $x->{sign} eq '-';        # 124 > 123, but -124 < -123
 549       return $cmp if $cmp;
 550
 551   }
 552
 553   # We must normalize each exponent by the length of the corresponding
 554   # mantissa. Life is a lot easier if we first make both exponents
 555   # non-negative. We do this by adding the same positive value to both
 556   # exponent. This is safe, because when comparing the exponents, only the
 557   # relative difference is important.
 558
 559   my $ex;
 560   my $ey;
 561
 562   if ($x->{_es} eq '+') {
 563
 564       # If the exponent of x is >= 0 and the exponent of y is >= 0, there is no
 565       # need to do anything special.
 566
 567       if ($y->{_es} eq '+') {
 568           $ex = $MBI->_copy($x->{_e});
 569           $ey = $MBI->_copy($y->{_e});
 570       }
 571
 572       # If the exponent of x is >= 0 and the exponent of y is < 0, add the
 573       # absolute value of the exponent of y to both.
 574
 575       else {
 576           $ex = $MBI->_copy($x->{_e});
 577           $ex = $MBI->_add($ex, $y->{_e});      # ex + |ey|
 578           $ey = $MBI->_zero();                  # -ex + |ey| = 0
 579       }
 580
 581   } else {
 582
 583       # If the exponent of x is < 0 and the exponent of y is >= 0, add the
 584       # absolute value of the exponent of x to both.
 585
 586       if ($y->{_es} eq '+') {
 587           $ex = $MBI->_zero();                  # -ex + |ex| = 0
 588           $ey = $MBI->_copy($y->{_e});
 589           $ey = $MBI->_add($ey, $x->{_e});      # ey + |ex|
 590       }
 591
 592       # If the exponent of x is < 0 and the exponent of y is < 0, add the
 593       # absolute values of both exponents to both exponents.
 594
 595       else {
 596           $ex = $MBI->_copy($y->{_e});          # -ex + |ey| + |ex| = |ey|
 597           $ey = $MBI->_copy($x->{_e});          # -ey + |ex| + |ey| = |ex|
 598       }
 599
 600   }
 601
 602   # Now we can normalize the exponents by adding lengths of the mantissas.
 603
 604   $MBI->_add($ex, $MBI->_new($mxl));
 605   $MBI->_add($ey, $MBI->_new($myl));
 606
 607   # We're done if the exponents are different.
 608
 609   $cmp = $MBI->_acmp($ex, $ey);
 610   $cmp = -$cmp if $x->{sign} eq '-';            # 124 > 123, but -124 < -123
 611   return $cmp if $cmp;
 612
 613   # Compare the mantissas, but first normalize them by padding the shorter
 614   # mantissa with zeros (shift left) until it has the same length as the longer
 615   # mantissa.
 616
 617   my $mx = $x->{_m};
 618   my $my = $y->{_m};
 619
 620   if ($mxl > $myl) {
 621       $my = $MBI->_lsft($MBI->_copy($my), $MBI->_new($mxl - $myl), 10);
 622   } elsif ($mxl < $myl) {
 623       $mx = $MBI->_lsft($MBI->_copy($mx), $MBI->_new($myl - $mxl), 10);
 624   }
 625
 626   $cmp = $MBI->_acmp($mx, $my);
 627   $cmp = -$cmp if $x->{sign} eq '-';            # 124 > 123, but -124 < -123
 628   return $cmp;
 629
 630   }
 631
 632 sub bacmp
 633   {
 634   # Compares 2 values, ignoring their signs.
 635   # Returns one of undef, <0, =0, >0. (suitable for sort)
 636
 637   # set up parameters
 638   my ($self,$x,$y) = (ref($_[0]),@_);
 639   # objectify is costly, so avoid it
 640   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 641     {
 642     ($self,$x,$y) = objectify(2,@_);
 643     }
 644
 645   return $upgrade->bacmp($x,$y) if defined $upgrade &&
 646     ((!$x->isa($self)) || (!$y->isa($self)));
 647
 648   # handle +-inf and NaN's
 649   if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
 650     {
 651     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 652     return 0 if ($x->is_inf() && $y->is_inf());
 653     return 1 if ($x->is_inf() && !$y->is_inf());
 654     return -1;
 655     }
 656
 657   # shortcut
 658   my $xz = $x->is_zero();
 659   my $yz = $y->is_zero();
 660   return 0 if $xz && $yz;                               # 0 <=> 0
 661   return -1 if $xz && !$yz;                             # 0 <=> +y
 662   return 1 if $yz && !$xz;                              # +x <=> 0
 663
 664   # adjust so that exponents are equal
 665   my $lxm = $MBI->_len($x->{_m});
 666   my $lym = $MBI->_len($y->{_m});
 667   my ($xes,$yes) = (1,1);
 668   $xes = -1 if $x->{_es} ne '+';
 669   $yes = -1 if $y->{_es} ne '+';
 670   # the numify somewhat limits our length, but makes it much faster
 671   my $lx = $lxm + $xes * $MBI->_num($x->{_e});
 672   my $ly = $lym + $yes * $MBI->_num($y->{_e});
 673   my $l = $lx - $ly;
 674   return $l <=> 0 if $l != 0;
 675
 676   # lengths (corrected by exponent) are equal
 677   # so make mantissa equal-length by padding with zero (shift left)
 678   my $diff = $lxm - $lym;
 679   my $xm = $x->{_m};            # not yet copy it
 680   my $ym = $y->{_m};
 681   if ($diff > 0)
 682     {
 683     $ym = $MBI->_copy($y->{_m});
 684     $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
 685     }
 686   elsif ($diff < 0)
 687     {
 688     $xm = $MBI->_copy($x->{_m});
 689     $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
 690     }
 691   $MBI->_acmp($xm,$ym);
 692   }
 693
 694 sub badd
 695   {
 696   # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
 697   # return result as BFLOAT
 698
 699   # set up parameters
 700   my ($self,$x,$y,@r) = (ref($_[0]),@_);
 701   # objectify is costly, so avoid it
 702   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 703     {
 704     ($self,$x,$y,@r) = objectify(2,@_);
 705     }
 706
 707   return $x if $x->modify('badd');
 708
 709   # inf and NaN handling
 710   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
 711     {
 712     # NaN first
 713     return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 714     # inf handling
 715     if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
 716       {
 717       # +inf++inf or -inf+-inf => same, rest is NaN
 718       return $x if $x->{sign} eq $y->{sign};
 719       return $x->bnan();
 720       }
 721     # +-inf + something => +inf; something +-inf => +-inf
 722     $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
 723     return $x;
 724     }
 725
 726   return $upgrade->badd($x,$y,@r) if defined $upgrade &&
 727    ((!$x->isa($self)) || (!$y->isa($self)));
 728
 729   $r[3] = $y;                                           # no push!
 730
 731   # speed: no add for 0+y or x+0
 732   return $x->bround(@r) if $y->is_zero();               # x+0
 733   if ($x->is_zero())                                    # 0+y
 734     {
 735     # make copy, clobbering up x (modify in place!)
 736     $x->{_e} = $MBI->_copy($y->{_e});
 737     $x->{_es} = $y->{_es};
 738     $x->{_m} = $MBI->_copy($y->{_m});
 739     $x->{sign} = $y->{sign} || $nan;
 740     return $x->round(@r);
 741     }
 742
 743   # take lower of the two e's and adapt m1 to it to match m2
 744   my $e = $y->{_e};
 745   $e = $MBI->_zero() if !defined $e;            # if no BFLOAT?
 746   $e = $MBI->_copy($e);                         # make copy (didn't do it yet)
 747
 748   my $es;
 749
 750   ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
 751
 752   my $add = $MBI->_copy($y->{_m});
 753
 754   if ($es eq '-')                               # < 0
 755     {
 756     $MBI->_lsft( $x->{_m}, $e, 10);
 757     ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
 758     }
 759   elsif (!$MBI->_is_zero($e))                   # > 0
 760     {
 761     $MBI->_lsft($add, $e, 10);
 762     }
 763   # else: both e are the same, so just leave them
 764
 765   if ($x->{sign} eq $y->{sign})
 766     {
 767     # add
 768     $x->{_m} = $MBI->_add($x->{_m}, $add);
 769     }
 770   else
 771     {
 772     ($x->{_m}, $x->{sign}) =
 773      _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
 774     }
 775
 776   # delete trailing zeros, then round
 777   $x->bnorm()->round(@r);
 778   }
 779
 780 # sub bsub is inherited from Math::BigInt!
 781
 782 sub binc
 783   {
 784   # increment arg by one
 785   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 786
 787   return $x if $x->modify('binc');
 788
 789   if ($x->{_es} eq '-')
 790     {
 791     return $x->badd($self->bone(),@r);  #  digits after dot
 792     }
 793
 794   if (!$MBI->_is_zero($x->{_e}))                # _e == 0 for NaN, inf, -inf
 795     {
 796     # 1e2 => 100, so after the shift below _m has a '0' as last digit
 797     $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);      # 1e2 => 100
 798     $x->{_e} = $MBI->_zero();                           # normalize
 799     $x->{_es} = '+';
 800     # we know that the last digit of $x will be '1' or '9', depending on the
 801     # sign
 802     }
 803   # now $x->{_e} == 0
 804   if ($x->{sign} eq '+')
 805     {
 806     $MBI->_inc($x->{_m});
 807     return $x->bnorm()->bround(@r);
 808     }
 809   elsif ($x->{sign} eq '-')
 810     {
 811     $MBI->_dec($x->{_m});
 812     $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
 813     return $x->bnorm()->bround(@r);
 814     }
 815   # inf, nan handling etc
 816   $x->badd($self->bone(),@r);                   # badd() does round
 817   }
 818
 819 sub bdec
 820   {
 821   # decrement arg by one
 822   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 823
 824   return $x if $x->modify('bdec');
 825
 826   if ($x->{_es} eq '-')
 827     {
 828     return $x->badd($self->bone('-'),@r);       #  digits after dot
 829     }
 830
 831   if (!$MBI->_is_zero($x->{_e}))
 832     {
 833     $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);      # 1e2 => 100
 834     $x->{_e} = $MBI->_zero();                           # normalize
 835     $x->{_es} = '+';
 836     }
 837   # now $x->{_e} == 0
 838   my $zero = $x->is_zero();
 839   # <= 0
 840   if (($x->{sign} eq '-') || $zero)
 841     {
 842     $MBI->_inc($x->{_m});
 843     $x->{sign} = '-' if $zero;                          # 0 => 1 => -1
 844     $x->{sign} = '+' if $MBI->_is_zero($x->{_m});       # -1 +1 => -0 => +0
 845     return $x->bnorm()->round(@r);
 846     }
 847   # > 0
 848   elsif ($x->{sign} eq '+')
 849     {
 850     $MBI->_dec($x->{_m});
 851     return $x->bnorm()->round(@r);
 852     }
 853   # inf, nan handling etc
 854   $x->badd($self->bone('-'),@r);                # does round
 855   }
 856
 857 sub DEBUG () { 0; }
 858
 859 sub blog
 860   {
 861   my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 862
 863   return $x if $x->modify('blog');
 864
 865   # $base > 0, $base != 1; if $base == undef default to $base == e
 866   # $x >= 0
 867
 868   # we need to limit the accuracy to protect against overflow
 869   my $fallback = 0;
 870   my ($scale,@params);
 871   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
 872
 873   # also takes care of the "error in _find_round_parameters?" case
 874   return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
 875
 876   # no rounding at all, so must use fallback
 877   if (scalar @params == 0)
 878     {
 879     # simulate old behaviour
 880     $params[0] = $self->div_scale();    # and round to it as accuracy
 881     $params[1] = undef;                 # P = undef
 882     $scale = $params[0]+4;              # at least four more for proper round
 883     $params[2] = $r;                    # round mode by caller or undef
 884     $fallback = 1;                      # to clear a/p afterwards
 885     }
 886   else
 887     {
 888     # the 4 below is empirical, and there might be cases where it is not
 889     # enough...
 890     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
 891     }
 892
 893   return $x->bzero(@params) if $x->is_one();
 894   # base not defined => base == Euler's number e
 895   if (defined $base)
 896     {
 897     # make object, since we don't feed it through objectify() to still get the
 898     # case of $base == undef
 899     $base = $self->new($base) unless ref($base);
 900     # $base > 0; $base != 1
 901     return $x->bnan() if $base->is_zero() || $base->is_one() ||
 902       $base->{sign} ne '+';
 903     # if $x == $base, we know the result must be 1.0
 904     if ($x->bcmp($base) == 0)
 905       {
 906       $x->bone('+',@params);
 907       if ($fallback)
 908         {
 909         # clear a/p after round, since user did not request it
 910         delete $x->{_a}; delete $x->{_p};
 911         }
 912       return $x;
 913       }
 914     }
 915
 916   # when user set globals, they would interfere with our calculation, so
 917   # disable them and later re-enable them
 918   no strict 'refs';
 919   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
 920   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
 921   # we also need to disable any set A or P on $x (_find_round_parameters took
 922   # them already into account), since these would interfere, too
 923   delete $x->{_a}; delete $x->{_p};
 924   # need to disable $upgrade in BigInt, to avoid deep recursion
 925   local $Math::BigInt::upgrade = undef;
 926   local $Math::BigFloat::downgrade = undef;
 927
 928   # upgrade $x if $x is not a BigFloat (handle BigInt input)
 929   # XXX TODO: rebless!
 930   if (!$x->isa('Math::BigFloat'))
 931     {
 932     $x = Math::BigFloat->new($x);
 933     $self = ref($x);
 934     }
 935
 936   my $done = 0;
 937
 938   # If the base is defined and an integer, try to calculate integer result
 939   # first. This is very fast, and in case the real result was found, we can
 940   # stop right here.
 941   if (defined $base && $base->is_int() && $x->is_int())
 942     {
 943     my $i = $MBI->_copy( $x->{_m} );
 944     $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
 945     my $int = Math::BigInt->bzero();
 946     $int->{value} = $i;
 947     $int->blog($base->as_number());
 948     # if ($exact)
 949     if ($base->as_number()->bpow($int) == $x)
 950       {
 951       # found result, return it
 952       $x->{_m} = $int->{value};
 953       $x->{_e} = $MBI->_zero();
 954       $x->{_es} = '+';
 955       $x->bnorm();
 956       $done = 1;
 957       }
 958     }
 959
 960   if ($done == 0)
 961     {
 962     # base is undef, so base should be e (Euler's number), so first calculate the
 963     # log to base e (using reduction by 10 (and probably 2)):
 964     $self->_log_10($x,$scale);
 965
 966     # and if a different base was requested, convert it
 967     if (defined $base)
 968       {
 969       $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
 970       # not ln, but some other base (don't modify $base)
 971       $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
 972       }
 973     }
 974
 975   # shortcut to not run through _find_round_parameters again
 976   if (defined $params[0])
 977     {
 978     $x->bround($params[0],$params[2]);          # then round accordingly
 979     }
 980   else
 981     {
 982     $x->bfround($params[1],$params[2]);         # then round accordingly
 983     }
 984   if ($fallback)
 985     {
 986     # clear a/p after round, since user did not request it
 987     delete $x->{_a}; delete $x->{_p};
 988     }
 989   # restore globals
 990   $$abr = $ab; $$pbr = $pb;
 991
 992   $x;
 993   }
 994
 995 sub _len_to_steps
 996   {
 997   # Given D (digits in decimal), compute N so that N! (N factorial) is
 998   # at least D digits long. D should be at least 50.
 999   my $d = shift;
1000
1001   # two constants for the Ramanujan estimate of ln(N!)
1002   my $lg2 = log(2 * 3.14159265) / 2;
1003   my $lg10 = log(10);
1004
1005   # D = 50 => N => 42, so L = 40 and R = 50
1006   my $l = 40; my $r = $d;
1007
1008   # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
1009   $l = $l->numify if ref($l);
1010   $r = $r->numify if ref($r);
1011   $lg2 = $lg2->numify if ref($lg2);
1012   $lg10 = $lg10->numify if ref($lg10);
1013
1014   # binary search for the right value (could this be written as the reverse of lg(n!)?)
1015   while ($r - $l > 1)
1016     {
1017     my $n = int(($r - $l) / 2) + $l;
1018     my $ramanujan =
1019       int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
1020     $ramanujan > $d ? $r = $n : $l = $n;
1021     }
1022   $l;
1023   }
1024
1025 sub bnok
1026   {
1027   # Calculate n over k (binomial coefficient or "choose" function) as integer.
1028   # set up parameters
1029   my ($self,$x,$y,@r) = (ref($_[0]),@_);
1030
1031   # objectify is costly, so avoid it
1032   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1033     {
1034     ($self,$x,$y,@r) = objectify(2,@_);
1035     }
1036
1037   return $x if $x->modify('bnok');
1038
1039   return $x->bnan() if $x->is_nan() || $y->is_nan();
1040   return $x->binf() if $x->is_inf();
1041
1042   my $u = $x->as_int();
1043   $u->bnok($y->as_int());
1044
1045   $x->{_m} = $u->{value};
1046   $x->{_e} = $MBI->_zero();
1047   $x->{_es} = '+';
1048   $x->{sign} = '+';
1049   $x->bnorm(@r);
1050   }
1051
1052 sub bexp
1053   {
1054   # Calculate e ** X (Euler's number to the power of X)
1055   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1056
1057   return $x if $x->modify('bexp');
1058
1059   return $x->binf() if $x->{sign} eq '+inf';
1060   return $x->bzero() if $x->{sign} eq '-inf';
1061
1062   # we need to limit the accuracy to protect against overflow
1063   my $fallback = 0;
1064   my ($scale,@params);
1065   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1066
1067   # also takes care of the "error in _find_round_parameters?" case
1068   return $x if $x->{sign} eq 'NaN';
1069
1070   # no rounding at all, so must use fallback
1071   if (scalar @params == 0)
1072     {
1073     # simulate old behaviour
1074     $params[0] = $self->div_scale();    # and round to it as accuracy
1075     $params[1] = undef;                 # P = undef
1076     $scale = $params[0]+4;              # at least four more for proper round
1077     $params[2] = $r;                    # round mode by caller or undef
1078     $fallback = 1;                      # to clear a/p afterwards
1079     }
1080   else
1081     {
1082     # the 4 below is empirical, and there might be cases where it's not enough...
1083     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1084     }
1085
1086   return $x->bone(@params) if $x->is_zero();
1087
1088   if (!$x->isa('Math::BigFloat'))
1089     {
1090     $x = Math::BigFloat->new($x);
1091     $self = ref($x);
1092     }
1093
1094   # when user set globals, they would interfere with our calculation, so
1095   # disable them and later re-enable them
1096   no strict 'refs';
1097   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1098   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1099   # we also need to disable any set A or P on $x (_find_round_parameters took
1100   # them already into account), since these would interfere, too
1101   delete $x->{_a}; delete $x->{_p};
1102   # need to disable $upgrade in BigInt, to avoid deep recursion
1103   local $Math::BigInt::upgrade = undef;
1104   local $Math::BigFloat::downgrade = undef;
1105
1106   my $x_org = $x->copy();
1107
1108   # We use the following Taylor series:
1109
1110   #           x    x^2   x^3   x^4
1111   #  e = 1 + --- + --- + --- + --- ...
1112   #           1!    2!    3!    4!
1113
1114   # The difference for each term is X and N, which would result in:
1115   # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term
1116
1117   # But it is faster to compute exp(1) and then raising it to the
1118   # given power, esp. if $x is really big and an integer because:
1119
1120   #  * The numerator is always 1, making the computation faster
1121   #  * the series converges faster in the case of x == 1
1122   #  * We can also easily check when we have reached our limit: when the
1123   #    term to be added is smaller than "1E$scale", we can stop - f.i.
1124   #    scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
1125   #  * we can compute the *exact* result by simulating bigrat math:
1126
1127   #  1   1    gcd(3,4) = 1    1*24 + 1*6    5
1128   #  - + -                  = ---------- =  --
1129   #  6   24                      6*24       24
1130
1131   # We do not compute the gcd() here, but simple do:
1132   #  1   1    1*24 + 1*6   30
1133   #  - + -  = --------- =  --
1134   #  6   24       6*24     144
1135
1136   # In general:
1137   #  a   c    a*d + c*b         and note that c is always 1 and d = (b*f)
1138   #  - + -  = ---------
1139   #  b   d       b*d
1140
1141   # This leads to:         which can be reduced by b to:
1142   #  a   1     a*b*f + b    a*f + 1
1143   #  - + -   = --------- =  -------
1144   #  b   b*f     b*b*f        b*f
1145
1146   # The first terms in the series are:
1147
1148   # 1     1    1    1    1    1     1     1     13700
1149   # -- + -- + -- + -- + -- + --- + --- + ---- = -----
1150   # 1     1    2    6   24   120   720   5040   5040
1151
1152   # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!
1153
1154   if ($scale <= 75)
1155     {
1156     # set $x directly from a cached string form
1157     $x->{_m} = $MBI->_new(
1158     "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
1159     $x->{sign} = '+';
1160     $x->{_es} = '-';
1161     $x->{_e} = $MBI->_new(79);
1162     }
1163   else
1164     {
1165     # compute A and B so that e = A / B.
1166
1167     # After some terms we end up with this, so we use it as a starting point:
1168     my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
1169     my $F = $MBI->_new(42); my $step = 42;
1170
1171     # Compute how many steps we need to take to get $A and $B sufficiently big
1172     my $steps = _len_to_steps($scale - 4);
1173 #    print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
1174     while ($step++ <= $steps)
1175       {
1176       # calculate $a * $f + 1
1177       $A = $MBI->_mul($A, $F);
1178       $A = $MBI->_inc($A);
1179       # increment f
1180       $F = $MBI->_inc($F);
1181       }
1182     # compute $B as factorial of $steps (this is faster than doing it manually)
1183     my $B = $MBI->_fac($MBI->_new($steps));
1184
1185 #  print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
1186
1187     # compute A/B with $scale digits in the result (truncate, not round)
1188     $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
1189     $A = $MBI->_div( $A, $B );
1190
1191     $x->{_m} = $A;
1192     $x->{sign} = '+';
1193     $x->{_es} = '-';
1194     $x->{_e} = $MBI->_new($scale);
1195     }
1196
1197   # $x contains now an estimate of e, with some surplus digits, so we can round
1198   if (!$x_org->is_one())
1199     {
1200     # raise $x to the wanted power and round it in one step:
1201     $x->bpow($x_org, @params);
1202     }
1203   else
1204     {
1205     # else just round the already computed result
1206     delete $x->{_a}; delete $x->{_p};
1207     # shortcut to not run through _find_round_parameters again
1208     if (defined $params[0])
1209       {
1210       $x->bround($params[0],$params[2]);                # then round accordingly
1211       }
1212     else
1213       {
1214       $x->bfround($params[1],$params[2]);               # then round accordingly
1215       }
1216     }
1217   if ($fallback)
1218     {
1219     # clear a/p after round, since user did not request it
1220     delete $x->{_a}; delete $x->{_p};
1221     }
1222   # restore globals
1223   $$abr = $ab; $$pbr = $pb;
1224
1225   $x;                                           # return modified $x
1226   }
1227
1228 sub _log
1229   {
1230   # internal log function to calculate ln() based on Taylor series.
1231   # Modifies $x in place.
1232   my ($self,$x,$scale) = @_;
1233
1234   # in case of $x == 1, result is 0
1235   return $x->bzero() if $x->is_one();
1236
1237   # XXX TODO: rewrite this in a similar manner to bexp()
1238
1239   # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
1240
1241   # u = x-1, v = x+1
1242   #              _                               _
1243   # Taylor:     |    u    1   u^3   1   u^5       |
1244   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
1245   #             |_   v    3   v^3   5   v^5      _|
1246
1247   # This takes much more steps to calculate the result and is thus not used
1248   # u = x-1
1249   #              _                               _
1250   # Taylor:     |    u    1   u^2   1   u^3       |
1251   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 1/2
1252   #             |_   x    2   x^2   3   x^3      _|
1253
1254   my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
1255
1256   $v = $x->copy(); $v->binc();          # v = x+1
1257   $x->bdec(); $u = $x->copy();          # u = x-1; x = x-1
1258   $x->bdiv($v,$scale);                  # first term: u/v
1259   $below = $v->copy();
1260   $over = $u->copy();
1261   $u *= $u; $v *= $v;                           # u^2, v^2
1262   $below->bmul($v);                             # u^3, v^3
1263   $over->bmul($u);
1264   $factor = $self->new(3); $f = $self->new(2);
1265
1266   my $steps = 0 if DEBUG;
1267   $limit = $self->new("1E-". ($scale-1));
1268   while (3 < 5)
1269     {
1270     # we calculate the next term, and add it to the last
1271     # when the next term is below our limit, it won't affect the outcome
1272     # anymore, so we stop
1273
1274     # calculating the next term simple from over/below will result in quite
1275     # a time hog if the input has many digits, since over and below will
1276     # accumulate more and more digits, and the result will also have many
1277     # digits, but in the end it is rounded to $scale digits anyway. So if we
1278     # round $over and $below first, we save a lot of time for the division
1279     # (not with log(1.2345), but try log (123**123) to see what I mean. This
1280     # can introduce a rounding error if the division result would be f.i.
1281     # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
1282     # if we truncated $over and $below we might get 0.12345. Does this matter
1283     # for the end result? So we give $over and $below 4 more digits to be
1284     # on the safe side (unscientific error handling as usual... :+D
1285
1286     $next = $over->copy->bround($scale+4)->bdiv(
1287       $below->copy->bmul($factor)->bround($scale+4),
1288       $scale);
1289
1290 ## old version:
1291 ##    $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
1292
1293     last if $next->bacmp($limit) <= 0;
1294
1295     delete $next->{_a}; delete $next->{_p};
1296     $x->badd($next);
1297     # calculate things for the next term
1298     $over *= $u; $below *= $v; $factor->badd($f);
1299     if (DEBUG)
1300       {
1301       $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
1302       }
1303     }
1304   print "took $steps steps\n" if DEBUG;
1305   $x->bmul($f);                                 # $x *= 2
1306   }
1307
1308 sub _log_10
1309   {
1310   # Internal log function based on reducing input to the range of 0.1 .. 9.99
1311   # and then "correcting" the result to the proper one. Modifies $x in place.
1312   my ($self,$x,$scale) = @_;
1313
1314   # Taking blog() from numbers greater than 10 takes a *very long* time, so we
1315   # break the computation down into parts based on the observation that:
1316   #  blog(X*Y) = blog(X) + blog(Y)
1317   # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
1318   # $x is the faster it gets. Since 2*$x takes about 10 times as
1319   # long, we make it faster by about a factor of 100 by dividing $x by 10.
1320
1321   # The same observation is valid for numbers smaller than 0.1, e.g. computing
1322   # log(1) is fastest, and the further away we get from 1, the longer it takes.
1323   # So we also 'break' this down by multiplying $x with 10 and subtract the
1324   # log(10) afterwards to get the correct result.
1325
1326   # To get $x even closer to 1, we also divide by 2 and then use log(2) to
1327   # correct for this. For instance if $x is 2.4, we use the formula:
1328   #  blog(2.4 * 2) == blog (1.2) + blog(2)
1329   # and thus calculate only blog(1.2) and blog(2), which is faster in total
1330   # than calculating blog(2.4).
1331
1332   # In addition, the values for blog(2) and blog(10) are cached.
1333
1334   # Calculate nr of digits before dot:
1335   my $dbd = $MBI->_num($x->{_e});
1336   $dbd = -$dbd if $x->{_es} eq '-';
1337   $dbd += $MBI->_len($x->{_m});
1338
1339   # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
1340   # infinite recursion
1341
1342   my $calc = 1;                                 # do some calculation?
1343
1344   # disable the shortcut for 10, since we need log(10) and this would recurse
1345   # infinitely deep
1346   if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
1347     {
1348     $dbd = 0;                                   # disable shortcut
1349     # we can use the cached value in these cases
1350     if ($scale <= $LOG_10_A)
1351       {
1352       $x->bzero(); $x->badd($LOG_10);           # modify $x in place
1353       $calc = 0;                                # no need to calc, but round
1354       }
1355     # if we can't use the shortcut, we continue normally
1356     }
1357   else
1358     {
1359     # disable the shortcut for 2, since we maybe have it cached
1360     if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
1361       {
1362       $dbd = 0;                                 # disable shortcut
1363       # we can use the cached value in these cases
1364       if ($scale <= $LOG_2_A)
1365         {
1366         $x->bzero(); $x->badd($LOG_2);          # modify $x in place
1367         $calc = 0;                              # no need to calc, but round
1368         }
1369       # if we can't use the shortcut, we continue normally
1370       }
1371     }
1372
1373   # if $x = 0.1, we know the result must be 0-log(10)
1374   if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
1375       $MBI->_is_one($x->{_m}))
1376     {
1377     $dbd = 0;                                   # disable shortcut
1378     # we can use the cached value in these cases
1379     if ($scale <= $LOG_10_A)
1380       {
1381       $x->bzero(); $x->bsub($LOG_10);
1382       $calc = 0;                                # no need to calc, but round
1383       }
1384     }
1385
1386   return if $calc == 0;                         # already have the result
1387
1388   # default: these correction factors are undef and thus not used
1389   my $l_10;                             # value of ln(10) to A of $scale
1390   my $l_2;                              # value of ln(2) to A of $scale
1391
1392   my $two = $self->new(2);
1393
1394   # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1395   # so don't do this shortcut for 1 or 0
1396   if (($dbd > 1) || ($dbd < 0))
1397     {
1398     # convert our cached value to an object if not already (avoid doing this
1399     # at import() time, since not everybody needs this)
1400     $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1401
1402     #print "x = $x, dbd = $dbd, calc = $calc\n";
1403     # got more than one digit before the dot, or more than one zero after the
1404     # dot, so do:
1405     #  log(123)    == log(1.23) + log(10) * 2
1406     #  log(0.0123) == log(1.23) - log(10) * 2
1407
1408     if ($scale <= $LOG_10_A)
1409       {
1410       # use cached value
1411       $l_10 = $LOG_10->copy();          # copy for mul
1412       }
1413     else
1414       {
1415       # else: slower, compute and cache result
1416       # also disable downgrade for this code path
1417       local $Math::BigFloat::downgrade = undef;
1418
1419       # shorten the time to calculate log(10) based on the following:
1420       # log(1.25 * 8) = log(1.25) + log(8)
1421       #               = log(1.25) + log(2) + log(2) + log(2)
1422
1423       # first get $l_2 (and possible compute and cache log(2))
1424       $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1425       if ($scale <= $LOG_2_A)
1426         {
1427         # use cached value
1428         $l_2 = $LOG_2->copy();                  # copy() for the mul below
1429         }
1430       else
1431         {
1432         # else: slower, compute and cache result
1433         $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1434         $LOG_2 = $l_2->copy();                  # cache the result for later
1435                                                 # the copy() is for mul below
1436         $LOG_2_A = $scale;
1437         }
1438
1439       # now calculate log(1.25):
1440       $l_10 = $self->new('1.25'); $self->_log($l_10, $scale); # scale+4, actually
1441
1442       # log(1.25) + log(2) + log(2) + log(2):
1443       $l_10->badd($l_2);
1444       $l_10->badd($l_2);
1445       $l_10->badd($l_2);
1446       $LOG_10 = $l_10->copy();          # cache the result for later
1447                                         # the copy() is for mul below
1448       $LOG_10_A = $scale;
1449       }
1450     $dbd-- if ($dbd > 1);               # 20 => dbd=2, so make it dbd=1
1451     $l_10->bmul( $self->new($dbd));     # log(10) * (digits_before_dot-1)
1452     my $dbd_sign = '+';
1453     if ($dbd < 0)
1454       {
1455       $dbd = -$dbd;
1456       $dbd_sign = '-';
1457       }
1458     ($x->{_e}, $x->{_es}) =
1459         _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
1460
1461     }
1462
1463   # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1464
1465   ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1466   ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1467
1468   $HALF = $self->new($HALF) unless ref($HALF);
1469
1470   my $twos = 0;                         # default: none (0 times)
1471   while ($x->bacmp($HALF) <= 0)         # X <= 0.5
1472     {
1473     $twos--; $x->bmul($two);
1474     }
1475   while ($x->bacmp($two) >= 0)          # X >= 2
1476     {
1477     $twos++; $x->bdiv($two,$scale+4);           # keep all digits
1478     }
1479   # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
1480   # So calculate correction factor based on ln(2):
1481   if ($twos != 0)
1482     {
1483     $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1484     if ($scale <= $LOG_2_A)
1485       {
1486       # use cached value
1487       $l_2 = $LOG_2->copy();                    # copy() for the mul below
1488       }
1489     else
1490       {
1491       # else: slower, compute and cache result
1492       # also disable downgrade for this code path
1493       local $Math::BigFloat::downgrade = undef;
1494       $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1495       $LOG_2 = $l_2->copy();                    # cache the result for later
1496                                                 # the copy() is for mul below
1497       $LOG_2_A = $scale;
1498       }
1499     $l_2->bmul($twos);          # * -2 => subtract, * 2 => add
1500     }
1501
1502   $self->_log($x,$scale);                       # need to do the "normal" way
1503   $x->badd($l_10) if defined $l_10;             # correct it by ln(10)
1504   $x->badd($l_2) if defined $l_2;               # and maybe by ln(2)
1505
1506   # all done, $x contains now the result
1507   $x;
1508   }
1509
1510 sub blcm
1511   {
1512   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1513   # does not modify arguments, but returns new object
1514   # Lowest Common Multiplicator
1515
1516   my ($self,@arg) = objectify(0,@_);
1517   my $x = $self->new(shift @arg);
1518   while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); }
1519   $x;
1520   }
1521
1522 sub bgcd
1523   {
1524   # (BINT or num_str, BINT or num_str) return BINT
1525   # does not modify arguments, but returns new object
1526
1527   my $y = shift;
1528   $y = __PACKAGE__->new($y) if !ref($y);
1529   my $self = ref($y);
1530   my $x = $y->copy()->babs();                   # keep arguments
1531
1532   return $x->bnan() if $x->{sign} !~ /^[+-]$/   # x NaN?
1533         || !$x->is_int();                       # only for integers now
1534
1535   while (@_)
1536     {
1537     my $t = shift; $t = $self->new($t) if !ref($t);
1538     $y = $t->copy()->babs();
1539
1540     return $x->bnan() if $y->{sign} !~ /^[+-]$/ # y NaN?
1541         || !$y->is_int();                       # only for integers now
1542
1543     # greatest common divisor
1544     while (! $y->is_zero())
1545       {
1546       ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1547       }
1548
1549     last if $x->is_one();
1550     }
1551   $x;
1552   }
1553
1554 ##############################################################################
1555
1556 sub _e_add
1557   {
1558   # Internal helper sub to take two positive integers and their signs and
1559   # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1560   # output ($CALC,('+'|'-'))
1561   my ($x,$y,$xs,$ys) = @_;
1562
1563   # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1564   if ($xs eq $ys)
1565     {
1566     $x = $MBI->_add ($x, $y );          # a+b
1567     # the sign follows $xs
1568     return ($x, $xs);
1569     }
1570
1571   my $a = $MBI->_acmp($x,$y);
1572   if ($a > 0)
1573     {
1574     $x = $MBI->_sub ($x , $y);                          # abs sub
1575     }
1576   elsif ($a == 0)
1577     {
1578     $x = $MBI->_zero();                                 # result is 0
1579     $xs = '+';
1580     }
1581   else # a < 0
1582     {
1583     $x = $MBI->_sub ( $y, $x, 1 );                      # abs sub
1584     $xs = $ys;
1585     }
1586   ($x,$xs);
1587   }
1588
1589 sub _e_sub
1590   {
1591   # Internal helper sub to take two positive integers and their signs and
1592   # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1593   # output ($CALC,('+'|'-'))
1594   my ($x,$y,$xs,$ys) = @_;
1595
1596   # flip sign
1597   $ys =~ tr/+-/-+/;
1598   _e_add($x,$y,$xs,$ys);                # call add (does subtract now)
1599   }
1600
1601 ###############################################################################
1602 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1603
1604 sub is_int
1605   {
1606   # return true if arg (BFLOAT or num_str) is an integer
1607   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1608
1609   (($x->{sign} =~ /^[+-]$/) &&                  # NaN and +-inf aren't
1610    ($x->{_es} eq '+')) ? 1 : 0;                 # 1e-1 => no integer
1611   }
1612
1613 sub is_zero
1614   {
1615   # return true if arg (BFLOAT or num_str) is zero
1616   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1617
1618   ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m})) ? 1 : 0;
1619   }
1620
1621 sub is_one
1622   {
1623   # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1624   my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1625
1626   $sign = '+' if !defined $sign || $sign ne '-';
1627
1628   ($x->{sign} eq $sign &&
1629    $MBI->_is_zero($x->{_e}) &&
1630    $MBI->_is_one($x->{_m}) ) ? 1 : 0;
1631   }
1632
1633 sub is_odd
1634   {
1635   # return true if arg (BFLOAT or num_str) is odd or false if even
1636   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1637
1638   (($x->{sign} =~ /^[+-]$/) &&          # NaN & +-inf aren't
1639    ($MBI->_is_zero($x->{_e})) &&
1640    ($MBI->_is_odd($x->{_m}))) ? 1 : 0;
1641   }
1642
1643 sub is_even
1644   {
1645   # return true if arg (BINT or num_str) is even or false if odd
1646   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1647
1648   (($x->{sign} =~ /^[+-]$/) &&                  # NaN & +-inf aren't
1649    ($x->{_es} eq '+') &&                        # 123.45 isn't
1650    ($MBI->_is_even($x->{_m}))) ? 1 : 0;         # but 1200 is
1651   }
1652
1653 sub bmul
1654   {
1655   # multiply two numbers
1656
1657   # set up parameters
1658   my ($self,$x,$y,@r) = (ref($_[0]),@_);
1659   # objectify is costly, so avoid it
1660   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1661     {
1662     ($self,$x,$y,@r) = objectify(2,@_);
1663     }
1664
1665   return $x if $x->modify('bmul');
1666
1667   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1668
1669   # inf handling
1670   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1671     {
1672     return $x->bnan() if $x->is_zero() || $y->is_zero();
1673     # result will always be +-inf:
1674     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1675     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1676     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1677     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1678     return $x->binf('-');
1679     }
1680
1681   return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1682    ((!$x->isa($self)) || (!$y->isa($self)));
1683
1684   # aEb * cEd = (a*c)E(b+d)
1685   $MBI->_mul($x->{_m},$y->{_m});
1686   ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1687
1688   $r[3] = $y;                           # no push!
1689
1690   # adjust sign:
1691   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1692   $x->bnorm->round(@r);
1693   }
1694
1695 sub bmuladd
1696   {
1697   # multiply two numbers and add the third to the result
1698
1699   # set up parameters
1700   my ($self,$x,$y,$z,@r) = (ref($_[0]),@_);
1701   # objectify is costly, so avoid it
1702   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1703     {
1704     ($self,$x,$y,$z,@r) = objectify(3,@_);
1705     }
1706
1707   return $x if $x->modify('bmuladd');
1708
1709   return $x->bnan() if (($x->{sign} eq $nan) ||
1710                         ($y->{sign} eq $nan) ||
1711                         ($z->{sign} eq $nan));
1712
1713   # inf handling
1714   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1715     {
1716     return $x->bnan() if $x->is_zero() || $y->is_zero();
1717     # result will always be +-inf:
1718     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1719     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1720     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1721     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1722     return $x->binf('-');
1723     }
1724
1725   return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1726    ((!$x->isa($self)) || (!$y->isa($self)));
1727
1728   # aEb * cEd = (a*c)E(b+d)
1729   $MBI->_mul($x->{_m},$y->{_m});
1730   ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1731
1732   $r[3] = $y;                           # no push!
1733
1734   # adjust sign:
1735   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1736
1737   # z=inf handling (z=NaN handled above)
1738   $x->{sign} = $z->{sign}, return $x if $z->{sign} =~ /^[+-]inf$/;
1739
1740   # take lower of the two e's and adapt m1 to it to match m2
1741   my $e = $z->{_e};
1742   $e = $MBI->_zero() if !defined $e;            # if no BFLOAT?
1743   $e = $MBI->_copy($e);                         # make copy (didn't do it yet)
1744
1745   my $es;
1746
1747   ($e,$es) = _e_sub($e, $x->{_e}, $z->{_es} || '+', $x->{_es});
1748
1749   my $add = $MBI->_copy($z->{_m});
1750
1751   if ($es eq '-')                               # < 0
1752     {
1753     $MBI->_lsft( $x->{_m}, $e, 10);
1754     ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
1755     }
1756   elsif (!$MBI->_is_zero($e))                   # > 0
1757     {
1758     $MBI->_lsft($add, $e, 10);
1759     }
1760   # else: both e are the same, so just leave them
1761
1762   if ($x->{sign} eq $z->{sign})
1763     {
1764     # add
1765     $x->{_m} = $MBI->_add($x->{_m}, $add);
1766     }
1767   else
1768     {
1769     ($x->{_m}, $x->{sign}) =
1770      _e_add($x->{_m}, $add, $x->{sign}, $z->{sign});
1771     }
1772
1773   # delete trailing zeros, then round
1774   $x->bnorm()->round(@r);
1775   }
1776
1777 sub bdiv
1778   {
1779   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1780   # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1781
1782   # set up parameters
1783   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1784   # objectify is costly, so avoid it
1785   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1786     {
1787     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1788     }
1789
1790   return $x if $x->modify('bdiv');
1791
1792   return $self->_div_inf($x,$y)
1793    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1794
1795   # x== 0 # also: or y == 1 or y == -1
1796   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1797
1798   # upgrade ?
1799   return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1800
1801   # we need to limit the accuracy to protect against overflow
1802   my $fallback = 0;
1803   my (@params,$scale);
1804   ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1805
1806   return $x if $x->is_nan();            # error in _find_round_parameters?
1807
1808   # no rounding at all, so must use fallback
1809   if (scalar @params == 0)
1810     {
1811     # simulate old behaviour
1812     $params[0] = $self->div_scale();    # and round to it as accuracy
1813     $scale = $params[0]+4;              # at least four more for proper round
1814     $params[2] = $r;                    # round mode by caller or undef
1815     $fallback = 1;                      # to clear a/p afterwards
1816     }
1817   else
1818     {
1819     # the 4 below is empirical, and there might be cases where it is not
1820     # enough...
1821     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1822     }
1823
1824   my $rem; $rem = $self->bzero() if wantarray;
1825
1826   $y = $self->new($y) unless $y->isa('Math::BigFloat');
1827
1828   my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
1829   $scale = $lx if $lx > $scale;
1830   $scale = $ly if $ly > $scale;
1831   my $diff = $ly - $lx;
1832   $scale += $diff if $diff > 0;         # if lx << ly, but not if ly << lx!
1833
1834   # already handled inf/NaN/-inf above:
1835
1836   # check that $y is not 1 nor -1 and cache the result:
1837   my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));
1838
1839   # flipping the sign of $y will also flip the sign of $x for the special
1840   # case of $x->bsub($x); so we can catch it below:
1841   my $xsign = $x->{sign};
1842   $y->{sign} =~ tr/+-/-+/;
1843
1844   if ($xsign ne $x->{sign})
1845     {
1846     # special case of $x /= $x results in 1
1847     $x->bone();                 # "fixes" also sign of $y, since $x is $y
1848     }
1849   else
1850     {
1851     # correct $y's sign again
1852     $y->{sign} =~ tr/+-/-+/;
1853     # continue with normal div code:
1854
1855     # make copy of $x in case of list context for later remainder calculation
1856     if (wantarray && $y_not_one)
1857       {
1858       $rem = $x->copy();
1859       }
1860
1861     $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1862
1863     # check for / +-1 ( +/- 1E0)
1864     if ($y_not_one)
1865       {
1866       # promote BigInts and it's subclasses (except when already a BigFloat)
1867       $y = $self->new($y) unless $y->isa('Math::BigFloat');
1868
1869       # calculate the result to $scale digits and then round it
1870       # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1871       $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1872       $MBI->_div ($x->{_m},$y->{_m});   # a/c
1873
1874       # correct exponent of $x
1875       ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1876       # correct for 10**scale
1877       ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1878       $x->bnorm();              # remove trailing 0's
1879       }
1880     } # ende else $x != $y
1881
1882   # shortcut to not run through _find_round_parameters again
1883   if (defined $params[0])
1884     {
1885     delete $x->{_a};                            # clear before round
1886     $x->bround($params[0],$params[2]);          # then round accordingly
1887     }
1888   else
1889     {
1890     delete $x->{_p};                            # clear before round
1891     $x->bfround($params[1],$params[2]);         # then round accordingly
1892     }
1893   if ($fallback)
1894     {
1895     # clear a/p after round, since user did not request it
1896     delete $x->{_a}; delete $x->{_p};
1897     }
1898
1899   if (wantarray)
1900     {
1901     if ($y_not_one)
1902       {
1903       $rem->bmod($y,@params);                   # copy already done
1904       }
1905     if ($fallback)
1906       {
1907       # clear a/p after round, since user did not request it
1908       delete $rem->{_a}; delete $rem->{_p};
1909       }
1910     return ($x,$rem);
1911     }
1912   $x;
1913   }
1914
1915 sub bmod
1916   {
1917   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return remainder
1918
1919   # set up parameters
1920   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1921   # objectify is costly, so avoid it
1922   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1923     {
1924     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1925     }
1926
1927   return $x if $x->modify('bmod');
1928
1929   # handle NaN, inf, -inf
1930   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1931     {
1932     my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1933     $x->{sign} = $re->{sign};
1934     $x->{_e} = $re->{_e};
1935     $x->{_m} = $re->{_m};
1936     return $x->round($a,$p,$r,$y);
1937     }
1938   if ($y->is_zero())
1939     {
1940     return $x->bnan() if $x->is_zero();
1941     return $x;
1942     }
1943
1944   return $x->bzero() if $x->is_zero()
1945  || ($x->is_int() &&
1946   # check that $y == +1 or $y == -1:
1947     ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));
1948
1949   my $cmp = $x->bacmp($y);                      # equal or $x < $y?
1950   return $x->bzero($a,$p) if $cmp == 0;         # $x == $y => result 0
1951
1952   # only $y of the operands negative?
1953   my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1954
1955   $x->{sign} = $y->{sign};                              # calc sign first
1956   return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0;  # $x < $y => result $x
1957
1958   my $ym = $MBI->_copy($y->{_m});
1959
1960   # 2e1 => 20
1961   $MBI->_lsft( $ym, $y->{_e}, 10)
1962    if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1963
1964   # if $y has digits after dot
1965   my $shifty = 0;                       # correct _e of $x by this
1966   if ($y->{_es} eq '-')                 # has digits after dot
1967     {
1968     # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1969     $shifty = $MBI->_num($y->{_e});     # no more digits after dot
1970     $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1971     }
1972   # $ym is now mantissa of $y based on exponent 0
1973
1974   my $shiftx = 0;                       # correct _e of $x by this
1975   if ($x->{_es} eq '-')                 # has digits after dot
1976     {
1977     # 123.4 % 20 => 1234 % 200
1978     $shiftx = $MBI->_num($x->{_e});     # no more digits after dot
1979     $MBI->_lsft($ym, $x->{_e}, 10);     # 123 => 1230
1980     }
1981   # 123e1 % 20 => 1230 % 20
1982   if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1983     {
1984     $MBI->_lsft( $x->{_m}, $x->{_e},10);        # es => '+' here
1985     }
1986
1987   $x->{_e} = $MBI->_new($shiftx);
1988   $x->{_es} = '+';
1989   $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1990   $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1991
1992   # now mantissas are equalized, exponent of $x is adjusted, so calc result
1993
1994   $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1995
1996   $x->{sign} = '+' if $MBI->_is_zero($x->{_m});         # fix sign for -0
1997   $x->bnorm();
1998
1999   if ($neg != 0)        # one of them negative => correct in place
2000     {
2001     my $r = $y - $x;
2002     $x->{_m} = $r->{_m};
2003     $x->{_e} = $r->{_e};
2004     $x->{_es} = $r->{_es};
2005     $x->{sign} = '+' if $MBI->_is_zero($x->{_m});       # fix sign for -0
2006     $x->bnorm();
2007     }
2008
2009   $x->round($a,$p,$r,$y);       # round and return
2010   }
2011
2012 sub broot
2013   {
2014   # calculate $y'th root of $x
2015
2016   # set up parameters
2017   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2018   # objectify is costly, so avoid it
2019   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2020     {
2021     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2022     }
2023
2024   return $x if $x->modify('broot');
2025
2026   # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
2027   return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
2028          $y->{sign} !~ /^\+$/;
2029
2030   return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
2031
2032   # we need to limit the accuracy to protect against overflow
2033   my $fallback = 0;
2034   my (@params,$scale);
2035   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2036
2037   return $x if $x->is_nan();            # error in _find_round_parameters?
2038
2039   # no rounding at all, so must use fallback
2040   if (scalar @params == 0)
2041     {
2042     # simulate old behaviour
2043     $params[0] = $self->div_scale();    # and round to it as accuracy
2044     $scale = $params[0]+4;              # at least four more for proper round
2045     $params[2] = $r;                    # iound mode by caller or undef
2046     $fallback = 1;                      # to clear a/p afterwards
2047     }
2048   else
2049     {
2050     # the 4 below is empirical, and there might be cases where it is not
2051     # enough...
2052     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2053     }
2054
2055   # when user set globals, they would interfere with our calculation, so
2056   # disable them and later re-enable them
2057   no strict 'refs';
2058   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2059   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2060   # we also need to disable any set A or P on $x (_find_round_parameters took
2061   # them already into account), since these would interfere, too
2062   delete $x->{_a}; delete $x->{_p};
2063   # need to disable $upgrade in BigInt, to avoid deep recursion
2064   local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2065
2066   # remember sign and make $x positive, since -4 ** (1/2) => -2
2067   my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
2068
2069   my $is_two = 0;
2070   if ($y->isa('Math::BigFloat'))
2071     {
2072     $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
2073     }
2074   else
2075     {
2076     $is_two = ($y == 2);
2077     }
2078
2079   # normal square root if $y == 2:
2080   if ($is_two)
2081     {
2082     $x->bsqrt($scale+4);
2083     }
2084   elsif ($y->is_one('-'))
2085     {
2086     # $x ** -1 => 1/$x
2087     my $u = $self->bone()->bdiv($x,$scale);
2088     # copy private parts over
2089     $x->{_m} = $u->{_m};
2090     $x->{_e} = $u->{_e};
2091     $x->{_es} = $u->{_es};
2092     }
2093   else
2094     {
2095     # calculate the broot() as integer result first, and if it fits, return
2096     # it rightaway (but only if $x and $y are integer):
2097
2098     my $done = 0;                               # not yet
2099     if ($y->is_int() && $x->is_int())
2100       {
2101       my $i = $MBI->_copy( $x->{_m} );
2102       $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2103       my $int = Math::BigInt->bzero();
2104       $int->{value} = $i;
2105       $int->broot($y->as_number());
2106       # if ($exact)
2107       if ($int->copy()->bpow($y) == $x)
2108         {
2109         # found result, return it
2110         $x->{_m} = $int->{value};
2111         $x->{_e} = $MBI->_zero();
2112         $x->{_es} = '+';
2113         $x->bnorm();
2114         $done = 1;
2115         }
2116       }
2117     if ($done == 0)
2118       {
2119       my $u = $self->bone()->bdiv($y,$scale+4);
2120       delete $u->{_a}; delete $u->{_p};         # otherwise it conflicts
2121       $x->bpow($u,$scale+4);                    # el cheapo
2122       }
2123     }
2124   $x->bneg() if $sign == 1;
2125
2126   # shortcut to not run through _find_round_parameters again
2127   if (defined $params[0])
2128     {
2129     $x->bround($params[0],$params[2]);          # then round accordingly
2130     }
2131   else
2132     {
2133     $x->bfround($params[1],$params[2]);         # then round accordingly
2134     }
2135   if ($fallback)
2136     {
2137     # clear a/p after round, since user did not request it
2138     delete $x->{_a}; delete $x->{_p};
2139     }
2140   # restore globals
2141   $$abr = $ab; $$pbr = $pb;
2142   $x;
2143   }
2144
2145 sub bsqrt
2146   {
2147   # calculate square root
2148   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2149
2150   return $x if $x->modify('bsqrt');
2151
2152   return $x->bnan() if $x->{sign} !~ /^[+]/;    # NaN, -inf or < 0
2153   return $x if $x->{sign} eq '+inf';            # sqrt(inf) == inf
2154   return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
2155
2156   # we need to limit the accuracy to protect against overflow
2157   my $fallback = 0;
2158   my (@params,$scale);
2159   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2160
2161   return $x if $x->is_nan();            # error in _find_round_parameters?
2162
2163   # no rounding at all, so must use fallback
2164   if (scalar @params == 0)
2165     {
2166     # simulate old behaviour
2167     $params[0] = $self->div_scale();    # and round to it as accuracy
2168     $scale = $params[0]+4;              # at least four more for proper round
2169     $params[2] = $r;                    # round mode by caller or undef
2170     $fallback = 1;                      # to clear a/p afterwards
2171     }
2172   else
2173     {
2174     # the 4 below is empirical, and there might be cases where it is not
2175     # enough...
2176     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2177     }
2178
2179   # when user set globals, they would interfere with our calculation, so
2180   # disable them and later re-enable them
2181   no strict 'refs';
2182   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2183   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2184   # we also need to disable any set A or P on $x (_find_round_parameters took
2185   # them already into account), since these would interfere, too
2186   delete $x->{_a}; delete $x->{_p};
2187   # need to disable $upgrade in BigInt, to avoid deep recursion
2188   local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2189
2190   my $i = $MBI->_copy( $x->{_m} );
2191   $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2192   my $xas = Math::BigInt->bzero();
2193   $xas->{value} = $i;
2194
2195   my $gs = $xas->copy()->bsqrt();       # some guess
2196
2197   if (($x->{_es} ne '-')                # guess can't be accurate if there are
2198                                         # digits after the dot
2199    && ($xas->bacmp($gs * $gs) == 0))    # guess hit the nail on the head?
2200     {
2201     # exact result, copy result over to keep $x
2202     $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
2203     $x->bnorm();
2204     # shortcut to not run through _find_round_parameters again
2205     if (defined $params[0])
2206       {
2207       $x->bround($params[0],$params[2]);        # then round accordingly
2208       }
2209     else
2210       {
2211       $x->bfround($params[1],$params[2]);       # then round accordingly
2212       }
2213     if ($fallback)
2214       {
2215       # clear a/p after round, since user did not request it
2216       delete $x->{_a}; delete $x->{_p};
2217       }
2218     # re-enable A and P, upgrade is taken care of by "local"
2219     ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
2220     return $x;
2221     }
2222
2223   # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
2224   # of the result by multiplying the input by 100 and then divide the integer
2225   # result of sqrt(input) by 10. Rounding afterwards returns the real result.
2226
2227   # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
2228   my $y1 = $MBI->_copy($x->{_m});
2229
2230   my $length = $MBI->_len($y1);
2231
2232   # Now calculate how many digits the result of sqrt(y1) would have
2233   my $digits = int($length / 2);
2234
2235   # But we need at least $scale digits, so calculate how many are missing
2236   my $shift = $scale - $digits;
2237
2238   # This happens if the input had enough digits
2239   # (we take care of integer guesses above)
2240   $shift = 0 if $shift < 0;
2241
2242   # Multiply in steps of 100, by shifting left two times the "missing" digits
2243   my $s2 = $shift * 2;
2244
2245   # We now make sure that $y1 has the same odd or even number of digits than
2246   # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
2247   # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
2248   # steps of 10. The length of $x does not count, since an even or odd number
2249   # of digits before the dot is not changed by adding an even number of digits
2250   # after the dot (the result is still odd or even digits long).
2251   $s2++ if $MBI->_is_odd($x->{_e});
2252
2253   $MBI->_lsft( $y1, $MBI->_new($s2), 10);
2254
2255   # now take the square root and truncate to integer
2256   $y1 = $MBI->_sqrt($y1);
2257
2258   # By "shifting" $y1 right (by creating a negative _e) we calculate the final
2259   # result, which is than later rounded to the desired scale.
2260
2261   # calculate how many zeros $x had after the '.' (or before it, depending
2262   # on sign of $dat, the result should have half as many:
2263   my $dat = $MBI->_num($x->{_e});
2264   $dat = -$dat if $x->{_es} eq '-';
2265   $dat += $length;
2266
2267   if ($dat > 0)
2268     {
2269     # no zeros after the dot (e.g. 1.23, 0.49 etc)
2270     # preserve half as many digits before the dot than the input had
2271     # (but round this "up")
2272     $dat = int(($dat+1)/2);
2273     }
2274   else
2275     {
2276     $dat = int(($dat)/2);
2277     }
2278   $dat -= $MBI->_len($y1);
2279   if ($dat < 0)
2280     {
2281     $dat = abs($dat);
2282     $x->{_e} = $MBI->_new( $dat );
2283     $x->{_es} = '-';
2284     }
2285   else
2286     {
2287     $x->{_e} = $MBI->_new( $dat );
2288     $x->{_es} = '+';
2289     }
2290   $x->{_m} = $y1;
2291   $x->bnorm();
2292
2293   # shortcut to not run through _find_round_parameters again
2294   if (defined $params[0])
2295     {
2296     $x->bround($params[0],$params[2]);          # then round accordingly
2297     }
2298   else
2299     {
2300     $x->bfround($params[1],$params[2]);         # then round accordingly
2301     }
2302   if ($fallback)
2303     {
2304     # clear a/p after round, since user did not request it
2305     delete $x->{_a}; delete $x->{_p};
2306     }
2307   # restore globals
2308   $$abr = $ab; $$pbr = $pb;
2309   $x;
2310   }
2311
2312 sub bfac
2313   {
2314   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2315   # compute factorial number, modifies first argument
2316
2317   # set up parameters
2318   my ($self,$x,@r) = (ref($_[0]),@_);
2319   # objectify is costly, so avoid it
2320   ($self,$x,@r) = objectify(1,@_) if !ref($x);
2321
2322   # inf => inf
2323   return $x if $x->modify('bfac') || $x->{sign} eq '+inf';
2324
2325   return $x->bnan()
2326     if (($x->{sign} ne '+') ||          # inf, NaN, <0 etc => NaN
2327      ($x->{_es} ne '+'));               # digits after dot?
2328
2329   # use BigInt's bfac() for faster calc
2330   if (! $MBI->_is_zero($x->{_e}))
2331     {
2332     $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
2333     $x->{_e} = $MBI->_zero();           # normalize
2334     $x->{_es} = '+';
2335     }
2336   $MBI->_fac($x->{_m});                 # calculate factorial
2337   $x->bnorm()->round(@r);               # norm again and round result
2338   }
2339
2340 sub _pow
2341   {
2342   # Calculate a power where $y is a non-integer, like 2 ** 0.3
2343   my ($x,$y,@r) = @_;
2344   my $self = ref($x);
2345
2346   # if $y == 0.5, it is sqrt($x)
2347   $HALF = $self->new($HALF) unless ref($HALF);
2348   return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0;
2349
2350   # Using:
2351   # a ** x == e ** (x * ln a)
2352
2353   # u = y * ln x
2354   #                _                         _
2355   # Taylor:       |   u    u^2    u^3         |
2356   # x ** y  = 1 + |  --- + --- + ----- + ...  |
2357   #               |_  1    1*2   1*2*3       _|
2358
2359   # we need to limit the accuracy to protect against overflow
2360   my $fallback = 0;
2361   my ($scale,@params);
2362   ($x,@params) = $x->_find_round_parameters(@r);
2363
2364   return $x if $x->is_nan();            # error in _find_round_parameters?
2365
2366   # no rounding at all, so must use fallback
2367   if (scalar @params == 0)
2368     {
2369     # simulate old behaviour
2370     $params[0] = $self->div_scale();    # and round to it as accuracy
2371     $params[1] = undef;                 # disable P
2372     $scale = $params[0]+4;              # at least four more for proper round
2373     $params[2] = $r[2];                 # round mode by caller or undef
2374     $fallback = 1;                      # to clear a/p afterwards
2375     }
2376   else
2377     {
2378     # the 4 below is empirical, and there might be cases where it is not
2379     # enough...
2380     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2381     }
2382
2383   # when user set globals, they would interfere with our calculation, so
2384   # disable them and later re-enable them
2385   no strict 'refs';
2386   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2387   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2388   # we also need to disable any set A or P on $x (_find_round_parameters took
2389   # them already into account), since these would interfere, too
2390   delete $x->{_a}; delete $x->{_p};
2391   # need to disable $upgrade in BigInt, to avoid deep recursion
2392   local $Math::BigInt::upgrade = undef;
2393
2394   my ($limit,$v,$u,$below,$factor,$next,$over);
2395
2396   $u = $x->copy()->blog(undef,$scale)->bmul($y);
2397   $v = $self->bone();                           # 1
2398   $factor = $self->new(2);                      # 2
2399   $x->bone();                                   # first term: 1
2400
2401   $below = $v->copy();
2402   $over = $u->copy();
2403
2404   $limit = $self->new("1E-". ($scale-1));
2405   #my $steps = 0;
2406   while (3 < 5)
2407     {
2408     # we calculate the next term, and add it to the last
2409     # when the next term is below our limit, it won't affect the outcome
2410     # anymore, so we stop:
2411     $next = $over->copy()->bdiv($below,$scale);
2412     last if $next->bacmp($limit) <= 0;
2413     $x->badd($next);
2414     # calculate things for the next term
2415     $over *= $u; $below *= $factor; $factor->binc();
2416
2417     last if $x->{sign} !~ /^[-+]$/;
2418
2419     #$steps++;
2420     }
2421
2422   # shortcut to not run through _find_round_parameters again
2423   if (defined $params[0])
2424     {
2425     $x->bround($params[0],$params[2]);          # then round accordingly
2426     }
2427   else
2428     {
2429     $x->bfround($params[1],$params[2]);         # then round accordingly
2430     }
2431   if ($fallback)
2432     {
2433     # clear a/p after round, since user did not request it
2434     delete $x->{_a}; delete $x->{_p};
2435     }
2436   # restore globals
2437   $$abr = $ab; $$pbr = $pb;
2438   $x;
2439   }
2440
2441 sub bpow
2442   {
2443   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2444   # compute power of two numbers, second arg is used as integer
2445   # modifies first argument
2446
2447   # set up parameters
2448   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2449   # objectify is costly, so avoid it
2450   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2451     {
2452     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2453     }
2454
2455   return $x if $x->modify('bpow');
2456
2457   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
2458   return $x if $x->{sign} =~ /^[+-]inf$/;
2459
2460   # cache the result of is_zero
2461   my $y_is_zero = $y->is_zero();
2462   return $x->bone() if $y_is_zero;
2463   return $x         if $x->is_one() || $y->is_one();
2464
2465   my $x_is_zero = $x->is_zero();
2466   return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int();         # non-integer power
2467
2468   my $y1 = $y->as_number()->{value};                    # make MBI part
2469
2470   # if ($x == -1)
2471   if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
2472     {
2473     # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
2474     return $MBI->_is_odd($y1) ? $x : $x->babs(1);
2475     }
2476   if ($x_is_zero)
2477     {
2478     return $x if $y->{sign} eq '+';     # 0**y => 0 (if not y <= 0)
2479     # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
2480     return $x->binf();
2481     }
2482
2483   my $new_sign = '+';
2484   $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
2485
2486   # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
2487   $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
2488   $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
2489
2490   $x->{sign} = $new_sign;
2491   $x->bnorm();
2492   if ($y->{sign} eq '-')
2493     {
2494     # modify $x in place!
2495     my $z = $x->copy(); $x->bone();
2496     return scalar $x->bdiv($z,$a,$p,$r);        # round in one go (might ignore y's A!)
2497     }
2498   $x->round($a,$p,$r,$y);
2499   }
2500
2501 sub bmodpow
2502   {
2503   # takes a very large number to a very large exponent in a given very
2504   # large modulus, quickly, thanks to binary exponentiation. Supports
2505   # negative exponents.
2506   my ($self,$num,$exp,$mod,@r) = objectify(3,@_);
2507
2508   return $num if $num->modify('bmodpow');
2509
2510   # check modulus for valid values
2511   return $num->bnan() if ($mod->{sign} ne '+'           # NaN, - , -inf, +inf
2512                        || $mod->is_zero());
2513
2514   # check exponent for valid values
2515   if ($exp->{sign} =~ /\w/)
2516     {
2517     # i.e., if it's NaN, +inf, or -inf...
2518     return $num->bnan();
2519     }
2520
2521   $num->bmodinv ($mod) if ($exp->{sign} eq '-');
2522
2523   # check num for valid values (also NaN if there was no inverse but $exp < 0)
2524   return $num->bnan() if $num->{sign} !~ /^[+-]$/;
2525
2526   # $mod is positive, sign on $exp is ignored, result also positive
2527
2528   # XXX TODO: speed it up when all three numbers are integers
2529   $num->bpow($exp)->bmod($mod);
2530   }
2531
2532 ###############################################################################
2533 # trigonometric functions
2534
2535 # helper function for bpi() and batan2(), calculates arcus tanges (1/x)
2536
2537 sub _atan_inv
2538   {
2539   # return a/b so that a/b approximates atan(1/x) to at least limit digits
2540   my ($self, $x, $limit) = @_;
2541
2542   # Taylor:       x^3   x^5   x^7   x^9
2543   #    atan = x - --- + --- - --- + --- - ...
2544   #                3     5     7     9
2545
2546   #               1      1         1        1
2547   #    atan 1/x = - - ------- + ------- - ------- + ...
2548   #               x   x^3 * 3   x^5 * 5   x^7 * 7
2549
2550   #               1      1         1            1
2551   #    atan 1/x = - - --------- + ---------- - ----------- + ...
2552   #               5    3 * 125     5 * 3125     7 * 78125
2553
2554   # Subtraction/addition of a rational:
2555
2556   #  5    7    5*3 +- 7*4
2557   #  - +- -  = ----------
2558   #  4    3       4*3
2559
2560   # Term:  N        N+1
2561   #
2562   #        a             1                  a * d * c +- b
2563   #        ----- +- ------------------  =  ----------------
2564   #        b           d * c                b * d * c
2565
2566   #  since b1 = b0 * (d-2) * c
2567
2568   #        a             1                  a * d +- b / c
2569   #        ----- +- ------------------  =  ----------------
2570   #        b           d * c                b * d
2571
2572   # and  d = d + 2
2573   # and  c = c * x * x
2574
2575   #        u = d * c
2576   #        stop if length($u) > limit
2577   #        a = a * u +- b
2578   #        b = b * u
2579   #        d = d + 2
2580   #        c = c * x * x
2581   #        sign = 1 - sign
2582
2583   my $a = $MBI->_one();
2584   my $b = $MBI->_copy($x);
2585
2586   my $x2  = $MBI->_mul( $MBI->_copy($x), $b);           # x2 = x * x
2587   my $d   = $MBI->_new( 3 );                            # d = 3
2588   my $c   = $MBI->_mul( $MBI->_copy($x), $x2);          # c = x ^ 3
2589   my $two = $MBI->_new( 2 );
2590
2591   # run the first step unconditionally
2592   my $u = $MBI->_mul( $MBI->_copy($d), $c);
2593   $a = $MBI->_mul($a, $u);
2594   $a = $MBI->_sub($a, $b);
2595   $b = $MBI->_mul($b, $u);
2596   $d = $MBI->_add($d, $two);
2597   $c = $MBI->_mul($c, $x2);
2598
2599   # a is now a * (d-3) * c
2600   # b is now b * (d-2) * c
2601
2602   # run the second step unconditionally
2603   $u = $MBI->_mul( $MBI->_copy($d), $c);
2604   $a = $MBI->_mul($a, $u);
2605   $a = $MBI->_add($a, $b);
2606   $b = $MBI->_mul($b, $u);
2607   $d = $MBI->_add($d, $two);
2608   $c = $MBI->_mul($c, $x2);
2609
2610   # a is now a * (d-3) * (d-5) * c * c
2611   # b is now b * (d-2) * (d-4) * c * c
2612
2613   # so we can remove c * c from both a and b to shorten the numbers involved:
2614   $a = $MBI->_div($a, $x2);
2615   $b = $MBI->_div($b, $x2);
2616   $a = $MBI->_div($a, $x2);
2617   $b = $MBI->_div($b, $x2);
2618
2619 #  my $step = 0;
2620   my $sign = 0;                                         # 0 => -, 1 => +
2621   while (3 < 5)
2622     {
2623 #    $step++;
2624 #    if (($i++ % 100) == 0)
2625 #      {
2626 #    print "a=",$MBI->_str($a),"\n";
2627 #    print "b=",$MBI->_str($b),"\n";
2628 #      }
2629 #    print "d=",$MBI->_str($d),"\n";
2630 #    print "x2=",$MBI->_str($x2),"\n";
2631 #    print "c=",$MBI->_str($c),"\n";
2632
2633     my $u = $MBI->_mul( $MBI->_copy($d), $c);
2634     # use _alen() for libs like GMP where _len() would be O(N^2)
2635     last if $MBI->_alen($u) > $limit;
2636     my ($bc,$r) = $MBI->_div( $MBI->_copy($b), $c);
2637     if ($MBI->_is_zero($r))
2638       {
2639       # b / c is an integer, so we can remove c from all terms
2640       # this happens almost every time:
2641       $a = $MBI->_mul($a, $d);
2642       $a = $MBI->_sub($a, $bc) if $sign == 0;
2643       $a = $MBI->_add($a, $bc) if $sign == 1;
2644       $b = $MBI->_mul($b, $d);
2645       }
2646     else
2647       {
2648       # b / c is not an integer, so we keep c in the terms
2649       # this happens very rarely, for instance for x = 5, this happens only
2650       # at the following steps:
2651       # 1, 5, 14, 32, 72, 157, 340, ...
2652       $a = $MBI->_mul($a, $u);
2653       $a = $MBI->_sub($a, $b) if $sign == 0;
2654       $a = $MBI->_add($a, $b) if $sign == 1;
2655       $b = $MBI->_mul($b, $u);
2656       }
2657     $d = $MBI->_add($d, $two);
2658     $c = $MBI->_mul($c, $x2);
2659     $sign = 1 - $sign;
2660
2661     }
2662
2663 #  print "Took $step steps for ", $MBI->_str($x),"\n";
2664 #  print "a=",$MBI->_str($a),"\n"; print "b=",$MBI->_str($b),"\n";
2665   # return a/b so that a/b approximates atan(1/x)
2666   ($a,$b);
2667   }
2668
2669 sub bpi
2670   {
2671   my ($self,$n) = @_;
2672   if (@_ == 0)
2673     {
2674     $self = $class;
2675     }
2676   if (@_ == 1)
2677     {
2678     # called like Math::BigFloat::bpi(10);
2679     $n = $self; $self = $class;
2680     # called like Math::BigFloat->bpi();
2681     $n = undef if $n eq 'Math::BigFloat';
2682     }
2683   $self = ref($self) if ref($self);
2684   my $fallback = defined $n ? 0 : 1;
2685   $n = 40 if !defined $n || $n < 1;
2686
2687   # after 黃見利 (Hwang Chien-Lih) (1997)
2688   # pi/4 = 183 * atan(1/239) + 32 * atan(1/1023) – 68 * atan(1/5832)
2689   #      + 12 * atan(1/110443) - 12 * atan(1/4841182) - 100 * atan(1/6826318)
2690
2691   # a few more to prevent rounding errors
2692   $n += 4;
2693
2694   my ($a,$b) = $self->_atan_inv( $MBI->_new(239),$n);
2695   my ($c,$d) = $self->_atan_inv( $MBI->_new(1023),$n);
2696   my ($e,$f) = $self->_atan_inv( $MBI->_new(5832),$n);
2697   my ($g,$h) = $self->_atan_inv( $MBI->_new(110443),$n);
2698   my ($i,$j) = $self->_atan_inv( $MBI->_new(4841182),$n);
2699   my ($k,$l) = $self->_atan_inv( $MBI->_new(6826318),$n);
2700
2701   $MBI->_mul($a, $MBI->_new(732));
2702   $MBI->_mul($c, $MBI->_new(128));
2703   $MBI->_mul($e, $MBI->_new(272));
2704   $MBI->_mul($g, $MBI->_new(48));
2705   $MBI->_mul($i, $MBI->_new(48));
2706   $MBI->_mul($k, $MBI->_new(400));
2707
2708   my $x = $self->bone(); $x->{_m} = $a; my $x_d = $self->bone(); $x_d->{_m} = $b;
2709   my $y = $self->bone(); $y->{_m} = $c; my $y_d = $self->bone(); $y_d->{_m} = $d;
2710   my $z = $self->bone(); $z->{_m} = $e; my $z_d = $self->bone(); $z_d->{_m} = $f;
2711   my $u = $self->bone(); $u->{_m} = $g; my $u_d = $self->bone(); $u_d->{_m} = $h;
2712   my $v = $self->bone(); $v->{_m} = $i; my $v_d = $self->bone(); $v_d->{_m} = $j;
2713   my $w = $self->bone(); $w->{_m} = $k; my $w_d = $self->bone(); $w_d->{_m} = $l;
2714   $x->bdiv($x_d, $n);
2715   $y->bdiv($y_d, $n);
2716   $z->bdiv($z_d, $n);
2717   $u->bdiv($u_d, $n);
2718   $v->bdiv($v_d, $n);
2719   $w->bdiv($w_d, $n);
2720
2721   delete $x->{_a}; delete $y->{_a}; delete $z->{_a};
2722   delete $u->{_a}; delete $v->{_a}; delete $w->{_a};
2723   $x->badd($y)->bsub($z)->badd($u)->bsub($v)->bsub($w);
2724
2725   $x->bround($n-4);
2726   delete $x->{_a} if $fallback == 1;
2727   $x;
2728   }
2729
2730 sub bcos
2731   {
2732   # Calculate a cosinus of x.
2733   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2734
2735   # Taylor:      x^2   x^4   x^6   x^8
2736   #    cos = 1 - --- + --- - --- + --- ...
2737   #               2!    4!    6!    8!
2738
2739   # we need to limit the accuracy to protect against overflow
2740   my $fallback = 0;
2741   my ($scale,@params);
2742   ($x,@params) = $x->_find_round_parameters(@r);
2743
2744   #         constant object       or error in _find_round_parameters?
2745   return $x if $x->modify('bcos') || $x->is_nan();
2746
2747   return $x->bone(@r) if $x->is_zero();
2748
2749   # no rounding at all, so must use fallback
2750   if (scalar @params == 0)
2751     {
2752     # simulate old behaviour
2753     $params[0] = $self->div_scale();    # and round to it as accuracy
2754     $params[1] = undef;                 # disable P
2755     $scale = $params[0]+4;              # at least four more for proper round
2756     $params[2] = $r[2];                 # round mode by caller or undef
2757     $fallback = 1;                      # to clear a/p afterwards
2758     }
2759   else
2760     {
2761     # the 4 below is empirical, and there might be cases where it is not
2762     # enough...
2763     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2764     }
2765
2766   # when user set globals, they would interfere with our calculation, so
2767   # disable them and later re-enable them
2768   no strict 'refs';
2769   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2770   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2771   # we also need to disable any set A or P on $x (_find_round_parameters took
2772   # them already into account), since these would interfere, too
2773   delete $x->{_a}; delete $x->{_p};
2774   # need to disable $upgrade in BigInt, to avoid deep recursion
2775   local $Math::BigInt::upgrade = undef;
2776
2777   my $last = 0;
2778   my $over = $x * $x;                   # X ^ 2
2779   my $x2 = $over->copy();               # X ^ 2; difference between terms
2780   my $sign = 1;                         # start with -=
2781   my $below = $self->new(2); my $factorial = $self->new(3);
2782   $x->bone(); delete $x->{_a}; delete $x->{_p};
2783
2784   my $limit = $self->new("1E-". ($scale-1));
2785   #my $steps = 0;
2786   while (3 < 5)
2787     {
2788     # we calculate the next term, and add it to the last
2789     # when the next term is below our limit, it won't affect the outcome
2790     # anymore, so we stop:
2791     my $next = $over->copy()->bdiv($below,$scale);
2792     last if $next->bacmp($limit) <= 0;
2793
2794     if ($sign == 0)
2795       {
2796       $x->badd($next);
2797       }
2798     else
2799       {
2800       $x->bsub($next);
2801       }
2802     $sign = 1-$sign;                                    # alternate
2803     # calculate things for the next term
2804     $over->bmul($x2);                                   # $x*$x
2805     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
2806     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
2807     }
2808
2809   # shortcut to not run through _find_round_parameters again
2810   if (defined $params[0])
2811     {
2812     $x->bround($params[0],$params[2]);          # then round accordingly
2813     }
2814   else
2815     {
2816     $x->bfround($params[1],$params[2]);         # then round accordingly
2817     }
2818   if ($fallback)
2819     {
2820     # clear a/p after round, since user did not request it
2821     delete $x->{_a}; delete $x->{_p};
2822     }
2823   # restore globals
2824   $$abr = $ab; $$pbr = $pb;
2825   $x;
2826   }
2827
2828 sub bsin
2829   {
2830   # Calculate a sinus of x.
2831   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2832
2833   # taylor:      x^3   x^5   x^7   x^9
2834   #    sin = x - --- + --- - --- + --- ...
2835   #               3!    5!    7!    9!
2836
2837   # we need to limit the accuracy to protect against overflow
2838   my $fallback = 0;
2839   my ($scale,@params);
2840   ($x,@params) = $x->_find_round_parameters(@r);
2841
2842   #         constant object       or error in _find_round_parameters?
2843   return $x if $x->modify('bsin') || $x->is_nan();
2844
2845   return $x->bzero(@r) if $x->is_zero();
2846
2847   # no rounding at all, so must use fallback
2848   if (scalar @params == 0)
2849     {
2850     # simulate old behaviour
2851     $params[0] = $self->div_scale();    # and round to it as accuracy
2852     $params[1] = undef;                 # disable P
2853     $scale = $params[0]+4;              # at least four more for proper round
2854     $params[2] = $r[2];                 # round mode by caller or undef
2855     $fallback = 1;                      # to clear a/p afterwards
2856     }
2857   else
2858     {
2859     # the 4 below is empirical, and there might be cases where it is not
2860     # enough...
2861     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2862     }
2863
2864   # when user set globals, they would interfere with our calculation, so
2865   # disable them and later re-enable them
2866   no strict 'refs';
2867   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2868   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2869   # we also need to disable any set A or P on $x (_find_round_parameters took
2870   # them already into account), since these would interfere, too
2871   delete $x->{_a}; delete $x->{_p};
2872   # need to disable $upgrade in BigInt, to avoid deep recursion
2873   local $Math::BigInt::upgrade = undef;
2874
2875   my $last = 0;
2876   my $over = $x * $x;                   # X ^ 2
2877   my $x2 = $over->copy();               # X ^ 2; difference between terms
2878   $over->bmul($x);                      # X ^ 3 as starting value
2879   my $sign = 1;                         # start with -=
2880   my $below = $self->new(6); my $factorial = $self->new(4);
2881   delete $x->{_a}; delete $x->{_p};
2882
2883   my $limit = $self->new("1E-". ($scale-1));
2884   #my $steps = 0;
2885   while (3 < 5)
2886     {
2887     # we calculate the next term, and add it to the last
2888     # when the next term is below our limit, it won't affect the outcome
2889     # anymore, so we stop:
2890     my $next = $over->copy()->bdiv($below,$scale);
2891     last if $next->bacmp($limit) <= 0;
2892
2893     if ($sign == 0)
2894       {
2895       $x->badd($next);
2896       }
2897     else
2898       {
2899       $x->bsub($next);
2900       }
2901     $sign = 1-$sign;                                    # alternate
2902     # calculate things for the next term
2903     $over->bmul($x2);                                   # $x*$x
2904     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
2905     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
2906     }
2907
2908   # shortcut to not run through _find_round_parameters again
2909   if (defined $params[0])
2910     {
2911     $x->bround($params[0],$params[2]);          # then round accordingly
2912     }
2913   else
2914     {
2915     $x->bfround($params[1],$params[2]);         # then round accordingly
2916     }
2917   if ($fallback)
2918     {
2919     # clear a/p after round, since user did not request it
2920     delete $x->{_a}; delete $x->{_p};
2921     }
2922   # restore globals
2923   $$abr = $ab; $$pbr = $pb;
2924   $x;
2925   }
2926
2927 sub batan2
2928   {
2929   # calculate arcus tangens of ($y/$x)
2930
2931   # set up parameters
2932   my ($self,$y,$x,@r) = (ref($_[0]),@_);
2933   # objectify is costly, so avoid it
2934   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2935     {
2936     ($self,$y,$x,@r) = objectify(2,@_);
2937     }
2938
2939   return $y if $y->modify('batan2');
2940
2941   return $y->bnan() if ($y->{sign} eq $nan) || ($x->{sign} eq $nan);
2942
2943   # Y X
2944   # 0 0 result is 0
2945   # 0 +x result is 0
2946   # ? inf result is 0
2947   return $y->bzero(@r) if ($x->is_inf('+') && !$y->is_inf()) || ($y->is_zero() && $x->{sign} eq '+');
2948
2949   # Y    X
2950   # != 0 -inf result is +- pi
2951   if ($x->is_inf() || $y->is_inf())
2952     {
2953     # calculate PI
2954     my $pi = $self->bpi(@r);
2955     if ($y->is_inf())
2956       {
2957       # upgrade to BigRat etc.
2958       return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
2959       if ($x->{sign} eq '-inf')
2960         {
2961         # calculate 3 pi/4
2962         $MBI->_mul($pi->{_m}, $MBI->_new(3));
2963         $MBI->_div($pi->{_m}, $MBI->_new(4));
2964         }
2965       elsif ($x->{sign} eq '+inf')
2966         {
2967         # calculate pi/4
2968         $MBI->_div($pi->{_m}, $MBI->_new(4));
2969         }
2970       else
2971         {
2972         # calculate pi/2
2973         $MBI->_div($pi->{_m}, $MBI->_new(2));
2974         }
2975       $y->{sign} = substr($y->{sign},0,1); # keep +/-
2976       }
2977     # modify $y in place
2978     $y->{_m} = $pi->{_m};
2979     $y->{_e} = $pi->{_e};
2980     $y->{_es} = $pi->{_es};
2981     # keep the sign of $y
2982     return $y;
2983     }
2984
2985   return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
2986
2987   # Y X
2988   # 0 -x result is PI
2989   if ($y->is_zero())
2990     {
2991     # calculate PI
2992     my $pi = $self->bpi(@r);
2993     # modify $y in place
2994     $y->{_m} = $pi->{_m};
2995     $y->{_e} = $pi->{_e};
2996     $y->{_es} = $pi->{_es};
2997     $y->{sign} = '+';
2998     return $y;
2999     }
3000
3001   # Y X
3002   # +y 0 result is PI/2
3003   # -y 0 result is -PI/2
3004   if ($x->is_zero())
3005     {
3006     # calculate PI/2
3007     my $pi = $self->bpi(@r);
3008     # modify $y in place
3009     $y->{_m} = $pi->{_m};
3010     $y->{_e} = $pi->{_e};
3011     $y->{_es} = $pi->{_es};
3012     # -y => -PI/2, +y => PI/2
3013     $MBI->_div($y->{_m}, $MBI->_new(2));
3014     return $y;
3015     }
3016
3017   # we need to limit the accuracy to protect against overflow
3018   my $fallback = 0;
3019   my ($scale,@params);
3020   ($y,@params) = $y->_find_round_parameters(@r);
3021
3022   # error in _find_round_parameters?
3023   return $y if $y->is_nan();
3024
3025   # no rounding at all, so must use fallback
3026   if (scalar @params == 0)
3027     {
3028     # simulate old behaviour
3029     $params[0] = $self->div_scale();    # and round to it as accuracy
3030     $params[1] = undef;                 # disable P
3031     $scale = $params[0]+4;              # at least four more for proper round
3032     $params[2] = $r[2];                 # round mode by caller or undef
3033     $fallback = 1;                      # to clear a/p afterwards
3034     }
3035   else
3036     {
3037     # the 4 below is empirical, and there might be cases where it is not
3038     # enough...
3039     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3040     }
3041
3042   # inlined is_one() && is_one('-')
3043   if ($MBI->_is_one($y->{_m}) && $MBI->_is_zero($y->{_e}))
3044     {
3045     # shortcut: 1 1 result is PI/4
3046     # inlined is_one() && is_one('-')
3047     if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3048       {
3049       # 1,1 => PI/4
3050       my $pi_4 = $self->bpi( $scale - 3);
3051       # modify $y in place
3052       $y->{_m} = $pi_4->{_m};
3053       $y->{_e} = $pi_4->{_e};
3054       $y->{_es} = $pi_4->{_es};
3055       # 1 1 => +
3056       # -1 1 => -
3057       # 1 -1 => -
3058       # -1 -1 => +
3059       $y->{sign} = $x->{sign} eq $y->{sign} ? '+' : '-';
3060       $MBI->_div($y->{_m}, $MBI->_new(4));
3061       return $y;
3062       }
3063     # shortcut: 1 int(X) result is _atan_inv(X)
3064
3065     # is integer
3066     if ($x->{_es} eq '+')
3067       {
3068       my $x1 = $MBI->_copy($x->{_m});
3069       $MBI->_lsft($x1, $x->{_e},10) unless $MBI->_is_zero($x->{_e});
3070
3071       my ($a,$b) = $self->_atan_inv($x1, $scale);
3072       my $y_sign = $y->{sign};
3073       # calculate A/B
3074       $y->bone(); $y->{_m} = $a; my $y_d = $self->bone(); $y_d->{_m} = $b;
3075       $y->bdiv($y_d, @r);
3076       $y->{sign} = $y_sign;
3077       return $y;
3078       }
3079     }
3080
3081   # handle all other cases
3082   #  X  Y
3083   # +x +y 0 to PI/2
3084   # -x +y PI/2 to PI
3085   # +x -y 0 to -PI/2
3086   # -x -y -PI/2 to -PI
3087
3088   my $y_sign = $y->{sign};
3089
3090   # divide $x by $y
3091   $y->bdiv($x, $scale) unless $x->is_one();
3092   $y->batan(@r);
3093
3094   # restore sign
3095   $y->{sign} = $y_sign;
3096
3097   $y;
3098   }
3099
3100 sub batan
3101   {
3102   # Calculate a arcus tangens of x.
3103   my ($x,@r) = @_;
3104   my $self = ref($x);
3105
3106   # taylor:       x^3   x^5   x^7   x^9
3107   #    atan = x - --- + --- - --- + --- ...
3108   #                3     5     7     9
3109
3110   # we need to limit the accuracy to protect against overflow
3111   my $fallback = 0;
3112   my ($scale,@params);
3113   ($x,@params) = $x->_find_round_parameters(@r);
3114
3115   #         constant object       or error in _find_round_parameters?
3116   return $x if $x->modify('batan') || $x->is_nan();
3117
3118   if ($x->{sign} =~ /^[+-]inf\z/)
3119     {
3120     # +inf result is PI/2
3121     # -inf result is -PI/2
3122     # calculate PI/2
3123     my $pi = $self->bpi(@r);
3124     # modify $x in place
3125     $x->{_m} = $pi->{_m};
3126     $x->{_e} = $pi->{_e};
3127     $x->{_es} = $pi->{_es};
3128     # -y => -PI/2, +y => PI/2
3129     $x->{sign} = substr($x->{sign},0,1);                # +inf => +
3130     $MBI->_div($x->{_m}, $MBI->_new(2));
3131     return $x;
3132     }
3133
3134   return $x->bzero(@r) if $x->is_zero();
3135
3136   # no rounding at all, so must use fallback
3137   if (scalar @params == 0)
3138     {
3139     # simulate old behaviour
3140     $params[0] = $self->div_scale();    # and round to it as accuracy
3141     $params[1] = undef;                 # disable P
3142     $scale = $params[0]+4;              # at least four more for proper round
3143     $params[2] = $r[2];                 # round mode by caller or undef
3144     $fallback = 1;                      # to clear a/p afterwards
3145     }
3146   else
3147     {
3148     # the 4 below is empirical, and there might be cases where it is not
3149     # enough...
3150     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3151     }
3152
3153   # 1 or -1 => PI/4
3154   # inlined is_one() && is_one('-')
3155   if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3156     {
3157     my $pi = $self->bpi($scale - 3);
3158     # modify $x in place
3159     $x->{_m} = $pi->{_m};
3160     $x->{_e} = $pi->{_e};
3161     $x->{_es} = $pi->{_es};
3162     # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4)
3163     $MBI->_div($x->{_m}, $MBI->_new(4));
3164     return $x;
3165     }
3166
3167   # This series is only valid if -1 < x < 1, so for other x we need to
3168   # to calculate PI/2 - atan(1/x):
3169   my $one = $MBI->_new(1);
3170   my $pi = undef;
3171   if ($x->{_es} eq '+' && ($MBI->_acmp($x->{_m},$one) >= 0))
3172     {
3173     # calculate PI/2
3174     $pi = $self->bpi($scale - 3);
3175     $MBI->_div($pi->{_m}, $MBI->_new(2));
3176     # calculate 1/$x:
3177     my $x_copy = $x->copy();
3178     # modify $x in place
3179     $x->bone(); $x->bdiv($x_copy,$scale);
3180     }
3181
3182   # when user set globals, they would interfere with our calculation, so
3183   # disable them and later re-enable them
3184   no strict 'refs';
3185   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
3186   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
3187   # we also need to disable any set A or P on $x (_find_round_parameters took
3188   # them already into account), since these would interfere, too
3189   delete $x->{_a}; delete $x->{_p};
3190   # need to disable $upgrade in BigInt, to avoid deep recursion
3191   local $Math::BigInt::upgrade = undef;
3192
3193   my $last = 0;
3194   my $over = $x * $x;                   # X ^ 2
3195   my $x2 = $over->copy();               # X ^ 2; difference between terms
3196   $over->bmul($x);                      # X ^ 3 as starting value
3197   my $sign = 1;                         # start with -=
3198   my $below = $self->new(3);
3199   my $two = $self->new(2);
3200   delete $x->{_a}; delete $x->{_p};
3201
3202   my $limit = $self->new("1E-". ($scale-1));
3203   #my $steps = 0;
3204   while (3 < 5)
3205     {
3206     # we calculate the next term, and add it to the last
3207     # when the next term is below our limit, it won't affect the outcome
3208     # anymore, so we stop:
3209     my $next = $over->copy()->bdiv($below,$scale);
3210     last if $next->bacmp($limit) <= 0;
3211
3212     if ($sign == 0)
3213       {
3214       $x->badd($next);
3215       }
3216     else
3217       {
3218       $x->bsub($next);
3219       }
3220     $sign = 1-$sign;                                    # alternate
3221     # calculate things for the next term
3222     $over->bmul($x2);                                   # $x*$x
3223     $below->badd($two);                                 # n += 2
3224     }
3225
3226   if (defined $pi)
3227     {
3228     my $x_copy = $x->copy();
3229     # modify $x in place
3230     $x->{_m} = $pi->{_m};
3231     $x->{_e} = $pi->{_e};
3232     $x->{_es} = $pi->{_es};
3233     # PI/2 - $x
3234     $x->bsub($x_copy);
3235     }
3236
3237   # shortcut to not run through _find_round_parameters again
3238   if (defined $params[0])
3239     {
3240     $x->bround($params[0],$params[2]);          # then round accordingly
3241     }
3242   else
3243     {
3244     $x->bfround($params[1],$params[2]);         # then round accordingly
3245     }
3246   if ($fallback)
3247     {
3248     # clear a/p after round, since user did not request it
3249     delete $x->{_a}; delete $x->{_p};
3250     }
3251   # restore globals
3252   $$abr = $ab; $$pbr = $pb;
3253   $x;
3254   }
3255
3256 ###############################################################################
3257 # rounding functions
3258
3259 sub bfround
3260   {
3261   # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
3262   # $n == 0 means round to integer
3263   # expects and returns normalized numbers!
3264   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
3265
3266   my ($scale,$mode) = $x->_scale_p(@_);
3267   return $x if !defined $scale || $x->modify('bfround'); # no-op
3268
3269   # never round a 0, +-inf, NaN
3270   if ($x->is_zero())
3271     {
3272     $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
3273     return $x;
3274     }
3275   return $x if $x->{sign} !~ /^[+-]$/;
3276
3277   # don't round if x already has lower precision
3278   return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
3279
3280   $x->{_p} = $scale;                    # remember round in any case
3281   delete $x->{_a};                      # and clear A
3282   if ($scale < 0)
3283     {
3284     # round right from the '.'
3285
3286     return $x if $x->{_es} eq '+';              # e >= 0 => nothing to round
3287
3288     $scale = -$scale;                           # positive for simplicity
3289     my $len = $MBI->_len($x->{_m});             # length of mantissa
3290
3291     # the following poses a restriction on _e, but if _e is bigger than a
3292     # scalar, you got other problems (memory etc) anyway
3293     my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e})));   # digits after dot
3294     my $zad = 0;                                # zeros after dot
3295     $zad = $dad - $len if (-$dad < -$len);      # for 0.00..00xxx style
3296
3297     # p rint "scale $scale dad $dad zad $zad len $len\n";
3298     # number  bsstr   len zad dad
3299     # 0.123   123e-3    3   0 3
3300     # 0.0123  123e-4    3   1 4
3301     # 0.001   1e-3      1   2 3
3302     # 1.23    123e-2    3   0 2
3303     # 1.2345  12345e-4  5   0 4
3304
3305     # do not round after/right of the $dad
3306     return $x if $scale > $dad;                 # 0.123, scale >= 3 => exit
3307
3308     # round to zero if rounding inside the $zad, but not for last zero like:
3309     # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
3310     return $x->bzero() if $scale < $zad;
3311     if ($scale == $zad)                 # for 0.006, scale -3 and trunc
3312       {
3313       $scale = -$len;
3314       }
3315     else
3316       {
3317       # adjust round-point to be inside mantissa
3318       if ($zad != 0)
3319         {
3320         $scale = $scale-$zad;
3321         }
3322       else
3323         {
3324         my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;    # digits before dot
3325         $scale = $dbd+$scale;
3326         }
3327       }
3328     }
3329   else
3330     {
3331     # round left from the '.'
3332
3333     # 123 => 100 means length(123) = 3 - $scale (2) => 1
3334
3335     my $dbt = $MBI->_len($x->{_m});
3336     # digits before dot
3337     my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
3338     # should be the same, so treat it as this
3339     $scale = 1 if $scale == 0;
3340     # shortcut if already integer
3341     return $x if $scale == 1 && $dbt <= $dbd;
3342     # maximum digits before dot
3343     ++$dbd;
3344
3345     if ($scale > $dbd)
3346        {
3347        # not enough digits before dot, so round to zero
3348        return $x->bzero;
3349        }
3350     elsif ( $scale == $dbd )
3351        {
3352        # maximum
3353        $scale = -$dbt;
3354        }
3355     else
3356        {
3357        $scale = $dbd - $scale;
3358        }
3359     }
3360   # pass sign to bround for rounding modes '+inf' and '-inf'
3361   my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
3362   $m->bround($scale,$mode);
3363   $x->{_m} = $m->{value};                       # get our mantissa back
3364   $x->bnorm();
3365   }
3366
3367 sub bround
3368   {
3369   # accuracy: preserve $N digits, and overwrite the rest with 0's
3370   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
3371
3372   if (($_[0] || 0) < 0)
3373     {
3374     require Carp; Carp::croak ('bround() needs positive accuracy');
3375     }
3376
3377   my ($scale,$mode) = $x->_scale_a(@_);
3378   return $x if !defined $scale || $x->modify('bround'); # no-op
3379
3380   # scale is now either $x->{_a}, $accuracy, or the user parameter
3381   # test whether $x already has lower accuracy, do nothing in this case
3382   # but do round if the accuracy is the same, since a math operation might
3383   # want to round a number with A=5 to 5 digits afterwards again
3384   return $x if defined $x->{_a} && $x->{_a} < $scale;
3385
3386   # scale < 0 makes no sense
3387   # scale == 0 => keep all digits
3388   # never round a +-inf, NaN
3389   return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
3390
3391   # 1: never round a 0
3392   # 2: if we should keep more digits than the mantissa has, do nothing
3393   if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
3394     {
3395     $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
3396     return $x;
3397     }
3398
3399   # pass sign to bround for '+inf' and '-inf' rounding modes
3400   my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
3401
3402   $m->bround($scale,$mode);             # round mantissa
3403   $x->{_m} = $m->{value};               # get our mantissa back
3404   $x->{_a} = $scale;                    # remember rounding
3405   delete $x->{_p};                      # and clear P
3406   $x->bnorm();                          # del trailing zeros gen. by bround()
3407   }
3408
3409 sub bfloor
3410   {
3411   # return integer less or equal then $x
3412   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3413
3414   return $x if $x->modify('bfloor');
3415
3416   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3417
3418   # if $x has digits after dot
3419   if ($x->{_es} eq '-')
3420     {
3421     $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3422     $x->{_e} = $MBI->_zero();                   # trunc/norm
3423     $x->{_es} = '+';                            # abs e
3424     $MBI->_inc($x->{_m}) if $x->{sign} eq '-';  # increment if negative
3425     }
3426   $x->round($a,$p,$r);
3427   }
3428
3429 sub bceil
3430   {
3431   # return integer greater or equal then $x
3432   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3433
3434   return $x if $x->modify('bceil');
3435   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3436
3437   # if $x has digits after dot
3438   if ($x->{_es} eq '-')
3439     {
3440     $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3441     $x->{_e} = $MBI->_zero();                   # trunc/norm
3442     $x->{_es} = '+';                            # abs e
3443     $MBI->_inc($x->{_m}) if $x->{sign} eq '+';  # increment if positive
3444     }
3445   $x->round($a,$p,$r);
3446   }
3447
3448 sub brsft
3449   {
3450   # shift right by $y (divide by power of $n)
3451
3452   # set up parameters
3453   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3454   # objectify is costly, so avoid it
3455   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3456     {
3457     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3458     }
3459
3460   return $x if $x->modify('brsft');
3461   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3462
3463   $n = 2 if !defined $n; $n = $self->new($n);
3464
3465   # negative amount?
3466   return $x->blsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3467
3468   # the following call to bdiv() will return either quo or (quo,remainder):
3469   $x->bdiv($n->bpow($y),$a,$p,$r,$y);
3470   }
3471
3472 sub blsft
3473   {
3474   # shift left by $y (multiply by power of $n)
3475
3476   # set up parameters
3477   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3478   # objectify is costly, so avoid it
3479   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3480     {
3481     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3482     }
3483
3484   return $x if $x->modify('blsft');
3485   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3486
3487   $n = 2 if !defined $n; $n = $self->new($n);
3488
3489   # negative amount?
3490   return $x->brsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3491
3492   $x->bmul($n->bpow($y),$a,$p,$r,$y);
3493   }
3494
3495 ###############################################################################
3496
3497 sub DESTROY
3498   {
3499   # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
3500   }
3501
3502 sub AUTOLOAD
3503   {
3504   # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
3505   # or falling back to MBI::bxxx()
3506   my $name = $AUTOLOAD;
3507
3508   $name =~ s/(.*):://;  # split package
3509   my $c = $1 || $class;
3510   no strict 'refs';
3511   $c->import() if $IMPORT == 0;
3512   if (!_method_alias($name))
3513     {
3514     if (!defined $name)
3515       {
3516       # delayed load of Carp and avoid recursion
3517       require Carp;
3518       Carp::croak ("$c: Can't call a method without name");
3519       }
3520     if (!_method_hand_up($name))
3521       {
3522       # delayed load of Carp and avoid recursion
3523       require Carp;
3524       Carp::croak ("Can't call $c\-\>$name, not a valid method");
3525       }
3526     # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
3527     $name =~ s/^f/b/;
3528     return &{"Math::BigInt"."::$name"}(@_);
3529     }
3530   my $bname = $name; $bname =~ s/^f/b/;
3531   $c .= "::$name";
3532   *{$c} = \&{$bname};
3533   &{$c};        # uses @_
3534   }
3535
3536 sub exponent
3537   {
3538   # return a copy of the exponent
3539   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3540
3541   if ($x->{sign} !~ /^[+-]$/)
3542     {
3543     my $s = $x->{sign}; $s =~ s/^[+-]//;
3544     return Math::BigInt->new($s);               # -inf, +inf => +inf
3545     }
3546   Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
3547   }
3548
3549 sub mantissa
3550   {
3551   # return a copy of the mantissa
3552   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3553
3554   if ($x->{sign} !~ /^[+-]$/)
3555     {
3556     my $s = $x->{sign}; $s =~ s/^[+]//;
3557     return Math::BigInt->new($s);               # -inf, +inf => +inf
3558     }
3559   my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
3560   $m->bneg() if $x->{sign} eq '-';
3561
3562   $m;
3563   }
3564
3565 sub parts
3566   {
3567   # return a copy of both the exponent and the mantissa
3568   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3569
3570   if ($x->{sign} !~ /^[+-]$/)
3571     {
3572     my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
3573     return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
3574     }
3575   my $m = Math::BigInt->bzero();
3576   $m->{value} = $MBI->_copy($x->{_m});
3577   $m->bneg() if $x->{sign} eq '-';
3578   ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
3579   }
3580
3581 ##############################################################################
3582 # private stuff (internal use only)
3583
3584 sub import
3585   {
3586   my $self = shift;
3587   my $l = scalar @_;
3588   my $lib = ''; my @a;
3589   my $lib_kind = 'try';
3590   $IMPORT=1;
3591   for ( my $i = 0; $i < $l ; $i++)
3592     {
3593     if ( $_[$i] eq ':constant' )
3594       {
3595       # This causes overlord er load to step in. 'binary' and 'integer'
3596       # are handled by BigInt.
3597       overload::constant float => sub { $self->new(shift); };
3598       }
3599     elsif ($_[$i] eq 'upgrade')
3600       {
3601       # this causes upgrading
3602       $upgrade = $_[$i+1];              # or undef to disable
3603       $i++;
3604       }
3605     elsif ($_[$i] eq 'downgrade')
3606       {
3607       # this causes downgrading
3608       $downgrade = $_[$i+1];            # or undef to disable
3609       $i++;
3610       }
3611     elsif ($_[$i] =~ /^(lib|try|only)\z/)
3612       {
3613       # alternative library
3614       $lib = $_[$i+1] || '';            # default Calc
3615       $lib_kind = $1;                   # lib, try or only
3616       $i++;
3617       }
3618     elsif ($_[$i] eq 'with')
3619       {
3620       # alternative class for our private parts()
3621       # XXX: no longer supported
3622       # $MBI = $_[$i+1] || 'Math::BigInt';
3623       $i++;
3624       }
3625     else
3626       {
3627       push @a, $_[$i];
3628       }
3629     }
3630
3631   $lib =~ tr/a-zA-Z0-9,://cd;           # restrict to sane characters
3632   # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
3633   my $mbilib = eval { Math::BigInt->config()->{lib} };
3634   if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
3635     {
3636     # MBI already loaded
3637     Math::BigInt->import( $lib_kind, "$lib,$mbilib", 'objectify');
3638     }
3639   else
3640     {
3641     # MBI not loaded, or with ne "Math::BigInt::Calc"
3642     $lib .= ",$mbilib" if defined $mbilib;
3643     $lib =~ s/^,//;                             # don't leave empty
3644
3645     # replacement library can handle lib statement, but also could ignore it
3646
3647     # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
3648     # used in the same script, or eval inside import(). So we require MBI:
3649     require Math::BigInt;
3650     Math::BigInt->import( $lib_kind => $lib, 'objectify' );
3651     }
3652   if ($@)
3653     {
3654     require Carp; Carp::croak ("Couldn't load $lib: $! $@");
3655     }
3656   # find out which one was actually loaded
3657   $MBI = Math::BigInt->config()->{lib};
3658
3659   # register us with MBI to get notified of future lib changes
3660   Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
3661
3662   $self->export_to_level(1,$self,@a);           # export wanted functions
3663   }
3664
3665 sub bnorm
3666   {
3667   # adjust m and e so that m is smallest possible
3668   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
3669
3670   return $x if $x->{sign} !~ /^[+-]$/;          # inf, nan etc
3671
3672   my $zeros = $MBI->_zeros($x->{_m});           # correct for trailing zeros
3673   if ($zeros != 0)
3674     {
3675     my $z = $MBI->_new($zeros);
3676     $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
3677     if ($x->{_es} eq '-')
3678       {
3679       if ($MBI->_acmp($x->{_e},$z) >= 0)
3680         {
3681         $x->{_e} = $MBI->_sub ($x->{_e}, $z);
3682         $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
3683         }
3684       else
3685         {
3686         $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
3687         $x->{_es} = '+';
3688         }
3689       }
3690     else
3691       {
3692       $x->{_e} = $MBI->_add ($x->{_e}, $z);
3693       }
3694     }
3695   else
3696     {
3697     # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
3698     # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
3699     $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
3700      if $MBI->_is_zero($x->{_m});
3701     }
3702
3703   $x;                                   # MBI bnorm is no-op, so dont call it
3704   }
3705
3706 ##############################################################################
3707
3708 sub as_hex
3709   {
3710   # return number as hexadecimal string (only for integers defined)
3711   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3712
3713   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3714   return '0x0' if $x->is_zero();
3715
3716   return $nan if $x->{_es} ne '+';              # how to do 1e-1 in hex!?
3717
3718   my $z = $MBI->_copy($x->{_m});
3719   if (! $MBI->_is_zero($x->{_e}))               # > 0
3720     {
3721     $MBI->_lsft( $z, $x->{_e},10);
3722     }
3723   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3724   $z->as_hex();
3725   }
3726
3727 sub as_bin
3728   {
3729   # return number as binary digit string (only for integers defined)
3730   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3731
3732   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3733   return '0b0' if $x->is_zero();
3734
3735   return $nan if $x->{_es} ne '+';              # how to do 1e-1 in hex!?
3736
3737   my $z = $MBI->_copy($x->{_m});
3738   if (! $MBI->_is_zero($x->{_e}))               # > 0
3739     {
3740     $MBI->_lsft( $z, $x->{_e},10);
3741     }
3742   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3743   $z->as_bin();
3744   }
3745
3746 sub as_oct
3747   {
3748   # return number as octal digit string (only for integers defined)
3749   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3750
3751   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3752   return '0' if $x->is_zero();
3753
3754   return $nan if $x->{_es} ne '+';              # how to do 1e-1 in hex!?
3755
3756   my $z = $MBI->_copy($x->{_m});
3757   if (! $MBI->_is_zero($x->{_e}))               # > 0
3758     {
3759     $MBI->_lsft( $z, $x->{_e},10);
3760     }
3761   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3762   $z->as_oct();
3763   }
3764
3765 sub as_number
3766   {
3767   # return copy as a bigint representation of this BigFloat number
3768   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3769
3770   return $x if $x->modify('as_number');
3771
3772   if (!$x->isa('Math::BigFloat'))
3773     {
3774     # if the object can as_number(), use it
3775     return $x->as_number() if $x->can('as_number');
3776     # otherwise, get us a float and then a number
3777     $x = $x->can('as_float') ? $x->as_float() : $self->new(0+"$x");
3778     }
3779
3780   return Math::BigInt->binf($x->sign()) if $x->is_inf();
3781   return Math::BigInt->bnan()           if $x->is_nan();
3782
3783   my $z = $MBI->_copy($x->{_m});
3784   if ($x->{_es} eq '-')                 # < 0
3785     {
3786     $MBI->_rsft( $z, $x->{_e},10);
3787     }
3788   elsif (! $MBI->_is_zero($x->{_e}))    # > 0
3789     {
3790     $MBI->_lsft( $z, $x->{_e},10);
3791     }
3792   $z = Math::BigInt->new( $x->{sign} . $MBI->_str($z));
3793   $z;
3794   }
3795
3796 sub length
3797   {
3798   my $x = shift;
3799   my $class = ref($x) || $x;
3800   $x = $class->new(shift) unless ref($x);
3801
3802   return 1 if $MBI->_is_zero($x->{_m});
3803
3804   my $len = $MBI->_len($x->{_m});
3805   $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
3806   if (wantarray())
3807     {
3808     my $t = 0;
3809     $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
3810     return ($len, $t);
3811     }
3812   $len;
3813   }
3814
3815 1;
3816 __END__
3817
3818 =head1 NAME
3819
3820 Math::BigFloat - Arbitrary size floating point math package
3821
3822 =head1 SYNOPSIS
3823
3824   use Math::BigFloat;
3825
3826   # Number creation
3827   my $x = Math::BigFloat->new($str);    # defaults to 0
3828   my $y = $x->copy();                   # make a true copy
3829   my $nan  = Math::BigFloat->bnan();    # create a NotANumber
3830   my $zero = Math::BigFloat->bzero();   # create a +0
3831   my $inf = Math::BigFloat->binf();     # create a +inf
3832   my $inf = Math::BigFloat->binf('-');  # create a -inf
3833   my $one = Math::BigFloat->bone();     # create a +1
3834   my $mone = Math::BigFloat->bone('-'); # create a -1
3835
3836   my $pi = Math::BigFloat->bpi(100);    # PI to 100 digits
3837
3838   # the following examples compute their result to 100 digits accuracy:
3839   my $cos  = Math::BigFloat->new(1)->bcos(100);         # cosinus(1)
3840   my $sin  = Math::BigFloat->new(1)->bsin(100);         # sinus(1)
3841   my $atan = Math::BigFloat->new(1)->batan(100);        # arcus tangens(1)
3842
3843   my $atan2 = Math::BigFloat->new(  1 )->batan2( 1 ,100); # batan(1)
3844   my $atan2 = Math::BigFloat->new(  1 )->batan2( 8 ,100); # batan(1/8)
3845   my $atan2 = Math::BigFloat->new( -2 )->batan2( 1 ,100); # batan(-2)
3846
3847   # Testing
3848   $x->is_zero();                # true if arg is +0
3849   $x->is_nan();                 # true if arg is NaN
3850   $x->is_one();                 # true if arg is +1
3851   $x->is_one('-');              # true if arg is -1
3852   $x->is_odd();                 # true if odd, false for even
3853   $x->is_even();                # true if even, false for odd
3854   $x->is_pos();                 # true if >= 0
3855   $x->is_neg();                 # true if <  0
3856   $x->is_inf(sign);             # true if +inf, or -inf (default is '+')
3857
3858   $x->bcmp($y);                 # compare numbers (undef,<0,=0,>0)
3859   $x->bacmp($y);                # compare absolutely (undef,<0,=0,>0)
3860   $x->sign();                   # return the sign, either +,- or NaN
3861   $x->digit($n);                # return the nth digit, counting from right
3862   $x->digit(-$n);               # return the nth digit, counting from left
3863
3864   # The following all modify their first argument. If you want to preserve
3865   # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
3866   # necessary when mixing $a = $b assignments with non-overloaded math.
3867
3868   # set
3869   $x->bzero();                  # set $i to 0
3870   $x->bnan();                   # set $i to NaN
3871   $x->bone();                   # set $x to +1
3872   $x->bone('-');                # set $x to -1
3873   $x->binf();                   # set $x to inf
3874   $x->binf('-');                # set $x to -inf
3875
3876   $x->bneg();                   # negation
3877   $x->babs();                   # absolute value
3878   $x->bnorm();                  # normalize (no-op)
3879   $x->bnot();                   # two's complement (bit wise not)
3880   $x->binc();                   # increment x by 1
3881   $x->bdec();                   # decrement x by 1
3882
3883   $x->badd($y);                 # addition (add $y to $x)
3884   $x->bsub($y);                 # subtraction (subtract $y from $x)
3885   $x->bmul($y);                 # multiplication (multiply $x by $y)
3886   $x->bdiv($y);                 # divide, set $x to quotient
3887                                 # return (quo,rem) or quo if scalar
3888
3889   $x->bmod($y);                 # modulus ($x % $y)
3890   $x->bpow($y);                 # power of arguments ($x ** $y)
3891   $x->bmodpow($exp,$mod);       # modular exponentiation (($num**$exp) % $mod))
3892   $x->blsft($y, $n);            # left shift by $y places in base $n
3893   $x->brsft($y, $n);            # right shift by $y places in base $n
3894                                 # returns (quo,rem) or quo if in scalar context
3895
3896   $x->blog();                   # logarithm of $x to base e (Euler's number)
3897   $x->blog($base);              # logarithm of $x to base $base (f.i. 2)
3898   $x->bexp();                   # calculate e ** $x where e is Euler's number
3899
3900   $x->band($y);                 # bit-wise and
3901   $x->bior($y);                 # bit-wise inclusive or
3902   $x->bxor($y);                 # bit-wise exclusive or
3903   $x->bnot();                   # bit-wise not (two's complement)
3904
3905   $x->bsqrt();                  # calculate square-root
3906   $x->broot($y);                # $y'th root of $x (e.g. $y == 3 => cubic root)
3907   $x->bfac();                   # factorial of $x (1*2*3*4*..$x)
3908
3909   $x->bround($N);               # accuracy: preserve $N digits
3910   $x->bfround($N);              # precision: round to the $Nth digit
3911
3912   $x->bfloor();                 # return integer less or equal than $x
3913   $x->bceil();                  # return integer greater or equal than $x
3914
3915   # The following do not modify their arguments:
3916
3917   bgcd(@values);                # greatest common divisor
3918   blcm(@values);                # lowest common multiplicator
3919
3920   $x->bstr();                   # return string
3921   $x->bsstr();                  # return string in scientific notation
3922
3923   $x->as_int();                 # return $x as BigInt
3924   $x->exponent();               # return exponent as BigInt
3925   $x->mantissa();               # return mantissa as BigInt
3926   $x->parts();                  # return (mantissa,exponent) as BigInt
3927
3928   $x->length();                 # number of digits (w/o sign and '.')
3929   ($l,$f) = $x->length();       # number of digits, and length of fraction
3930
3931   $x->precision();              # return P of $x (or global, if P of $x undef)
3932   $x->precision($n);            # set P of $x to $n
3933   $x->accuracy();               # return A of $x (or global, if A of $x undef)
3934   $x->accuracy($n);             # set A $x to $n
3935
3936   # these get/set the appropriate global value for all BigFloat objects
3937   Math::BigFloat->precision();  # Precision
3938   Math::BigFloat->accuracy();   # Accuracy
3939   Math::BigFloat->round_mode(); # rounding mode
3940
3941 =head1 DESCRIPTION
3942
3943 All operators (including basic math operations) are overloaded if you
3944 declare your big floating point numbers as
3945
3946   $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
3947
3948 Operations with overloaded operators preserve the arguments, which is
3949 exactly what you expect.
3950
3951 =head2 Canonical notation
3952
3953 Input to these routines are either BigFloat objects, or strings of the
3954 following four forms:
3955
3956 =over 2
3957
3958 =item *
3959
3960 C</^[+-]\d+$/>
3961
3962 =item *
3963
3964 C</^[+-]\d+\.\d*$/>
3965
3966 =item *
3967
3968 C</^[+-]\d+E[+-]?\d+$/>
3969
3970 =item *
3971
3972 C</^[+-]\d*\.\d+E[+-]?\d+$/>
3973
3974 =back
3975
3976 all with optional leading and trailing zeros and/or spaces. Additionally,
3977 numbers are allowed to have an underscore between any two digits.
3978
3979 Empty strings as well as other illegal numbers results in 'NaN'.
3980
3981 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
3982 are always stored in normalized form. On a string, it creates a BigFloat
3983 object.
3984
3985 =head2 Output
3986
3987 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
3988
3989 The string output will always have leading and trailing zeros stripped and drop
3990 a plus sign. C<bstr()> will give you always the form with a decimal point,
3991 while C<bsstr()> (s for scientific) gives you the scientific notation.
3992
3993         Input                   bstr()          bsstr()
3994         '-0'                    '0'             '0E1'
3995         '  -123 123 123'        '-123123123'    '-123123123E0'
3996         '00.0123'               '0.0123'        '123E-4'
3997         '123.45E-2'             '1.2345'        '12345E-4'
3998         '10E+3'                 '10000'         '1E4'
3999
4000 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
4001 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
4002 return either undef, <0, 0 or >0 and are suited for sort.
4003
4004 Actual math is done by using the class defined with C<< with => Class; >> (which
4005 defaults to BigInts) to represent the mantissa and exponent.
4006
4007 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
4008 represent the result when input arguments are not numbers, as well as
4009 the result of dividing by zero.
4010
4011 =head2 C<mantissa()>, C<exponent()> and C<parts()>
4012
4013 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
4014 as BigInts such that:
4015
4016         $m = $x->mantissa();
4017         $e = $x->exponent();
4018         $y = $m * ( 10 ** $e );
4019         print "ok\n" if $x == $y;
4020
4021 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
4022
4023 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
4024
4025 Currently the mantissa is reduced as much as possible, favouring higher
4026 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
4027 This might change in the future, so do not depend on it.
4028
4029 =head2 Accuracy vs. Precision
4030
4031 See also: L<Rounding|Rounding>.
4032
4033 Math::BigFloat supports both precision (rounding to a certain place before or
4034 after the dot) and accuracy (rounding to a certain number of digits). For a
4035 full documentation, examples and tips on these topics please see the large
4036 section about rounding in L<Math::BigInt>.
4037
4038 Since things like C<sqrt(2)> or C<1 / 3> must presented with a limited
4039 accuracy lest a operation consumes all resources, each operation produces
4040 no more than the requested number of digits.
4041
4042 If there is no global precision or accuracy set, B<and> the operation in
4043 question was not called with a requested precision or accuracy, B<and> the
4044 input $x has no accuracy or precision set, then a fallback parameter will
4045 be used. For historical reasons, it is called C<div_scale> and can be accessed
4046 via:
4047
4048         $d = Math::BigFloat->div_scale();               # query
4049         Math::BigFloat->div_scale($n);                  # set to $n digits
4050
4051 The default value for C<div_scale> is 40.
4052
4053 In case the result of one operation has more digits than specified,
4054 it is rounded. The rounding mode taken is either the default mode, or the one
4055 supplied to the operation after the I<scale>:
4056
4057         $x = Math::BigFloat->new(2);
4058         Math::BigFloat->accuracy(5);            # 5 digits max
4059         $y = $x->copy()->bdiv(3);               # will give 0.66667
4060         $y = $x->copy()->bdiv(3,6);             # will give 0.666667
4061         $y = $x->copy()->bdiv(3,6,undef,'odd'); # will give 0.666667
4062         Math::BigFloat->round_mode('zero');
4063         $y = $x->copy()->bdiv(3,6);             # will also give 0.666667
4064
4065 Note that C<< Math::BigFloat->accuracy() >> and C<< Math::BigFloat->precision() >>
4066 set the global variables, and thus B<any> newly created number will be subject
4067 to the global rounding B<immediately>. This means that in the examples above, the
4068 C<3> as argument to C<bdiv()> will also get an accuracy of B<5>.
4069
4070 It is less confusing to either calculate the result fully, and afterwards
4071 round it explicitly, or use the additional parameters to the math
4072 functions like so:
4073
4074         use Math::BigFloat;
4075         $x = Math::BigFloat->new(2);
4076         $y = $x->copy()->bdiv(3);
4077         print $y->bround(5),"\n";               # will give 0.66667
4078
4079         or
4080
4081         use Math::BigFloat;
4082         $x = Math::BigFloat->new(2);
4083         $y = $x->copy()->bdiv(3,5);             # will give 0.66667
4084         print "$y\n";
4085
4086 =head2 Rounding
4087
4088 =over 2
4089
4090 =item ffround ( +$scale )
4091
4092 Rounds to the $scale'th place left from the '.', counting from the dot.
4093 The first digit is numbered 1.
4094
4095 =item ffround ( -$scale )
4096
4097 Rounds to the $scale'th place right from the '.', counting from the dot.
4098
4099 =item ffround ( 0 )
4100
4101 Rounds to an integer.
4102
4103 =item fround  ( +$scale )
4104
4105 Preserves accuracy to $scale digits from the left (aka significant digits)
4106 and pads the rest with zeros. If the number is between 1 and -1, the
4107 significant digits count from the first non-zero after the '.'
4108
4109 =item fround  ( -$scale ) and fround ( 0 )
4110
4111 These are effectively no-ops.
4112
4113 =back
4114
4115 All rounding functions take as a second parameter a rounding mode from one of
4116 the following: 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'.
4117
4118 The default rounding mode is 'even'. By using
4119 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
4120 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
4121 no longer supported.
4122 The second parameter to the round functions then overrides the default
4123 temporarily.
4124
4125 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
4126 'trunc' as rounding mode to make it equivalent to:
4127
4128         $x = 2.5;
4129         $y = int($x) + 2;
4130
4131 You can override this by passing the desired rounding mode as parameter to
4132 C<as_number()>:
4133
4134         $x = Math::BigFloat->new(2.5);
4135         $y = $x->as_number('odd');      # $y = 3
4136
4137 =head1 METHODS
4138
4139 Math::BigFloat supports all methods that Math::BigInt supports, except it
4140 calculates non-integer results when possible. Please see L<Math::BigInt>
4141 for a full description of each method. Below are just the most important
4142 differences:
4143
4144 =head2 accuracy
4145
4146         $x->accuracy(5);                # local for $x
4147         CLASS->accuracy(5);             # global for all members of CLASS
4148                                         # Note: This also applies to new()!
4149
4150         $A = $x->accuracy();            # read out accuracy that affects $x
4151         $A = CLASS->accuracy();         # read out global accuracy
4152
4153 Set or get the global or local accuracy, aka how many significant digits the
4154 results have. If you set a global accuracy, then this also applies to new()!
4155
4156 Warning! The accuracy I<sticks>, e.g. once you created a number under the
4157 influence of C<< CLASS->accuracy($A) >>, all results from math operations with
4158 that number will also be rounded.
4159
4160 In most cases, you should probably round the results explicitly using one of
4161 L<round()>, L<bround()> or L<bfround()> or by passing the desired accuracy
4162 to the math operation as additional parameter:
4163
4164         my $x = Math::BigInt->new(30000);
4165         my $y = Math::BigInt->new(7);
4166         print scalar $x->copy()->bdiv($y, 2);           # print 4300
4167         print scalar $x->copy()->bdiv($y)->bround(2);   # print 4300
4168
4169 =head2 precision()
4170
4171         $x->precision(-2);      # local for $x, round at the second digit right of the dot
4172         $x->precision(2);       # ditto, round at the second digit left of the dot
4173
4174         CLASS->precision(5);    # Global for all members of CLASS
4175                                 # This also applies to new()!
4176         CLASS->precision(-5);   # ditto
4177
4178         $P = CLASS->precision();        # read out global precision
4179         $P = $x->precision();           # read out precision that affects $x
4180
4181 Note: You probably want to use L<accuracy()> instead. With L<accuracy> you
4182 set the number of digits each result should have, with L<precision> you
4183 set the place where to round!
4184
4185 =head2 bexp()
4186
4187         $x->bexp($accuracy);            # calculate e ** X
4188
4189 Calculates the expression C<e ** $x> where C<e> is Euler's number.
4190
4191 This method was added in v1.82 of Math::BigInt (April 2007).
4192
4193 =head2 bnok()
4194
4195         $x->bnok($y);              # x over y (binomial coefficient n over k)
4196
4197 Calculates the binomial coefficient n over k, also called the "choose"
4198 function. The result is equivalent to:
4199
4200         ( n )      n!
4201         | - |  = -------
4202         ( k )    k!(n-k)!
4203
4204 This method was added in v1.84 of Math::BigInt (April 2007).
4205
4206 =head2 bpi()
4207
4208         print Math::BigFloat->bpi(100), "\n";
4209
4210 Calculate PI to N digits (including the 3 before the dot). The result is
4211 rounded according to the current rounding mode, which defaults to "even".
4212
4213 This method was added in v1.87 of Math::BigInt (June 2007).
4214
4215 =head2 bcos()
4216
4217         my $x = Math::BigFloat->new(1);
4218         print $x->bcos(100), "\n";
4219
4220 Calculate the cosinus of $x, modifying $x in place.
4221
4222 This method was added in v1.87 of Math::BigInt (June 2007).
4223
4224 =head2 bsin()
4225
4226         my $x = Math::BigFloat->new(1);
4227         print $x->bsin(100), "\n";
4228
4229 Calculate the sinus of $x, modifying $x in place.
4230
4231 This method was added in v1.87 of Math::BigInt (June 2007).
4232
4233 =head2 batan2()
4234
4235         my $y = Math::BigFloat->new(2);
4236         my $x = Math::BigFloat->new(3);
4237         print $y->batan2($x), "\n";
4238
4239 Calculate the arcus tanges of C<$y> divided by C<$x>, modifying $y in place.
4240 See also L<batan()>.
4241
4242 This method was added in v1.87 of Math::BigInt (June 2007).
4243
4244 =head2 batan()
4245
4246         my $x = Math::BigFloat->new(1);
4247         print $x->batan(100), "\n";
4248
4249 Calculate the arcus tanges of $x, modifying $x in place. See also L<batan2()>.
4250
4251 This method was added in v1.87 of Math::BigInt (June 2007).
4252
4253 =head2 bmuladd()
4254
4255         $x->bmuladd($y,$z);
4256
4257 Multiply $x by $y, and then add $z to the result.
4258
4259 This method was added in v1.87 of Math::BigInt (June 2007).
4260
4261 =head1 Autocreating constants
4262
4263 After C<use Math::BigFloat ':constant'> all the floating point constants
4264 in the given scope are converted to C<Math::BigFloat>. This conversion
4265 happens at compile time.
4266
4267 In particular
4268
4269   perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
4270
4271 prints the value of C<2E-100>. Note that without conversion of
4272 constants the expression 2E-100 will be calculated as normal floating point
4273 number.
4274
4275 Please note that ':constant' does not affect integer constants, nor binary
4276 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
4277 work.
4278
4279 =head2 Math library
4280
4281 Math with the numbers is done (by default) by a module called
4282 Math::BigInt::Calc. This is equivalent to saying:
4283
4284         use Math::BigFloat lib => 'Calc';
4285
4286 You can change this by using:
4287
4288         use Math::BigFloat lib => 'GMP';
4289
4290 B<Note>: General purpose packages should not be explicit about the library
4291 to use; let the script author decide which is best.
4292
4293 Note: The keyword 'lib' will warn when the requested library could not be
4294 loaded. To suppress the warning use 'try' instead:
4295
4296         use Math::BigFloat try => 'GMP';
4297
4298 If your script works with huge numbers and Calc is too slow for them,
4299 you can also for the loading of one of these libraries and if none
4300 of them can be used, the code will die:
4301
4302         use Math::BigFloat only => 'GMP,Pari';
4303
4304 The following would first try to find Math::BigInt::Foo, then
4305 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
4306
4307         use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
4308
4309 See the respective low-level library documentation for further details.
4310
4311 Please note that Math::BigFloat does B<not> use the denoted library itself,
4312 but it merely passes the lib argument to Math::BigInt. So, instead of the need
4313 to do:
4314
4315         use Math::BigInt lib => 'GMP';
4316         use Math::BigFloat;
4317
4318 you can roll it all into one line:
4319
4320         use Math::BigFloat lib => 'GMP';
4321
4322 It is also possible to just require Math::BigFloat:
4323
4324         require Math::BigFloat;
4325
4326 This will load the necessary things (like BigInt) when they are needed, and
4327 automatically.
4328
4329 See L<Math::BigInt> for more details than you ever wanted to know about using
4330 a different low-level library.
4331
4332 =head2 Using Math::BigInt::Lite
4333
4334 For backwards compatibility reasons it is still possible to
4335 request a different storage class for use with Math::BigFloat:
4336
4337         use Math::BigFloat with => 'Math::BigInt::Lite';
4338
4339 However, this request is ignored, as the current code now uses the low-level
4340 math library for directly storing the number parts.
4341
4342 =head1 EXPORTS
4343
4344 C<Math::BigFloat> exports nothing by default, but can export the C<bpi()> method:
4345
4346         use Math::BigFloat qw/bpi/;
4347
4348         print bpi(10), "\n";
4349
4350 =head1 BUGS
4351
4352 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
4353
4354 =head1 CAVEATS
4355
4356 Do not try to be clever to insert some operations in between switching
4357 libraries:
4358
4359         require Math::BigFloat;
4360         my $matter = Math::BigFloat->bone() + 4;        # load BigInt and Calc
4361         Math::BigFloat->import( lib => 'Pari' );        # load Pari, too
4362         my $anti_matter = Math::BigFloat->bone()+4;     # now use Pari
4363
4364 This will create objects with numbers stored in two different backend libraries,
4365 and B<VERY BAD THINGS> will happen when you use these together:
4366
4367         my $flash_and_bang = $matter + $anti_matter;    # Don't do this!
4368
4369 =over 1
4370
4371 =item stringify, bstr()
4372
4373 Both stringify and bstr() now drop the leading '+'. The old code would return
4374 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
4375 reasoning and details.
4376
4377 =item bdiv
4378
4379 The following will probably not print what you expect:
4380
4381         print $c->bdiv(123.456),"\n";
4382
4383 It prints both quotient and remainder since print works in list context. Also,
4384 bdiv() will modify $c, so be careful. You probably want to use
4385
4386         print $c / 123.456,"\n";
4387         print scalar $c->bdiv(123.456),"\n";  # or if you want to modify $c
4388
4389 instead.
4390
4391 =item brsft
4392
4393 The following will probably not print what you expect:
4394
4395         my $c = Math::BigFloat->new('3.14159');
4396         print $c->brsft(3,10),"\n";     # prints 0.00314153.1415
4397
4398 It prints both quotient and remainder, since print calls C<brsft()> in list
4399 context. Also, C<< $c->brsft() >> will modify $c, so be careful.
4400 You probably want to use
4401
4402         print scalar $c->copy()->brsft(3,10),"\n";
4403         # or if you really want to modify $c
4404         print scalar $c->brsft(3,10),"\n";
4405
4406 instead.
4407
4408 =item Modifying and =
4409
4410 Beware of:
4411
4412         $x = Math::BigFloat->new(5);
4413         $y = $x;
4414
4415 It will not do what you think, e.g. making a copy of $x. Instead it just makes
4416 a second reference to the B<same> object and stores it in $y. Thus anything
4417 that modifies $x will modify $y (except overloaded math operators), and vice
4418 versa. See L<Math::BigInt> for details and how to avoid that.
4419
4420 =item bpow
4421
4422 C<bpow()> now modifies the first argument, unlike the old code which left
4423 it alone and only returned the result. This is to be consistent with
4424 C<badd()> etc. The first will modify $x, the second one won't:
4425
4426         print bpow($x,$i),"\n";         # modify $x
4427         print $x->bpow($i),"\n";        # ditto
4428         print $x ** $i,"\n";            # leave $x alone
4429
4430 =item precision() vs. accuracy()
4431
4432 A common pitfall is to use L<precision()> when you want to round a result to
4433 a certain number of digits:
4434
4435         use Math::BigFloat;
4436
4437         Math::BigFloat->precision(4);           # does not do what you think it does
4438         my $x = Math::BigFloat->new(12345);     # rounds $x to "12000"!
4439         print "$x\n";                           # print "12000"
4440         my $y = Math::BigFloat->new(3);         # rounds $y to "0"!
4441         print "$y\n";                           # print "0"
4442         $z = $x / $y;                           # 12000 / 0 => NaN!
4443         print "$z\n";
4444         print $z->precision(),"\n";             # 4
4445
4446 Replacing L<precision> with L<accuracy> is probably not what you want, either:
4447
4448         use Math::BigFloat;
4449
4450         Math::BigFloat->accuracy(4);            # enables global rounding:
4451         my $x = Math::BigFloat->new(123456);    # rounded immediately to "12350"
4452         print "$x\n";                           # print "123500"
4453         my $y = Math::BigFloat->new(3);         # rounded to "3
4454         print "$y\n";                           # print "3"
4455         print $z = $x->copy()->bdiv($y),"\n";   # 41170
4456         print $z->accuracy(),"\n";              # 4
4457
4458 What you want to use instead is:
4459
4460         use Math::BigFloat;
4461
4462         my $x = Math::BigFloat->new(123456);    # no rounding
4463         print "$x\n";                           # print "123456"
4464         my $y = Math::BigFloat->new(3);         # no rounding
4465         print "$y\n";                           # print "3"
4466         print $z = $x->copy()->bdiv($y,4),"\n"; # 41150
4467         print $z->accuracy(),"\n";              # undef
4468
4469 In addition to computing what you expected, the last example also does B<not>
4470 "taint" the result with an accuracy or precision setting, which would
4471 influence any further operation.
4472
4473 =back
4474
4475 =head1 SEE ALSO
4476
4477 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
4478 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and  L<Math::BigInt::GMP>.
4479
4480 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
4481 because they solve the autoupgrading/downgrading issue, at least partly.
4482
4483 The package at L<http://search.cpan.org/~tels/Math-BigInt> contains
4484 more documentation including a full version history, testcases, empty
4485 subclass files and benchmarks.
4486
4487 =head1 LICENSE
4488
4489 This program is free software; you may redistribute it and/or modify it under
4490 the same terms as Perl itself.
4491
4492 =head1 AUTHORS
4493
4494 Mark Biggar, overloaded interface by Ilya Zakharevich.
4495 Completely rewritten by Tels L<http://bloodgate.com> in 2001 - 2006, and still
4496 at it in 2007.
4497
4498 =cut