lib/Math/BigFloat.pm

   1 # The following hash values are internally used:
   2 #   _e: exponent (BigInt)
   3 #   _m: mantissa (absolute BigInt)
   4 # sign: +,-,"NaN" if not a number
   5 #   _a: accuracy
   6 #   _p: precision
   7 #   _f: flags, used to signal MBI not to touch our private parts
   8
   9 package Math::BigFloat;
  10
  11 $VERSION = '1.28';
  12 require 5.005;
  13 use Exporter;
  14 use Math::BigInt qw/objectify/;
  15 @ISA =       qw( Exporter Math::BigInt);
  16
  17 use strict;
  18 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode/;
  19 use vars qw/$upgrade $downgrade/;
  20 my $class = "Math::BigFloat";
  21
  22 use overload
  23 '<=>'   =>      sub { $_[2] ?
  24                       ref($_[0])->bcmp($_[1],$_[0]) :
  25                       ref($_[0])->bcmp($_[0],$_[1])},
  26 'int'   =>      sub { $_[0]->as_number() },             # 'trunc' to bigint
  27 'log'   =>      sub { $_[0]->blog() },
  28 ;
  29
  30 ##############################################################################
  31 # global constants, flags and accessory
  32
  33 use constant MB_NEVER_ROUND => 0x0001;
  34
  35 # are NaNs ok?
  36 my $NaNOK=1;
  37 # constant for easier life
  38 my $nan = 'NaN';
  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
  49 ##############################################################################
  50 # the old code had $rnd_mode, so we need to support it, too
  51
  52 $rnd_mode   = 'even';
  53 sub TIESCALAR   { my ($class) = @_; bless \$round_mode, $class; }
  54 sub FETCH       { return $round_mode; }
  55 sub STORE       { $rnd_mode = $_[0]->round_mode($_[1]); }
  56
  57 BEGIN { tie $rnd_mode, 'Math::BigFloat'; }
  58
  59 ##############################################################################
  60
  61 # in case we call SUPER::->foo() and this wants to call modify()
  62 # sub modify () { 0; }
  63
  64 {
  65   # valid method aliases for AUTOLOAD
  66   my %methods = map { $_ => 1 }
  67    qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
  68         fint facmp fcmp fzero fnan finf finc fdec flog ffac
  69         fceil ffloor frsft flsft fone flog
  70       /;
  71   # valid method's that can be hand-ed up (for AUTOLOAD)
  72   my %hand_ups = map { $_ => 1 }
  73    qw / is_nan is_inf is_negative is_positive
  74         accuracy precision div_scale round_mode fneg fabs babs fnot
  75       /;
  76
  77   sub method_alias { return exists $methods{$_[0]||''}; }
  78   sub method_hand_up { return exists $hand_ups{$_[0]||''}; }
  79 }
  80
  81 ##############################################################################
  82 # constructors
  83
  84 sub new
  85   {
  86   # create a new BigFloat object from a string or another bigfloat object.
  87   # _e: exponent
  88   # _m: mantissa
  89   # sign  => sign (+/-), or "NaN"
  90
  91   my ($class,$wanted,@r) = @_;
  92
  93   # avoid numify-calls by not using || on $wanted!
  94   return $class->bzero() if !defined $wanted;   # default to 0
  95   return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
  96
  97   my $self = {}; bless $self, $class;
  98   # shortcut for bigints and its subclasses
  99   if ((ref($wanted)) && (ref($wanted) ne $class))
 100     {
 101     $self->{_m} = $wanted->as_number();         # get us a bigint copy
 102     $self->{_e} = Math::BigInt->bzero();
 103     $self->{_m}->babs();
 104     $self->{sign} = $wanted->sign();
 105     return $self->bnorm();
 106     }
 107   # got string
 108   # handle '+inf', '-inf' first
 109   if ($wanted =~ /^[+-]?inf$/)
 110     {
 111     $self->{_e} = Math::BigInt->bzero();
 112     $self->{_m} = Math::BigInt->bzero();
 113     $self->{sign} = $wanted;
 114     $self->{sign} = '+inf' if $self->{sign} eq 'inf';
 115     return $self->bnorm();
 116     }
 117   #print "new string '$wanted'\n";
 118   my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
 119   if (!ref $mis)
 120     {
 121     die "$wanted is not a number initialized to $class" if !$NaNOK;
 122     $self->{_e} = Math::BigInt->bzero();
 123     $self->{_m} = Math::BigInt->bzero();
 124     $self->{sign} = $nan;
 125     }
 126   else
 127     {
 128     # make integer from mantissa by adjusting exp, then convert to bigint
 129     # undef,undef to signal MBI that we don't need no bloody rounding
 130     $self->{_e} = Math::BigInt->new("$$es$$ev",undef,undef);    # exponent
 131     $self->{_m} = Math::BigInt->new("$$miv$$mfv",undef,undef);  # create mant.
 132     # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
 133     $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;
 134     $self->{sign} = $$mis;
 135     }
 136   # print "mbf new $self->{sign} $self->{_m} e $self->{_e}\n";
 137   $self->bnorm()->round(@r);            # first normalize, then round
 138   }
 139
 140 sub bnan
 141   {
 142   # create a bigfloat 'NaN', if given a BigFloat, set it to 'NaN'
 143   my $self = shift;
 144   $self = $class if !defined $self;
 145   if (!ref($self))
 146     {
 147     my $c = $self; $self = {}; bless $self, $c;
 148     }
 149   $self->{_m} = Math::BigInt->bzero();
 150   $self->{_e} = Math::BigInt->bzero();
 151   $self->{sign} = $nan;
 152   $self->{_a} = undef; $self->{_p} = undef;
 153   $self;
 154   }
 155
 156 sub binf
 157   {
 158   # create a bigfloat '+-inf', if given a BigFloat, set it to '+-inf'
 159   my $self = shift;
 160   my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
 161
 162   $self = $class if !defined $self;
 163   if (!ref($self))
 164     {
 165     my $c = $self; $self = {}; bless $self, $c;
 166     }
 167   $self->{_m} = Math::BigInt->bzero();
 168   $self->{_e} = Math::BigInt->bzero();
 169   $self->{sign} = $sign.'inf';
 170   $self->{_a} = undef; $self->{_p} = undef;
 171   $self;
 172   }
 173
 174 sub bone
 175   {
 176   # create a bigfloat '+-1', if given a BigFloat, set it to '+-1'
 177   my $self = shift;
 178   my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
 179
 180   $self = $class if !defined $self;
 181   if (!ref($self))
 182     {
 183     my $c = $self; $self = {}; bless $self, $c;
 184     }
 185   $self->{_m} = Math::BigInt->bone();
 186   $self->{_e} = Math::BigInt->bzero();
 187   $self->{sign} = $sign;
 188   if (@_ > 0)
 189     {
 190     $self->{_a} = $_[0]
 191      if (defined $self->{_a} && defined $_[0] && $_[0] > $self->{_a});
 192     $self->{_p} = $_[1]
 193      if (defined $self->{_p} && defined $_[1] && $_[1] < $self->{_p});
 194     }
 195   return $self;
 196   }
 197
 198 sub bzero
 199   {
 200   # create a bigfloat '+0', if given a BigFloat, set it to 0
 201   my $self = shift;
 202   $self = $class if !defined $self;
 203   if (!ref($self))
 204     {
 205     my $c = $self; $self = {}; bless $self, $c;
 206     }
 207   $self->{_m} = Math::BigInt->bzero();
 208   $self->{_e} = Math::BigInt->bone();
 209   $self->{sign} = '+';
 210   if (@_ > 0)
 211     {
 212     $self->{_a} = $_[0]
 213      if (defined $self->{_a} && defined $_[0] && $_[0] > $self->{_a});
 214     $self->{_p} = $_[1]
 215      if (defined $self->{_p} && defined $_[1] && $_[1] < $self->{_p});
 216     }
 217   return $self;
 218   }
 219
 220 ##############################################################################
 221 # string conversation
 222
 223 sub bstr
 224   {
 225   # (ref to BFLOAT or num_str ) return num_str
 226   # Convert number from internal format to (non-scientific) string format.
 227   # internal format is always normalized (no leading zeros, "-0" => "+0")
 228   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 229   #my $x = shift; my $class = ref($x) || $x;
 230   #$x = $class->new(shift) unless ref($x);
 231
 232   #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
 233   #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
 234   if ($x->{sign} !~ /^[+-]$/)
 235     {
 236     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
 237     return 'inf';                                       # +inf
 238     }
 239
 240   my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
 241
 242   my $not_zero = !$x->is_zero();
 243   if ($not_zero)
 244     {
 245     $es = $x->{_m}->bstr();
 246     $len = CORE::length($es);
 247     if (!$x->{_e}->is_zero())
 248       {
 249       if ($x->{_e}->sign() eq '-')
 250         {
 251         $dot = '';
 252         if ($x->{_e} <= -$len)
 253           {
 254           # print "style: 0.xxxx\n";
 255           my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
 256           $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
 257           }
 258         else
 259           {
 260           # print "insert '.' at $x->{_e} in '$es'\n";
 261           substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
 262           }
 263         }
 264       else
 265         {
 266         # expand with zeros
 267         $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
 268         }
 269       }
 270     } # if not zero
 271   $es = $x->{sign}.$es if $x->{sign} eq '-';
 272   # if set accuracy or precision, pad with zeros
 273   if ((defined $x->{_a}) && ($not_zero))
 274     {
 275     # 123400 => 6, 0.1234 => 4, 0.001234 => 4
 276     my $zeros = $x->{_a} - $cad;                # cad == 0 => 12340
 277     $zeros = $x->{_a} - $len if $cad != $len;
 278     $es .= $dot.'0' x $zeros if $zeros > 0;
 279     }
 280   elsif ($x->{_p} || 0 < 0)
 281     {
 282     # 123400 => 6, 0.1234 => 4, 0.001234 => 6
 283     my $zeros = -$x->{_p} + $cad;
 284     $es .= $dot.'0' x $zeros if $zeros > 0;
 285     }
 286   return $es;
 287   }
 288
 289 sub bsstr
 290   {
 291   # (ref to BFLOAT or num_str ) return num_str
 292   # Convert number from internal format to scientific string format.
 293   # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
 294   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 295   #my $x = shift; my $class = ref($x) || $x;
 296   #$x = $class->new(shift) unless ref($x);
 297
 298   #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
 299   #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
 300   if ($x->{sign} !~ /^[+-]$/)
 301     {
 302     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
 303     return 'inf';                                       # +inf
 304     }
 305   my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
 306   my $sep = 'e'.$sign;
 307   return $x->{_m}->bstr().$sep.$x->{_e}->bstr();
 308   }
 309
 310 sub numify
 311   {
 312   # Make a number from a BigFloat object
 313   # simple return string and let Perl's atoi()/atof() handle the rest
 314   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 315   return $x->bsstr();
 316   }
 317
 318 ##############################################################################
 319 # public stuff (usually prefixed with "b")
 320
 321 # tels 2001-08-04
 322 # todo: this must be overwritten and return NaN for non-integer values
 323 # band(), bior(), bxor(), too
 324 #sub bnot
 325 #  {
 326 #  $class->SUPER::bnot($class,@_);
 327 #  }
 328
 329 sub bcmp
 330   {
 331   # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
 332   # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
 333   my ($self,$x,$y) = objectify(2,@_);
 334
 335   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
 336     {
 337     # handle +-inf and NaN
 338     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 339     return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
 340     return +1 if $x->{sign} eq '+inf';
 341     return -1 if $x->{sign} eq '-inf';
 342     return -1 if $y->{sign} eq '+inf';
 343     return +1;
 344     }
 345
 346   # check sign for speed first
 347   return 1 if $x->{sign} eq '+' && $y->{sign} eq '-';   # does also 0 <=> -y
 348   return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';  # does also -x <=> 0
 349
 350   # shortcut
 351   my $xz = $x->is_zero();
 352   my $yz = $y->is_zero();
 353   return 0 if $xz && $yz;                               # 0 <=> 0
 354   return -1 if $xz && $y->{sign} eq '+';                # 0 <=> +y
 355   return 1 if $yz && $x->{sign} eq '+';                 # +x <=> 0
 356
 357   # adjust so that exponents are equal
 358   my $lxm = $x->{_m}->length();
 359   my $lym = $y->{_m}->length();
 360   my $lx = $lxm + $x->{_e};
 361   my $ly = $lym + $y->{_e};
 362   my $l = $lx - $ly; $l->bneg() if $x->{sign} eq '-';
 363   return $l <=> 0 if $l != 0;
 364
 365   # lengths (corrected by exponent) are equal
 366   # so make mantissa euqal length by padding with zero (shift left)
 367   my $diff = $lxm - $lym;
 368   my $xm = $x->{_m};            # not yet copy it
 369   my $ym = $y->{_m};
 370   if ($diff > 0)
 371     {
 372     $ym = $y->{_m}->copy()->blsft($diff,10);
 373     }
 374   elsif ($diff < 0)
 375     {
 376     $xm = $x->{_m}->copy()->blsft(-$diff,10);
 377     }
 378   my $rc = $xm->bcmp($ym);
 379   $rc = -$rc if $x->{sign} eq '-';              # -124 < -123
 380   $rc <=> 0;
 381   }
 382
 383 sub bacmp
 384   {
 385   # Compares 2 values, ignoring their signs.
 386   # Returns one of undef, <0, =0, >0. (suitable for sort)
 387   # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
 388   my ($self,$x,$y) = objectify(2,@_);
 389
 390   # handle +-inf and NaN's
 391   if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
 392     {
 393     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 394     return 0 if ($x->is_inf() && $y->is_inf());
 395     return 1 if ($x->is_inf() && !$y->is_inf());
 396     return -1;
 397     }
 398
 399   # shortcut
 400   my $xz = $x->is_zero();
 401   my $yz = $y->is_zero();
 402   return 0 if $xz && $yz;                               # 0 <=> 0
 403   return -1 if $xz && !$yz;                             # 0 <=> +y
 404   return 1 if $yz && !$xz;                              # +x <=> 0
 405
 406   # adjust so that exponents are equal
 407   my $lxm = $x->{_m}->length();
 408   my $lym = $y->{_m}->length();
 409   my $lx = $lxm + $x->{_e};
 410   my $ly = $lym + $y->{_e};
 411   my $l = $lx - $ly;
 412   return $l <=> 0 if $l != 0;
 413
 414   # lengths (corrected by exponent) are equal
 415   # so make mantissa equal-length by padding with zero (shift left)
 416   my $diff = $lxm - $lym;
 417   my $xm = $x->{_m};            # not yet copy it
 418   my $ym = $y->{_m};
 419   if ($diff > 0)
 420     {
 421     $ym = $y->{_m}->copy()->blsft($diff,10);
 422     }
 423   elsif ($diff < 0)
 424     {
 425     $xm = $x->{_m}->copy()->blsft(-$diff,10);
 426     }
 427   $xm->bcmp($ym) <=> 0;
 428   }
 429
 430 sub badd
 431   {
 432   # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
 433   # return result as BFLOAT
 434   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 435
 436   #print "mbf badd $x $y\n";
 437   # inf and NaN handling
 438   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
 439     {
 440     # NaN first
 441     return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 442     # inf handline
 443     if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
 444       {
 445       # + and + => +, - and - => -, + and - => 0, - and + => 0
 446       return $x->bzero() if $x->{sign} ne $y->{sign};
 447       return $x;
 448       }
 449     # +-inf + something => +inf
 450     # something +-inf => +-inf
 451     $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
 452     return $x;
 453     }
 454
 455   # speed: no add for 0+y or x+0
 456   return $x if $y->is_zero();                           # x+0
 457   if ($x->is_zero())                                    # 0+y
 458     {
 459     # make copy, clobbering up x (modify in place!)
 460     $x->{_e} = $y->{_e}->copy();
 461     $x->{_m} = $y->{_m}->copy();
 462     $x->{sign} = $y->{sign} || $nan;
 463     return $x->round($a,$p,$r,$y);
 464     }
 465
 466   # take lower of the two e's and adapt m1 to it to match m2
 467   my $e = $y->{_e}; $e = Math::BigInt::bzero() if !defined $e;  # if no BFLOAT
 468   $e = $e - $x->{_e};
 469   my $add = $y->{_m}->copy();
 470   if ($e < 0)
 471     {
 472     my $e1 = $e->copy()->babs();
 473     $x->{_m} *= (10 ** $e1);
 474     $x->{_e} += $e;                     # need the sign of e
 475     }
 476   elsif ($e > 0)
 477     {
 478     $add *= (10 ** $e);
 479     }
 480   # else: both e are the same, so just leave them
 481   $x->{_m}->{sign} = $x->{sign};                # fiddle with signs
 482   $add->{sign} = $y->{sign};
 483   $x->{_m} += $add;                             # finally do add/sub
 484   $x->{sign} = $x->{_m}->{sign};                # re-adjust signs
 485   $x->{_m}->{sign} = '+';                       # mantissa always positiv
 486   # delete trailing zeros, then round
 487   return $x->bnorm()->round($a,$p,$r,$y);
 488   }
 489
 490 sub bsub
 491   {
 492   # (BigFloat or num_str, BigFloat or num_str) return BigFloat
 493   # subtract second arg from first, modify first
 494   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 495
 496   if (!$y->is_zero())           # don't need to do anything if $y is 0
 497     {
 498     $y->{sign} =~ tr/+\-/-+/;   # does nothing for NaN
 499     $x->badd($y,$a,$p,$r);      # badd does not leave internal zeros
 500     $y->{sign} =~ tr/+\-/-+/;   # refix $y (does nothing for NaN)
 501     }
 502   $x;                           # already rounded by badd()
 503   }
 504
 505 sub binc
 506   {
 507   # increment arg by one
 508   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 509
 510   if ($x->{_e}->sign() eq '-')
 511     {
 512     return $x->badd($self->bone(),$a,$p,$r);    #  digits after dot
 513     }
 514
 515   if (!$x->{_e}->is_zero())
 516     {
 517     $x->{_m}->blsft($x->{_e},10);               # 1e2 => 100
 518     $x->{_e}->bzero();
 519     }
 520   # now $x->{_e} == 0
 521   if ($x->{sign} eq '+')
 522     {
 523     $x->{_m}->binc();
 524     return $x->bnorm()->bround($a,$p,$r);
 525     }
 526   elsif ($x->{sign} eq '-')
 527     {
 528     $x->{_m}->bdec();
 529     $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
 530     return $x->bnorm()->bround($a,$p,$r);
 531     }
 532   # inf, nan handling etc
 533   $x->badd($self->__one(),$a,$p,$r);            # does round
 534   }
 535
 536 sub bdec
 537   {
 538   # decrement arg by one
 539   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 540
 541   if ($x->{_e}->sign() eq '-')
 542     {
 543     return $x->badd($self->bone('-'),$a,$p,$r); #  digits after dot
 544     }
 545
 546   if (!$x->{_e}->is_zero())
 547     {
 548     $x->{_m}->blsft($x->{_e},10);               # 1e2 => 100
 549     $x->{_e}->bzero();
 550     }
 551   # now $x->{_e} == 0
 552   my $zero = $x->is_zero();
 553   # <= 0
 554   if (($x->{sign} eq '-') || $zero)
 555     {
 556     $x->{_m}->binc();
 557     $x->{sign} = '-' if $zero;                  # 0 => 1 => -1
 558     $x->{sign} = '+' if $x->{_m}->is_zero();    # -1 +1 => -0 => +0
 559     return $x->bnorm()->round($a,$p,$r);
 560     }
 561   # > 0
 562   elsif ($x->{sign} eq '+')
 563     {
 564     $x->{_m}->bdec();
 565     return $x->bnorm()->round($a,$p,$r);
 566     }
 567   # inf, nan handling etc
 568   $x->badd($self->bone('-'),$a,$p,$r);          # does round
 569   }
 570
 571 sub blog
 572   {
 573   my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(2,@_);
 574
 575   # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
 576
 577   # u = x-1, v = x +1
 578   #              _                               _
 579   # taylor:     |    u    1   u^3   1   u^5       |
 580   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
 581   #             |_   v    3   v^3   5   v^5      _|
 582
 583   # we need to limit the accuracy to protect against overflow
 584   my $fallback = 0;
 585   my $scale = 0;
 586   my @params = $x->_find_round_parameters($a,$p,$r);
 587
 588   # no rounding at all, so must use fallback
 589   if (scalar @params == 1)
 590     {
 591     # simulate old behaviour
 592     $params[1] = $self->div_scale();    # and round to it as accuracy
 593     $scale = $params[1]+4;              # at least four more for proper round
 594     $params[3] = $r;                    # round mode by caller or undef
 595     $fallback = 1;                      # to clear a/p afterwards
 596     }
 597   else
 598     {
 599     # the 4 below is empirical, and there might be cases where it is not
 600     # enough...
 601     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
 602     }
 603
 604   return $x->bzero(@params) if $x->is_one();
 605   return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
 606   #return $x->bone('+',@params) if $x->bcmp($base) == 0;
 607
 608   # when user set globals, they would interfere with our calculation, so
 609   # disable then and later re-enable them
 610   no strict 'refs';
 611   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
 612   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
 613   # we also need to disable any set A or P on $x (_find_round_parameters took
 614   # them already into account), since these would interfere, too
 615   delete $x->{_a}; delete $x->{_p};
 616   # need to disable $upgrade in BigInt, to aoid deep recursion
 617   local $Math::BigInt::upgrade = undef;
 618
 619   my $v = $x->copy(); $v->binc();               # v = x+1
 620   $x->bdec(); my $u = $x->copy();               # u = x-1; x = x-1
 621
 622   $x->bdiv($v,$scale);                          # first term: u/v
 623
 624   my $below = $v->copy();
 625   my $over = $u->copy();
 626   $u *= $u; $v *= $v;                           # u^2, v^2
 627   $below->bmul($v);                             # u^3, v^3
 628   $over->bmul($u);
 629   my $factor = $self->new(3); my $two = $self->new(2);
 630
 631   my $diff = $self->bone();
 632   my $limit = $self->new("1E-". ($scale-1)); my $last;
 633   # print "diff $diff limit $limit\n";
 634   while ($diff->bcmp($limit) > 0)
 635     {
 636     #print "$x $over $below $factor\n";
 637     $diff = $x->copy()->bsub($last)->babs();
 638     #print "diff $diff $limit\n";
 639     $last = $x->copy();
 640     $x += $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
 641     $over *= $u; $below *= $v; $factor->badd($two);
 642     }
 643   $x->bmul($two);
 644
 645   # shortcut to not run trough _find_round_parameters again
 646   if (defined $params[1])
 647     {
 648     $x->bround($params[1],$params[3]);          # then round accordingly
 649     }
 650   else
 651     {
 652     $x->bfround($params[2],$params[3]);         # then round accordingly
 653     }
 654   if ($fallback)
 655     {
 656     # clear a/p after round, since user did not request it
 657     $x->{_a} = undef; $x->{_p} = undef;
 658     }
 659   # restore globals
 660   $$abr = $ab; $$pbr = $pb;
 661
 662   $x;
 663   }
 664
 665 sub blcm
 666   {
 667   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
 668   # does not modify arguments, but returns new object
 669   # Lowest Common Multiplicator
 670
 671   my ($self,@arg) = objectify(0,@_);
 672   my $x = $self->new(shift @arg);
 673   while (@arg) { $x = _lcm($x,shift @arg); }
 674   $x;
 675   }
 676
 677 sub bgcd
 678   {
 679   # (BFLOAT or num_str, BFLOAT or num_str) return BINT
 680   # does not modify arguments, but returns new object
 681   # GCD -- Euclids algorithm Knuth Vol 2 pg 296
 682
 683   my ($self,@arg) = objectify(0,@_);
 684   my $x = $self->new(shift @arg);
 685   while (@arg) { $x = _gcd($x,shift @arg); }
 686   $x;
 687   }
 688
 689 ###############################################################################
 690 # is_foo methods (is_negative, is_positive are inherited from BigInt)
 691
 692 sub is_int
 693   {
 694   # return true if arg (BFLOAT or num_str) is an integer
 695   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 696
 697   return 1 if ($x->{sign} =~ /^[+-]$/) &&       # NaN and +-inf aren't
 698     $x->{_e}->{sign} eq '+';                    # 1e-1 => no integer
 699   0;
 700   }
 701
 702 sub is_zero
 703   {
 704   # return true if arg (BFLOAT or num_str) is zero
 705   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 706
 707   return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
 708   0;
 709   }
 710
 711 sub is_one
 712   {
 713   # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
 714   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 715
 716   my $sign = shift || ''; $sign = '+' if $sign ne '-';
 717   return 1
 718    if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
 719   0;
 720   }
 721
 722 sub is_odd
 723   {
 724   # return true if arg (BFLOAT or num_str) is odd or false if even
 725   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 726
 727   return 1 if ($x->{sign} =~ /^[+-]$/) &&               # NaN & +-inf aren't
 728     ($x->{_e}->is_zero() && $x->{_m}->is_odd());
 729   0;
 730   }
 731
 732 sub is_even
 733   {
 734   # return true if arg (BINT or num_str) is even or false if odd
 735   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 736
 737   return 0 if $x->{sign} !~ /^[+-]$/;                   # NaN & +-inf aren't
 738 #  return 1 if $x->{_m}->is_zero();                     # 0e1 is even
 739   return 1 if ($x->{_e}->{sign} eq '+'                  # 123.45 is never
 740      && $x->{_m}->is_even());                           # but 1200 is
 741   0;
 742   }
 743
 744 sub bmul
 745   {
 746   # multiply two numbers -- stolen from Knuth Vol 2 pg 233
 747   # (BINT or num_str, BINT or num_str) return BINT
 748   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 749
 750   # print "mbf bmul $x->{_m}e$x->{_e} $y->{_m}e$y->{_e}\n";
 751   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 752
 753   # handle result = 0
 754   return $x->bzero() if $x->is_zero() || $y->is_zero();
 755   # inf handling
 756   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
 757     {
 758     # result will always be +-inf:
 759     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
 760     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
 761     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
 762     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
 763     return $x->binf('-');
 764     }
 765
 766   # aEb * cEd = (a*c)E(b+d)
 767   $x->{_m}->bmul($y->{_m});
 768   $x->{_e}->badd($y->{_e});
 769   # adjust sign:
 770   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
 771   return $x->bnorm()->round($a,$p,$r,$y);
 772   }
 773
 774 sub bdiv
 775   {
 776   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
 777   # (BFLOAT,BFLOAT) (quo,rem) or BINT (only rem)
 778   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 779
 780   # x / +-inf => 0, reminder x
 781   return wantarray ? ($x->bzero(),$x->copy()) : $x->bzero()
 782    if $y->{sign} =~ /^[+-]inf$/;
 783
 784   # NaN if x == NaN or y == NaN or x==y==0
 785   return wantarray ? ($x->bnan(),bnan()) : $x->bnan()
 786    if (($x->is_nan() || $y->is_nan()) ||
 787       ($x->is_zero() && $y->is_zero()));
 788
 789   # 5 / 0 => +inf, -6 / 0 => -inf
 790   return wantarray
 791    ? ($x->binf($x->{sign}),$self->bnan()) : $x->binf($x->{sign})
 792    if ($x->{sign} =~ /^[+-]$/ && $y->is_zero());
 793
 794   # x== 0 or y == 1 or y == -1
 795   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
 796
 797   # we need to limit the accuracy to protect against overflow
 798   my $fallback = 0;
 799   my $scale = 0;
 800   my @params = $x->_find_round_parameters($a,$p,$r,$y);
 801
 802   # no rounding at all, so must use fallback
 803   if (scalar @params == 1)
 804     {
 805     # simulate old behaviour
 806     $params[1] = $self->div_scale();    # and round to it as accuracy
 807     $scale = $params[1]+4;              # at least four more for proper round
 808     $params[3] = $r;                    # round mode by caller or undef
 809     $fallback = 1;                      # to clear a/p afterwards
 810     }
 811   else
 812     {
 813     # the 4 below is empirical, and there might be cases where it is not
 814     # enough...
 815     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
 816     }
 817   my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
 818   $scale = $lx if $lx > $scale;
 819   $scale = $ly if $ly > $scale;
 820   my $diff = $ly - $lx;
 821   $scale += $diff if $diff > 0;         # if lx << ly, but not if ly << lx!
 822
 823   # make copy of $x in case of list context for later reminder calculation
 824   my $rem;
 825   if (wantarray && !$y->is_one())
 826     {
 827     $rem = $x->copy();
 828     }
 829
 830   $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
 831
 832   # check for / +-1 ( +/- 1E0)
 833   if (!$y->is_one())
 834     {
 835     # promote BigInts and it's subclasses (except when already a BigFloat)
 836     $y = $self->new($y) unless $y->isa('Math::BigFloat');
 837
 838     # calculate the result to $scale digits and then round it
 839     # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
 840     $x->{_m}->blsft($scale,10);
 841     $x->{_m}->bdiv( $y->{_m} ); # a/c
 842     $x->{_e}->bsub( $y->{_e} ); # b-d
 843     $x->{_e}->bsub($scale);     # correct for 10**scale
 844     $x->bnorm();                # remove trailing 0's
 845     }
 846
 847   # shortcut to not run trough _find_round_parameters again
 848   if (defined $params[1])
 849     {
 850     $x->bround($params[1],$params[3]);          # then round accordingly
 851     }
 852   else
 853     {
 854     $x->bfround($params[2],$params[3]);         # then round accordingly
 855     }
 856   if ($fallback)
 857     {
 858     # clear a/p after round, since user did not request it
 859     $x->{_a} = undef; $x->{_p} = undef;
 860     }
 861
 862   if (wantarray)
 863     {
 864     if (!$y->is_one())
 865       {
 866       $rem->bmod($y,$params[1],$params[2],$params[3]);  # copy already done
 867       }
 868     else
 869       {
 870       $rem = $self->bzero();
 871       }
 872     if ($fallback)
 873       {
 874       # clear a/p after round, since user did not request it
 875       $rem->{_a} = undef; $rem->{_p} = undef;
 876       }
 877     return ($x,$rem);
 878     }
 879   return $x;
 880   }
 881
 882 sub bmod
 883   {
 884   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
 885   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 886
 887   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
 888     {
 889     my ($d,$re) = $self->SUPER::_div_inf($x,$y);
 890     return $re->round($a,$p,$r,$y);
 891     }
 892   return $x->bnan() if $x->is_zero() && $y->is_zero();
 893   return $x if $y->is_zero();
 894   return $x->bnan() if $x->is_nan() || $y->is_nan();
 895   return $x->bzero() if $y->is_one() || $x->is_zero();
 896
 897   # inf handling is missing here
 898
 899   my $cmp = $x->bacmp($y);                      # equal or $x < $y?
 900   return $x->bzero($a,$p) if $cmp == 0;         # $x == $y => result 0
 901
 902   # only $y of the operands negative?
 903   my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
 904
 905   $x->{sign} = $y->{sign};                              # calc sign first
 906   return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0;  # $x < $y => result $x
 907
 908   my $ym = $y->{_m}->copy();
 909
 910   # 2e1 => 20
 911   $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
 912
 913   # if $y has digits after dot
 914   my $shifty = 0;                       # correct _e of $x by this
 915   if ($y->{_e}->{sign} eq '-')          # has digits after dot
 916     {
 917     # 123 % 2.5 => 1230 % 25 => 5 => 0.5
 918     $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
 919     $x->blsft($shifty,10);              # 123 => 1230, $y->{_m} is already 25
 920     }
 921   # $ym is now mantissa of $y based on exponent 0
 922
 923   my $shiftx = 0;                       # correct _e of $x by this
 924   if ($x->{_e}->{sign} eq '-')          # has digits after dot
 925     {
 926     # 123.4 % 20 => 1234 % 200
 927     $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
 928     $ym->blsft($shiftx,10);
 929     }
 930   # 123e1 % 20 => 1230 % 20
 931   if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
 932     {
 933     $x->{_m}->blsft($x->{_e},10);
 934     }
 935   $x->{_e} = Math::BigInt->bzero() unless $x->{_e}->is_zero();
 936
 937   $x->{_e}->bsub($shiftx) if $shiftx != 0;
 938   $x->{_e}->bsub($shifty) if $shifty != 0;
 939
 940   # now mantissas are equalized, exponent of $x is adjusted, so calc result
 941 #  $ym->{sign} = '-' if $neg; # bmod() will make the correction for us
 942
 943   $x->{_m}->bmod($ym);
 944
 945   $x->{sign} = '+' if $x->{_m}->is_zero();              # fix sign for -0
 946   $x->bnorm();
 947
 948   if ($neg != 0)        # one of them negative => correct in place
 949     {
 950     my $r = $y - $x;
 951     $x->{_m} = $r->{_m};
 952     $x->{_e} = $r->{_e};
 953     $x->{sign} = '+' if $x->{_m}->is_zero();            # fix sign for -0
 954     $x->bnorm();
 955     }
 956
 957   $x->round($a,$p,$r,$y);       # round and return
 958   }
 959
 960 sub bsqrt
 961   {
 962   # calculate square root; this should probably
 963   # use a different test to see whether the accuracy we want is...
 964   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 965
 966   return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
 967   return $x if $x->{sign} eq '+inf';                              # +inf
 968   return $x if $x->is_zero() || $x->is_one();
 969
 970   # we need to limit the accuracy to protect against overflow
 971   my $fallback = 0;
 972   my $scale = 0;
 973   my @params = $x->_find_round_parameters($a,$p,$r);
 974
 975   # no rounding at all, so must use fallback
 976   if (scalar @params == 1)
 977     {
 978     # simulate old behaviour
 979     $params[1] = $self->div_scale();    # and round to it as accuracy
 980     $scale = $params[1]+4;              # at least four more for proper round
 981     $params[3] = $r;                    # round mode by caller or undef
 982     $fallback = 1;                      # to clear a/p afterwards
 983     }
 984   else
 985     {
 986     # the 4 below is empirical, and there might be cases where it is not
 987     # enough...
 988     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
 989     }
 990
 991   # when user set globals, they would interfere with our calculation, so
 992   # disable then and later re-enable them
 993   no strict 'refs';
 994   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
 995   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
 996   # we also need to disable any set A or P on $x (_find_round_parameters took
 997   # them already into account), since these would interfere, too
 998   delete $x->{_a}; delete $x->{_p};
 999   # need to disable $upgrade in BigInt, to aoid deep recursion
1000   local $Math::BigInt::upgrade = undef;
1001
1002   my $xas = $x->as_number();
1003   my $gs = $xas->copy()->bsqrt();       # some guess
1004
1005   if (($x->{_e}->{sign} ne '-')         # guess can't be accurate if there are
1006                                         # digits after the dot
1007    && ($xas->bacmp($gs * $gs) == 0))    # guess hit the nail on the head?
1008     {
1009     # exact result
1010     $x->{_m} = $gs; $x->{_e} = Math::BigInt->bzero(); $x->bnorm();
1011     # shortcut to not run trough _find_round_parameters again
1012     if (defined $params[1])
1013       {
1014       $x->bround($params[1],$params[3]);        # then round accordingly
1015       }
1016     else
1017       {
1018       $x->bfround($params[2],$params[3]);       # then round accordingly
1019       }
1020     if ($fallback)
1021       {
1022       # clear a/p after round, since user did not request it
1023       $x->{_a} = undef; $x->{_p} = undef;
1024       }
1025     ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1026     return $x;
1027     }
1028   $gs = $self->new( $gs );              # BigInt to BigFloat
1029
1030   my $lx = $x->{_m}->length();
1031   $scale = $lx if $scale < $lx;
1032   my $e = $self->new("1E-$scale");      # make test variable
1033 #  return $x->bnan() if $e->sign() eq 'NaN';
1034
1035   my $y = $x->copy();
1036   my $two = $self->new(2);
1037   my $diff = $e;
1038   # promote BigInts and it's subclasses (except when already a BigFloat)
1039   $y = $self->new($y) unless $y->isa('Math::BigFloat');
1040
1041   my $rem;
1042   while ($diff >= $e)
1043     {
1044     $rem = $y->copy()->bdiv($gs,$scale)->badd($gs)->bdiv($two,$scale);
1045     $diff = $rem->copy()->bsub($gs)->babs();
1046     $gs = $rem->copy();
1047     }
1048   # copy over to modify $x
1049   $x->{_m} = $rem->{_m}; $x->{_e} = $rem->{_e};
1050
1051   # shortcut to not run trough _find_round_parameters again
1052   if (defined $params[1])
1053     {
1054     $x->bround($params[1],$params[3]);          # then round accordingly
1055     }
1056   else
1057     {
1058     $x->bfround($params[2],$params[3]);         # then round accordingly
1059     }
1060   if ($fallback)
1061     {
1062     # clear a/p after round, since user did not request it
1063     $x->{_a} = undef; $x->{_p} = undef;
1064     }
1065   # restore globals
1066   $$abr = $ab; $$pbr = $pb;
1067   $x;
1068   }
1069
1070 sub bfac
1071   {
1072   # (BINT or num_str, BINT or num_str) return BINT
1073   # compute factorial numbers
1074   # modifies first argument
1075   my ($self,$x,@r) = objectify(1,@_);
1076
1077   return $x->bnan() if $x->{sign} ne '+';       # inf, NnN, <0 etc => NaN
1078   return $x->bone(@r) if $x->is_zero() || $x->is_one();         # 0 or 1 => 1
1079
1080   return $x->bnan() if $x->{_e}->{sign} ne '+'; # digits after dot?
1081
1082   # use BigInt's bfac() for faster calc
1083   $x->{_m}->blsft($x->{_e},10);         # un-norm m
1084   $x->{_e}->bzero();                    # norm $x again
1085   $x->{_m}->bfac();                     # factorial
1086   $x->bnorm();
1087
1088   #my $n = $x->copy();
1089   #$x->bone();
1090   #my $f = $self->new(2);
1091   #while ($f->bacmp($n) < 0)
1092   #  {
1093   #  $x->bmul($f); $f->binc();
1094   #  }
1095   #$x->bmul($f);                                        # last step
1096   $x->round(@r);                                # round
1097   }
1098
1099 sub bpow
1100   {
1101   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1102   # compute power of two numbers, second arg is used as integer
1103   # modifies first argument
1104
1105   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1106
1107   return $x if $x->{sign} =~ /^[+-]inf$/;
1108   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1109   return $x->bone() if $y->is_zero();
1110   return $x         if $x->is_one() || $y->is_one();
1111   my $y1 = $y->as_number();             # make bigint (trunc)
1112   # if ($x == -1)
1113   if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1114     {
1115     # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
1116     return $y1->is_odd() ? $x : $x->babs(1);
1117     }
1118   return $x if $x->is_zero() && $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1119   # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1120   return $x->binf() if $x->is_zero() && $y->{sign} eq '-';
1121
1122   # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1123   $y1->babs();
1124   $x->{_m}->bpow($y1);
1125   $x->{_e}->bmul($y1);
1126   $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1127   $x->bnorm();
1128   if ($y->{sign} eq '-')
1129     {
1130     # modify $x in place!
1131     my $z = $x->copy(); $x->bzero()->binc();
1132     return $x->bdiv($z,$a,$p,$r);       # round in one go (might ignore y's A!)
1133     }
1134   return $x->round($a,$p,$r,$y);
1135   }
1136
1137 ###############################################################################
1138 # rounding functions
1139
1140 sub bfround
1141   {
1142   # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1143   # $n == 0 means round to integer
1144   # expects and returns normalized numbers!
1145   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1146
1147   return $x if $x->modify('bfround');
1148
1149   my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1150   return $x if !defined $scale;                 # no-op
1151
1152   # never round a 0, +-inf, NaN
1153   if ($x->is_zero())
1154     {
1155     $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1156     return $x;
1157     }
1158   return $x if $x->{sign} !~ /^[+-]$/;
1159   # print "MBF bfround $x to scale $scale mode $mode\n";
1160
1161   # don't round if x already has lower precision
1162   return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1163
1164   $x->{_p} = $scale;                    # remember round in any case
1165   $x->{_a} = undef;                     # and clear A
1166   if ($scale < 0)
1167     {
1168     # print "bfround scale $scale e $x->{_e}\n";
1169     # round right from the '.'
1170     return $x if $x->{_e} >= 0;                 # nothing to round
1171     $scale = -$scale;                           # positive for simplicity
1172     my $len = $x->{_m}->length();               # length of mantissa
1173     my $dad = -$x->{_e};                        # digits after dot
1174     my $zad = 0;                                # zeros after dot
1175     $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
1176     #print "scale $scale dad $dad zad $zad len $len\n";
1177
1178     # number  bsstr   len zad dad
1179     # 0.123   123e-3    3   0 3
1180     # 0.0123  123e-4    3   1 4
1181     # 0.001   1e-3      1   2 3
1182     # 1.23    123e-2    3   0 2
1183     # 1.2345  12345e-4  5   0 4
1184
1185     # do not round after/right of the $dad
1186     return $x if $scale > $dad;                 # 0.123, scale >= 3 => exit
1187
1188     # round to zero if rounding inside the $zad, but not for last zero like:
1189     # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1190     return $x->bzero() if $scale < $zad;
1191     if ($scale == $zad)                 # for 0.006, scale -3 and trunc
1192       {
1193       $scale = -$len;
1194       }
1195     else
1196       {
1197       # adjust round-point to be inside mantissa
1198       if ($zad != 0)
1199         {
1200         $scale = $scale-$zad;
1201         }
1202       else
1203         {
1204         my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;    # digits before dot
1205         $scale = $dbd+$scale;
1206         }
1207       }
1208     # print "round to $x->{_m} to $scale\n";
1209     }
1210   else
1211     {
1212     # 123 => 100 means length(123) = 3 - $scale (2) => 1
1213
1214     my $dbt = $x->{_m}->length();
1215     # digits before dot
1216     my $dbd = $dbt + $x->{_e};
1217     # should be the same, so treat it as this
1218     $scale = 1 if $scale == 0;
1219     # shortcut if already integer
1220     return $x if $scale == 1 && $dbt <= $dbd;
1221     # maximum digits before dot
1222     ++$dbd;
1223
1224     if ($scale > $dbd)
1225        {
1226        # not enough digits before dot, so round to zero
1227        return $x->bzero;
1228        }
1229     elsif ( $scale == $dbd )
1230        {
1231        # maximum
1232        $scale = -$dbt;
1233        }
1234     else
1235        {
1236        $scale = $dbd - $scale;
1237        }
1238
1239     }
1240   # print "using $scale for $x->{_m} with '$mode'\n";
1241   # pass sign to bround for rounding modes '+inf' and '-inf'
1242   $x->{_m}->{sign} = $x->{sign};
1243   $x->{_m}->bround($scale,$mode);
1244   $x->{_m}->{sign} = '+';               # fix sign back
1245   $x->bnorm();
1246   }
1247
1248 sub bround
1249   {
1250   # accuracy: preserve $N digits, and overwrite the rest with 0's
1251   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1252
1253   die ('bround() needs positive accuracy') if ($_[0] || 0) < 0;
1254
1255   my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1256   return $x if !defined $scale;                         # no-op
1257
1258   return $x if $x->modify('bround');
1259
1260   # scale is now either $x->{_a}, $accuracy, or the user parameter
1261   # test whether $x already has lower accuracy, do nothing in this case
1262   # but do round if the accuracy is the same, since a math operation might
1263   # want to round a number with A=5 to 5 digits afterwards again
1264   return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1265
1266   # scale < 0 makes no sense
1267   # never round a +-inf, NaN
1268   return $x if ($scale < 0) ||  $x->{sign} !~ /^[+-]$/;
1269
1270   # 1: $scale == 0 => keep all digits
1271   # 2: never round a 0
1272   # 3: if we should keep more digits than the mantissa has, do nothing
1273   if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1274     {
1275     $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1276     return $x;
1277     }
1278
1279   # pass sign to bround for '+inf' and '-inf' rounding modes
1280   $x->{_m}->{sign} = $x->{sign};
1281   $x->{_m}->bround($scale,$mode);       # round mantissa
1282   $x->{_m}->{sign} = '+';               # fix sign back
1283   # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1284   $x->{_a} = $scale;                    # remember rounding
1285   $x->{_p} = undef;                     # and clear P
1286   $x->bnorm();                          # del trailing zeros gen. by bround()
1287   }
1288
1289 sub bfloor
1290   {
1291   # return integer less or equal then $x
1292   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1293
1294   return $x if $x->modify('bfloor');
1295
1296   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1297
1298   # if $x has digits after dot
1299   if ($x->{_e}->{sign} eq '-')
1300     {
1301     $x->{_m}->brsft(-$x->{_e},10);
1302     $x->{_e}->bzero();
1303     $x-- if $x->{sign} eq '-';
1304     }
1305   $x->round($a,$p,$r);
1306   }
1307
1308 sub bceil
1309   {
1310   # return integer greater or equal then $x
1311   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1312
1313   return $x if $x->modify('bceil');
1314   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1315
1316   # if $x has digits after dot
1317   if ($x->{_e}->{sign} eq '-')
1318     {
1319     $x->{_m}->brsft(-$x->{_e},10);
1320     $x->{_e}->bzero();
1321     $x++ if $x->{sign} eq '+';
1322     }
1323   $x->round($a,$p,$r);
1324   }
1325
1326 sub brsft
1327   {
1328   # shift right by $y (divide by power of 2)
1329   my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1330
1331   return $x if $x->modify('brsft');
1332   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1333
1334   $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1335   $x->bdiv($n ** $y,$a,$p,$r,$y);
1336   }
1337
1338 sub blsft
1339   {
1340   # shift right by $y (divide by power of 2)
1341   my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1342
1343   return $x if $x->modify('brsft');
1344   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1345
1346   $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1347   $x->bmul($n ** $y,$a,$p,$r,$y);
1348   }
1349
1350 ###############################################################################
1351
1352 sub DESTROY
1353   {
1354   # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1355   }
1356
1357 sub AUTOLOAD
1358   {
1359   # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
1360   # or falling back to MBI::bxxx()
1361   my $name = $AUTOLOAD;
1362
1363   $name =~ s/.*:://;    # split package
1364   no strict 'refs';
1365   if (!method_alias($name))
1366     {
1367     if (!defined $name)
1368       {
1369       # delayed load of Carp and avoid recursion
1370       require Carp;
1371       Carp::croak ("Can't call a method without name");
1372       }
1373     if (!method_hand_up($name))
1374       {
1375       # delayed load of Carp and avoid recursion
1376       require Carp;
1377       Carp::croak ("Can't call $class\-\>$name, not a valid method");
1378       }
1379     # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1380     $name =~ s/^f/b/;
1381     return &{'Math::BigInt'."::$name"}(@_);
1382     }
1383   my $bname = $name; $bname =~ s/^f/b/;
1384   *{$class."::$name"} = \&$bname;
1385   &$bname;      # uses @_
1386   }
1387
1388 sub exponent
1389   {
1390   # return a copy of the exponent
1391   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1392
1393   if ($x->{sign} !~ /^[+-]$/)
1394     {
1395     my $s = $x->{sign}; $s =~ s/^[+-]//;
1396     return $self->new($s);                      # -inf, +inf => +inf
1397     }
1398   return $x->{_e}->copy();
1399   }
1400
1401 sub mantissa
1402   {
1403   # return a copy of the mantissa
1404   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1405
1406   if ($x->{sign} !~ /^[+-]$/)
1407     {
1408     my $s = $x->{sign}; $s =~ s/^[+]//;
1409     return $self->new($s);                      # -inf, +inf => +inf
1410     }
1411   my $m = $x->{_m}->copy();             # faster than going via bstr()
1412   $m->bneg() if $x->{sign} eq '-';
1413
1414   $m;
1415   }
1416
1417 sub parts
1418   {
1419   # return a copy of both the exponent and the mantissa
1420   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1421
1422   if ($x->{sign} !~ /^[+-]$/)
1423     {
1424     my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
1425     return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
1426     }
1427   my $m = $x->{_m}->copy();     # faster than going via bstr()
1428   $m->bneg() if $x->{sign} eq '-';
1429   return ($m,$x->{_e}->copy());
1430   }
1431
1432 ##############################################################################
1433 # private stuff (internal use only)
1434
1435 sub import
1436   {
1437   my $self = shift;
1438   my $l = scalar @_; my $j = 0; my @a = @_;
1439   for ( my $i = 0; $i < $l ; $i++, $j++)
1440     {
1441     if ( $_[$i] eq ':constant' )
1442       {
1443       # this rest causes overlord er load to step in
1444       # print "overload @_\n";
1445       overload::constant float => sub { $self->new(shift); };
1446       splice @a, $j, 1; $j--;
1447       }
1448     elsif ($_[$i] eq 'upgrade')
1449       {
1450       # this causes upgrading
1451       $upgrade = $_[$i+1];              # or undef to disable
1452       my $s = 2; $s = 1 if @a-$j < 2;   # avoid "can not modify non-existant..."
1453       splice @a, $j, $s; $j -= $s;
1454       }
1455     }
1456   # any non :constant stuff is handled by our parent, Exporter
1457   # even if @_ is empty, to give it a chance
1458   $self->SUPER::import(@a);             # for subclasses
1459   $self->export_to_level(1,$self,@a);   # need this, too
1460   }
1461
1462 sub bnorm
1463   {
1464   # adjust m and e so that m is smallest possible
1465   # round number according to accuracy and precision settings
1466   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1467
1468   return $x if $x->{sign} !~ /^[+-]$/;          # inf, nan etc
1469
1470   my $zeros = $x->{_m}->_trailing_zeros();      # correct for trailing zeros
1471   if ($zeros != 0)
1472     {
1473     $x->{_m}->brsft($zeros,10); $x->{_e} += $zeros;
1474     }
1475   # for something like 0Ey, set y to 1, and -0 => +0
1476   $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
1477   # this is to prevent automatically rounding when MBI's globals are set
1478   $x->{_m}->{_f} = MB_NEVER_ROUND;
1479   $x->{_e}->{_f} = MB_NEVER_ROUND;
1480   # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
1481   $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
1482   $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
1483   $x;                                   # MBI bnorm is no-op, so dont call it
1484   }
1485
1486 ##############################################################################
1487 # internal calculation routines
1488
1489 sub as_number
1490   {
1491   # return copy as a bigint representation of this BigFloat number
1492   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1493
1494   my $z;
1495   if ($x->{_e}->is_zero())
1496     {
1497     $z = $x->{_m}->copy();
1498     $z->{sign} = $x->{sign};
1499     return $z;
1500     }
1501   $z = $x->{_m}->copy();
1502   if ($x->{_e} < 0)
1503     {
1504     $z->brsft(-$x->{_e},10);
1505     }
1506   else
1507     {
1508     $z->blsft($x->{_e},10);
1509     }
1510   $z->{sign} = $x->{sign};
1511   $z;
1512   }
1513
1514 sub length
1515   {
1516   my $x = shift;
1517   my $class = ref($x) || $x;
1518   $x = $class->new(shift) unless ref($x);
1519
1520   return 1 if $x->{_m}->is_zero();
1521   my $len = $x->{_m}->length();
1522   $len += $x->{_e} if $x->{_e}->sign() eq '+';
1523   if (wantarray())
1524     {
1525     my $t = Math::BigInt::bzero();
1526     $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1527     return ($len,$t);
1528     }
1529   $len;
1530   }
1531
1532 1;
1533 __END__
1534
1535 =head1 NAME
1536
1537 Math::BigFloat - Arbitrary size floating point math package
1538
1539 =head1 SYNOPSIS
1540
1541   use Math::BigFloat;
1542
1543   # Number creation
1544   $x = Math::BigFloat->new($str);       # defaults to 0
1545   $nan  = Math::BigFloat->bnan();       # create a NotANumber
1546   $zero = Math::BigFloat->bzero();      # create a +0
1547   $inf = Math::BigFloat->binf();        # create a +inf
1548   $inf = Math::BigFloat->binf('-');     # create a -inf
1549   $one = Math::BigFloat->bone();        # create a +1
1550   $one = Math::BigFloat->bone('-');     # create a -1
1551
1552   # Testing
1553   $x->is_zero();                # true if arg is +0
1554   $x->is_nan();                 # true if arg is NaN
1555   $x->is_one();                 # true if arg is +1
1556   $x->is_one('-');              # true if arg is -1
1557   $x->is_odd();                 # true if odd, false for even
1558   $x->is_even();                # true if even, false for odd
1559   $x->is_positive();            # true if >= 0
1560   $x->is_negative();            # true if <  0
1561   $x->is_inf(sign);             # true if +inf, or -inf (default is '+')
1562
1563   $x->bcmp($y);                 # compare numbers (undef,<0,=0,>0)
1564   $x->bacmp($y);                # compare absolutely (undef,<0,=0,>0)
1565   $x->sign();                   # return the sign, either +,- or NaN
1566   $x->digit($n);                # return the nth digit, counting from right
1567   $x->digit(-$n);               # return the nth digit, counting from left
1568
1569   # The following all modify their first argument:
1570
1571   # set
1572   $x->bzero();                  # set $i to 0
1573   $x->bnan();                   # set $i to NaN
1574   $x->bone();                   # set $x to +1
1575   $x->bone('-');                # set $x to -1
1576   $x->binf();                   # set $x to inf
1577   $x->binf('-');                # set $x to -inf
1578
1579   $x->bneg();                   # negation
1580   $x->babs();                   # absolute value
1581   $x->bnorm();                  # normalize (no-op)
1582   $x->bnot();                   # two's complement (bit wise not)
1583   $x->binc();                   # increment x by 1
1584   $x->bdec();                   # decrement x by 1
1585
1586   $x->badd($y);                 # addition (add $y to $x)
1587   $x->bsub($y);                 # subtraction (subtract $y from $x)
1588   $x->bmul($y);                 # multiplication (multiply $x by $y)
1589   $x->bdiv($y);                 # divide, set $i to quotient
1590                                 # return (quo,rem) or quo if scalar
1591
1592   $x->bmod($y);                 # modulus
1593   $x->bpow($y);                 # power of arguments (a**b)
1594   $x->blsft($y);                # left shift
1595   $x->brsft($y);                # right shift
1596                                 # return (quo,rem) or quo if scalar
1597
1598   $x->blog($base);              # logarithm of $x, base defaults to e
1599                                 # (other bases than e not supported yet)
1600
1601   $x->band($y);                 # bit-wise and
1602   $x->bior($y);                 # bit-wise inclusive or
1603   $x->bxor($y);                 # bit-wise exclusive or
1604   $x->bnot();                   # bit-wise not (two's complement)
1605
1606   $x->bsqrt();                  # calculate square-root
1607   $x->bfac();                   # factorial of $x (1*2*3*4*..$x)
1608
1609   $x->bround($N);               # accuracy: preserver $N digits
1610   $x->bfround($N);              # precision: round to the $Nth digit
1611
1612   # The following do not modify their arguments:
1613   bgcd(@values);                # greatest common divisor
1614   blcm(@values);                # lowest common multiplicator
1615
1616   $x->bstr();                   # return string
1617   $x->bsstr();                  # return string in scientific notation
1618
1619   $x->bfloor();                 # return integer less or equal than $x
1620   $x->bceil();                  # return integer greater or equal than $x
1621
1622   $x->exponent();               # return exponent as BigInt
1623   $x->mantissa();               # return mantissa as BigInt
1624   $x->parts();                  # return (mantissa,exponent) as BigInt
1625
1626   $x->length();                 # number of digits (w/o sign and '.')
1627   ($l,$f) = $x->length();       # number of digits, and length of fraction
1628
1629 =head1 DESCRIPTION
1630
1631 All operators (inlcuding basic math operations) are overloaded if you
1632 declare your big floating point numbers as
1633
1634   $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1635
1636 Operations with overloaded operators preserve the arguments, which is
1637 exactly what you expect.
1638
1639 =head2 Canonical notation
1640
1641 Input to these routines are either BigFloat objects, or strings of the
1642 following four forms:
1643
1644 =over 2
1645
1646 =item *
1647
1648 C</^[+-]\d+$/>
1649
1650 =item *
1651
1652 C</^[+-]\d+\.\d*$/>
1653
1654 =item *
1655
1656 C</^[+-]\d+E[+-]?\d+$/>
1657
1658 =item *
1659
1660 C</^[+-]\d*\.\d+E[+-]?\d+$/>
1661
1662 =back
1663
1664 all with optional leading and trailing zeros and/or spaces. Additonally,
1665 numbers are allowed to have an underscore between any two digits.
1666
1667 Empty strings as well as other illegal numbers results in 'NaN'.
1668
1669 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1670 are always stored in normalized form. On a string, it creates a BigFloat
1671 object.
1672
1673 =head2 Output
1674
1675 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1676
1677 The string output will always have leading and trailing zeros stripped and drop
1678 a plus sign. C<bstr()> will give you always the form with a decimal point,
1679 while C<bsstr()> (for scientific) gives you the scientific notation.
1680
1681         Input                   bstr()          bsstr()
1682         '-0'                    '0'             '0E1'
1683         '  -123 123 123'        '-123123123'    '-123123123E0'
1684         '00.0123'               '0.0123'        '123E-4'
1685         '123.45E-2'             '1.2345'        '12345E-4'
1686         '10E+3'                 '10000'         '1E4'
1687
1688 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1689 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1690 return either undef, <0, 0 or >0 and are suited for sort.
1691
1692 Actual math is done by using BigInts to represent the mantissa and exponent.
1693 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1694 represent the result when input arguments are not numbers, as well as
1695 the result of dividing by zero.
1696
1697 =head2 C<mantissa()>, C<exponent()> and C<parts()>
1698
1699 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1700 as BigInts such that:
1701
1702         $m = $x->mantissa();
1703         $e = $x->exponent();
1704         $y = $m * ( 10 ** $e );
1705         print "ok\n" if $x == $y;
1706
1707 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
1708
1709 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
1710
1711 Currently the mantissa is reduced as much as possible, favouring higher
1712 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
1713 This might change in the future, so do not depend on it.
1714
1715 =head2 Accuracy vs. Precision
1716
1717 See also: L<Rounding|Rounding>.
1718
1719 Math::BigFloat supports both precision and accuracy. For a full documentation,
1720 examples and tips on these topics please see the large section in
1721 L<Math::BigInt>.
1722
1723 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
1724 a operation consumes all resources, each operation produces no more than
1725 C<Math::BigFloat::precision()> digits.
1726
1727 In case the result of one operation has more precision than specified,
1728 it is rounded. The rounding mode taken is either the default mode, or the one
1729 supplied to the operation after the I<scale>:
1730
1731         $x = Math::BigFloat->new(2);
1732         Math::BigFloat::precision(5);           # 5 digits max
1733         $y = $x->copy()->bdiv(3);               # will give 0.66666
1734         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
1735         $y = $x->copy()->bdiv(3,6,'odd');       # will give 0.666667
1736         Math::BigFloat::round_mode('zero');
1737         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
1738
1739 =head2 Rounding
1740
1741 =over 2
1742
1743 =item ffround ( +$scale )
1744
1745 Rounds to the $scale'th place left from the '.', counting from the dot.
1746 The first digit is numbered 1.
1747
1748 =item ffround ( -$scale )
1749
1750 Rounds to the $scale'th place right from the '.', counting from the dot.
1751
1752 =item ffround ( 0 )
1753
1754 Rounds to an integer.
1755
1756 =item fround  ( +$scale )
1757
1758 Preserves accuracy to $scale digits from the left (aka significant digits)
1759 and pads the rest with zeros. If the number is between 1 and -1, the
1760 significant digits count from the first non-zero after the '.'
1761
1762 =item fround  ( -$scale ) and fround ( 0 )
1763
1764 These are effetively no-ops.
1765
1766 =back
1767
1768 All rounding functions take as a second parameter a rounding mode from one of
1769 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
1770
1771 The default rounding mode is 'even'. By using
1772 C<< Math::BigFloat::round_mode($round_mode); >> you can get and set the default
1773 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
1774 no longer supported.
1775 The second parameter to the round functions then overrides the default
1776 temporarily.
1777
1778 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
1779 'trunc' as rounding mode to make it equivalent to:
1780
1781         $x = 2.5;
1782         $y = int($x) + 2;
1783
1784 You can override this by passing the desired rounding mode as parameter to
1785 C<as_number()>:
1786
1787         $x = Math::BigFloat->new(2.5);
1788         $y = $x->as_number('odd');      # $y = 3
1789
1790 =head1 EXAMPLES
1791
1792   # not ready yet
1793
1794 =head1 Autocreating constants
1795
1796 After C<use Math::BigFloat ':constant'> all the floating point constants
1797 in the given scope are converted to C<Math::BigFloat>. This conversion
1798 happens at compile time.
1799
1800 In particular
1801
1802   perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
1803
1804 prints the value of C<2E-100>.  Note that without conversion of
1805 constants the expression 2E-100 will be calculated as normal floating point
1806 number.
1807
1808 =head1 BUGS
1809
1810 =over 2
1811
1812 =item *
1813
1814 The following does not work yet:
1815
1816         $m = $x->mantissa();
1817         $e = $x->exponent();
1818         $y = $m * ( 10 ** $e );
1819         print "ok\n" if $x == $y;
1820
1821 =item *
1822
1823 There is no fmod() function yet.
1824
1825 =back
1826
1827 =head1 CAVEAT
1828
1829 =over 1
1830
1831 =item stringify, bstr()
1832
1833 Both stringify and bstr() now drop the leading '+'. The old code would return
1834 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
1835 reasoning and details.
1836
1837 =item bdiv
1838
1839 The following will probably not do what you expect:
1840
1841         print $c->bdiv(123.456),"\n";
1842
1843 It prints both quotient and reminder since print works in list context. Also,
1844 bdiv() will modify $c, so be carefull. You probably want to use
1845
1846         print $c / 123.456,"\n";
1847         print scalar $c->bdiv(123.456),"\n";  # or if you want to modify $c
1848
1849 instead.
1850
1851 =item Modifying and =
1852
1853 Beware of:
1854
1855         $x = Math::BigFloat->new(5);
1856         $y = $x;
1857
1858 It will not do what you think, e.g. making a copy of $x. Instead it just makes
1859 a second reference to the B<same> object and stores it in $y. Thus anything
1860 that modifies $x will modify $y, and vice versa.
1861
1862         $x->bmul(2);
1863         print "$x, $y\n";       # prints '10, 10'
1864
1865 If you want a true copy of $x, use:
1866
1867         $y = $x->copy();
1868
1869 See also the documentation in L<overload> regarding C<=>.
1870
1871 =item bpow
1872
1873 C<bpow()> now modifies the first argument, unlike the old code which left
1874 it alone and only returned the result. This is to be consistent with
1875 C<badd()> etc. The first will modify $x, the second one won't:
1876
1877         print bpow($x,$i),"\n";         # modify $x
1878         print $x->bpow($i),"\n";        # ditto
1879         print $x ** $i,"\n";            # leave $x alone
1880
1881 =back
1882
1883 =head1 LICENSE
1884
1885 This program is free software; you may redistribute it and/or modify it under
1886 the same terms as Perl itself.
1887
1888 =head1 AUTHORS
1889
1890 Mark Biggar, overloaded interface by Ilya Zakharevich.
1891 Completely rewritten by Tels http://bloodgate.com in 2001.
1892
1893 =cut