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