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