lib/Math/BigFloat.pm

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