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