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