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