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