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