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