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
 716   # upgrade $x if $x is not a BigFloat (handle BigInt input)
 717   if (!$x->isa('Math::BigFloat'))
 718     {
 719     $x = Math::BigFloat->new($x);
 720     $self = ref($x);
 721     }
 722   # first calculate the log to base e (using reduction by 10 (and probably 2))
 723   $self->_log_10($x,$scale);
 724
 725   # and if a different base was requested, convert it
 726   if (defined $base)
 727     {
 728     # not ln, but some other base
 729     $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
 730     }
 731
 732   # shortcut to not run trough _find_round_parameters again
 733   if (defined $params[0])
 734     {
 735     $x->bround($params[0],$params[2]);          # then round accordingly
 736     }
 737   else
 738     {
 739     $x->bfround($params[1],$params[2]);         # then round accordingly
 740     }
 741   if ($fallback)
 742     {
 743     # clear a/p after round, since user did not request it
 744     $x->{_a} = undef; $x->{_p} = undef;
 745     }
 746   # restore globals
 747   $$abr = $ab; $$pbr = $pb;
 748
 749   $x;
 750   }
 751
 752 sub _log
 753   {
 754   # internal log function to calculate log based on Taylor.
 755   # Modifies $x in place.
 756   my ($self,$x,$scale) = @_;
 757
 758   # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
 759
 760   # u = x-1, v = x+1
 761   #              _                               _
 762   # Taylor:     |    u    1   u^3   1   u^5       |
 763   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
 764   #             |_   v    3   v^3   5   v^5      _|
 765
 766   # This takes much more steps to calculate the result and is thus not used
 767   # u = x-1
 768   #              _                               _
 769   # Taylor:     |    u    1   u^2   1   u^3       |
 770   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 1/2
 771   #             |_   x    2   x^2   3   x^3      _|
 772
 773   # "normal" log algorithmn
 774
 775   my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
 776
 777   $v = $x->copy(); $v->binc();          # v = x+1
 778   $x->bdec(); $u = $x->copy();          # u = x-1; x = x-1
 779   $x->bdiv($v,$scale);                  # first term: u/v
 780   $below = $v->copy();
 781   $over = $u->copy();
 782   $u *= $u; $v *= $v;                           # u^2, v^2
 783   $below->bmul($v);                             # u^3, v^3
 784   $over->bmul($u);
 785   $factor = $self->new(3); $f = $self->new(2);
 786
 787   my $steps = 0 if DEBUG;
 788   $limit = $self->new("1E-". ($scale-1));
 789   while (3 < 5)
 790     {
 791     # we calculate the next term, and add it to the last
 792     # when the next term is below our limit, it won't affect the outcome
 793     # anymore, so we stop
 794
 795     # calculating the next term simple from over/below will result in quite
 796     # a time hog if the input has many digits, since over and below will
 797     # accumulate more and more digits, and the result will also have many
 798     # digits, but in the end it is rounded to $scale digits anyway. So if we
 799     # round $over and $below first, we save a lot of time for the division
 800     # (not with log(1.2345), but try log (123**123) to see what I mean. This
 801     # can introduce a rounding error if the division result would be f.i.
 802     # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
 803     # if we truncated the $over and $below we might get 0.12345. Does this
 804     # matter for the end result? So we give over and below 4 more digits to be
 805     # on the safe side (unscientific error handling as usual...)
 806     # Makes blog(1.23) *slightly* slower, but try blog(123*123) w/o it :o)
 807
 808     $next = $over->copy->bround($scale+4)->bdiv(
 809       $below->copy->bmul($factor)->bround($scale+4),
 810       $scale);
 811
 812 ## old version:
 813 ##    $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
 814
 815     last if $next->bacmp($limit) <= 0;
 816
 817     delete $next->{_a}; delete $next->{_p};
 818     $x->badd($next);
 819     #print "step  $x\n  ($next - $limit = ",$next - $limit,")\n";
 820     # calculate things for the next term
 821     $over *= $u; $below *= $v; $factor->badd($f);
 822     if (DEBUG)
 823       {
 824       $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
 825       }
 826     }
 827   $x->bmul($f);                                 # $x *= 2
 828   print "took $steps steps\n" if DEBUG;
 829   }
 830
 831 sub _log_10
 832   {
 833   # internal log function based on reducing input to the range of 0.1 .. 9.99
 834   my ($self,$x,$scale) = @_;
 835
 836   # taking blog() from numbers greater than 10 takes a *very long* time, so we
 837   # break the computation down into parts based on the observation that:
 838   #  blog(x*y) = blog(x) + blog(y)
 839   # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
 840   # the faster it get's, especially because 2*$x takes about 10 times as long,
 841   # so by dividing $x by 10 we make it at least factor 100 faster...)
 842
 843   # The same observation is valid for numbers smaller than 0.1 (e.g. computing
 844   # log(1) is fastest, and the farther away we get from 1, the longer it takes)
 845   # so we also 'break' this down by multiplying $x with 10 and subtract the
 846   # log(10) afterwards to get the correct result.
 847
 848   # calculate nr of digits before dot
 849   my $dbd = $x->{_m}->length() + $x->{_e}->numify();
 850
 851   # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
 852   # infinite recursion
 853
 854   my $calc = 1;                                 # do some calculation?
 855
 856   # disable the shortcut for 10, since we need log(10) and this would recurse
 857   # infinitely deep
 858   if ($x->{_e}->is_one() && $x->{_m}->is_one())
 859     {
 860     $dbd = 0;                                   # disable shortcut
 861     # we can use the cached value in these cases
 862     if ($scale <= $LOG_10_A)
 863       {
 864       $x->bzero(); $x->badd($LOG_10);
 865       $calc = 0;                                # no need to calc, but round
 866       }
 867     }
 868   # disable the shortcut for 2, since we maybe have it cached
 869   my $two = $self->new(2);                      # also used later on
 870   if ($x->{_e}->is_zero() && $x->{_m}->bcmp($two) == 0)
 871     {
 872     $dbd = 0;                                   # disable shortcut
 873     # we can use the cached value in these cases
 874     if ($scale <= $LOG_2_A)
 875       {
 876       $x->bzero(); $x->badd($LOG_2);
 877       $calc = 0;                                # no need to calc, but round
 878       }
 879     }
 880
 881   # if $x = 0.1, we know the result must be 0-log(10)
 882   if ($x->{_e}->is_one('-') && $x->{_m}->is_one())
 883     {
 884     $dbd = 0;                                   # disable shortcut
 885     # we can use the cached value in these cases
 886     if ($scale <= $LOG_10_A)
 887       {
 888       $x->bzero(); $x->bsub($LOG_10);
 889       $calc = 0;                                # no need to calc, but round
 890       }
 891     }
 892
 893   # default: these correction factors are undef and thus not used
 894   my $l_10;                             # value of ln(10) to A of $scale
 895   my $l_2;                              # value of ln(2) to A of $scale
 896
 897   # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
 898   # so don't do this shortcut for 1 or 0
 899   if (($dbd > 1) || ($dbd < 0))
 900     {
 901     # convert our cached value to an object if not already (avoid doing this
 902     # at import() time, since not everybody needs this)
 903     $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
 904
 905     #print "x = $x, dbd = $dbd, calc = $calc\n";
 906     # got more than one digit before the dot, or more than one zero after the
 907     # dot, so do:
 908     #  log(123)    == log(1.23) + log(10) * 2
 909     #  log(0.0123) == log(1.23) - log(10) * 2
 910
 911     if ($scale <= $LOG_10_A)
 912       {
 913       # use cached value
 914       #print "using cached value for l_10\n";
 915       $l_10 = $LOG_10->copy();          # copy for mul
 916       }
 917     else
 918       {
 919       # else: slower, compute it (but don't cache it, because it could be big)
 920       # also disable downgrade for this code path
 921       local $Math::BigFloat::downgrade = undef;
 922       #print "l_10 = $l_10 (self = $self',
 923       #  ", ref(l_10) = ",ref($l_10)," scale $scale)\n";
 924       #print "calculating value for l_10, scale $scale\n";
 925       $l_10 = $self->new(10)->blog(undef,$scale);       # scale+4, actually
 926       }
 927     $dbd-- if ($dbd > 1);               # 20 => dbd=2, so make it dbd=1
 928     # make object
 929     $dbd = $self->new($dbd);
 930     #print "dbd $dbd\n";
 931     $l_10->bmul($dbd);                  # log(10) * (digits_before_dot-1)
 932     #print "l_10 = $l_10\n";
 933     #print "x = $x";
 934     $x->{_e}->bsub($dbd);               # 123 => 1.23
 935     #print " => $x\n";
 936     #print "calculating log($x) with scale=$scale\n";
 937
 938     }
 939
 940   # Now: 0.1 <= $x < 10 (and possible correction in l_10)
 941
 942   ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
 943   ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
 944
 945   if ($calc != 0)
 946     {
 947     my $half = $self->new('0.5');
 948     my $twos = 0;                               # default: none (0 times)
 949     while ($x->bacmp($half) < 0)
 950       {
 951       #print "$x\n";
 952       $twos--; $x->bmul($two);
 953       }
 954     while ($x->bacmp($two) > 0)
 955       {
 956       #print "$x\n";
 957       $twos++; $x->bdiv($two,$scale+4);         # keep all digits
 958       }
 959     #print "$twos\n";
 960     # $twos > 0 => did mul 2, < 0 => did div 2 (never both)
 961     # calculate correction factor based on ln(2)
 962     if ($twos != 0)
 963       {
 964       $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
 965       if ($scale <= $LOG_2_A)
 966         {
 967         # use cached value
 968         #print "using cached value for l_10\n";
 969         $l_2 = $LOG_2->copy();                  # copy for mul
 970         }
 971       else
 972         {
 973         # else: slower, compute it (but don't cache it, because it could be big)
 974         # also disable downgrade for this code path
 975         local $Math::BigFloat::downgrade = undef;
 976         #print "calculating value for l_2, scale $scale\n";
 977         $l_2 = $two->blog(undef,$scale);        # scale+4, actually
 978         }
 979       #print "$l_2 => \n";
 980       $l_2->bmul($twos);                # * -2 => subtract, * 2 => add
 981       #print "$l_2\n";
 982       }
 983     }
 984
 985   if ($calc != 0)
 986     {
 987     $self->_log($x,$scale);                     # need to do the "normal" way
 988     #print "log(x) = $x\n";
 989     $x->badd($l_10) if defined $l_10;           # correct it by ln(10)
 990     #print "result = $x\n";
 991     $x->badd($l_2) if defined $l_2;             # and maybe by ln(2)
 992     #print "result = $x\n";
 993     }
 994   # all done, $x contains now the result
 995   }
 996
 997 sub blcm
 998   {
 999   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1000   # does not modify arguments, but returns new object
1001   # Lowest Common Multiplicator
1002
1003   my ($self,@arg) = objectify(0,@_);
1004   my $x = $self->new(shift @arg);
1005   while (@arg) { $x = _lcm($x,shift @arg); }
1006   $x;
1007   }
1008
1009 sub bgcd
1010   {
1011   # (BFLOAT or num_str, BFLOAT or num_str) return BINT
1012   # does not modify arguments, but returns new object
1013   # GCD -- Euclids algorithm Knuth Vol 2 pg 296
1014
1015   my ($self,@arg) = objectify(0,@_);
1016   my $x = $self->new(shift @arg);
1017   while (@arg) { $x = _gcd($x,shift @arg); }
1018   $x;
1019   }
1020
1021 ###############################################################################
1022 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1023
1024 sub is_int
1025   {
1026   # return true if arg (BFLOAT or num_str) is an integer
1027   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1028
1029   return 1 if ($x->{sign} =~ /^[+-]$/) &&       # NaN and +-inf aren't
1030     $x->{_e}->{sign} eq '+';                    # 1e-1 => no integer
1031   0;
1032   }
1033
1034 sub is_zero
1035   {
1036   # return true if arg (BFLOAT or num_str) is zero
1037   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1038
1039   return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
1040   0;
1041   }
1042
1043 sub is_one
1044   {
1045   # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1046   my ($self,$x,$sign) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1047
1048   $sign = '+' if !defined $sign || $sign ne '-';
1049   return 1
1050    if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
1051   0;
1052   }
1053
1054 sub is_odd
1055   {
1056   # return true if arg (BFLOAT or num_str) is odd or false if even
1057   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1058
1059   return 1 if ($x->{sign} =~ /^[+-]$/) &&               # NaN & +-inf aren't
1060     ($x->{_e}->is_zero() && $x->{_m}->is_odd());
1061   0;
1062   }
1063
1064 sub is_even
1065   {
1066   # return true if arg (BINT or num_str) is even or false if odd
1067   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1068
1069   return 0 if $x->{sign} !~ /^[+-]$/;                   # NaN & +-inf aren't
1070   return 1 if ($x->{_e}->{sign} eq '+'                  # 123.45 is never
1071      && $x->{_m}->is_even());                           # but 1200 is
1072   0;
1073   }
1074
1075 sub bmul
1076   {
1077   # multiply two numbers -- stolen from Knuth Vol 2 pg 233
1078   # (BINT or num_str, BINT or num_str) return BINT
1079
1080   # set up parameters
1081   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1082   # objectify is costly, so avoid it
1083   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1084     {
1085     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1086     }
1087
1088   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1089
1090   # inf handling
1091   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1092     {
1093     return $x->bnan() if $x->is_zero() || $y->is_zero();
1094     # result will always be +-inf:
1095     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1096     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1097     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1098     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1099     return $x->binf('-');
1100     }
1101   # handle result = 0
1102   return $x->bzero() if $x->is_zero() || $y->is_zero();
1103
1104   return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
1105    ((!$x->isa($self)) || (!$y->isa($self)));
1106
1107   # aEb * cEd = (a*c)E(b+d)
1108   $x->{_m}->bmul($y->{_m});
1109   $x->{_e}->badd($y->{_e});
1110   # adjust sign:
1111   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1112   return $x->bnorm()->round($a,$p,$r,$y);
1113   }
1114
1115 sub bdiv
1116   {
1117   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1118   # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1119
1120   # set up parameters
1121   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1122   # objectify is costly, so avoid it
1123   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1124     {
1125     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1126     }
1127
1128   return $self->_div_inf($x,$y)
1129    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1130
1131   # x== 0 # also: or y == 1 or y == -1
1132   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1133
1134   # upgrade ?
1135   return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1136
1137   # we need to limit the accuracy to protect against overflow
1138   my $fallback = 0;
1139   my (@params,$scale);
1140   ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1141
1142   return $x if $x->is_nan();            # error in _find_round_parameters?
1143
1144   # no rounding at all, so must use fallback
1145   if (scalar @params == 0)
1146     {
1147     # simulate old behaviour
1148     $params[0] = $self->div_scale();    # and round to it as accuracy
1149     $scale = $params[0]+4;              # at least four more for proper round
1150     $params[2] = $r;                    # round mode by caller or undef
1151     $fallback = 1;                      # to clear a/p afterwards
1152     }
1153   else
1154     {
1155     # the 4 below is empirical, and there might be cases where it is not
1156     # enough...
1157     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1158     }
1159   my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
1160   $scale = $lx if $lx > $scale;
1161   $scale = $ly if $ly > $scale;
1162   my $diff = $ly - $lx;
1163   $scale += $diff if $diff > 0;         # if lx << ly, but not if ly << lx!
1164
1165   # make copy of $x in case of list context for later reminder calculation
1166   my $rem;
1167   if (wantarray && !$y->is_one())
1168     {
1169     $rem = $x->copy();
1170     }
1171
1172   $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1173
1174   # check for / +-1 ( +/- 1E0)
1175   if (!$y->is_one())
1176     {
1177     # promote BigInts and it's subclasses (except when already a BigFloat)
1178     $y = $self->new($y) unless $y->isa('Math::BigFloat');
1179
1180     # need to disable $upgrade in BigInt, to avoid deep recursion
1181     local $Math::BigInt::upgrade = undef;       # should be parent class vs MBI
1182
1183     # calculate the result to $scale digits and then round it
1184     # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1185     $x->{_m}->blsft($scale,10);
1186     $x->{_m}->bdiv( $y->{_m} ); # a/c
1187     $x->{_e}->bsub( $y->{_e} ); # b-d
1188     $x->{_e}->bsub($scale);     # correct for 10**scale
1189     $x->bnorm();                # remove trailing 0's
1190     }
1191
1192   # shortcut to not run trough _find_round_parameters again
1193   if (defined $params[0])
1194     {
1195     $x->{_a} = undef;                           # clear before round
1196     $x->bround($params[0],$params[2]);          # then round accordingly
1197     }
1198   else
1199     {
1200     $x->{_p} = undef;                           # clear before round
1201     $x->bfround($params[1],$params[2]);         # then round accordingly
1202     }
1203   if ($fallback)
1204     {
1205     # clear a/p after round, since user did not request it
1206     $x->{_a} = undef; $x->{_p} = undef;
1207     }
1208
1209   if (wantarray)
1210     {
1211     if (!$y->is_one())
1212       {
1213       $rem->bmod($y,@params);                   # copy already done
1214       }
1215     else
1216       {
1217       $rem = $self->bzero();
1218       }
1219     if ($fallback)
1220       {
1221       # clear a/p after round, since user did not request it
1222       $rem->{_a} = undef; $rem->{_p} = undef;
1223       }
1224     return ($x,$rem);
1225     }
1226   $x;
1227   }
1228
1229 sub bmod
1230   {
1231   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1232
1233   # set up parameters
1234   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1235   # objectify is costly, so avoid it
1236   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1237     {
1238     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1239     }
1240
1241   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1242     {
1243     my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1244     $x->{sign} = $re->{sign};
1245     $x->{_e} = $re->{_e};
1246     $x->{_m} = $re->{_m};
1247     return $x->round($a,$p,$r,$y);
1248     }
1249   return $x->bnan() if $x->is_zero() && $y->is_zero();
1250   return $x if $y->is_zero();
1251   return $x->bnan() if $x->is_nan() || $y->is_nan();
1252   return $x->bzero() if $y->is_one() || $x->is_zero();
1253
1254   # inf handling is missing here
1255
1256   my $cmp = $x->bacmp($y);                      # equal or $x < $y?
1257   return $x->bzero($a,$p) if $cmp == 0;         # $x == $y => result 0
1258
1259   # only $y of the operands negative?
1260   my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1261
1262   $x->{sign} = $y->{sign};                              # calc sign first
1263   return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0;  # $x < $y => result $x
1264
1265   my $ym = $y->{_m}->copy();
1266
1267   # 2e1 => 20
1268   $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
1269
1270   # if $y has digits after dot
1271   my $shifty = 0;                       # correct _e of $x by this
1272   if ($y->{_e}->{sign} eq '-')          # has digits after dot
1273     {
1274     # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1275     $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
1276     $x->blsft($shifty,10);              # 123 => 1230, $y->{_m} is already 25
1277     }
1278   # $ym is now mantissa of $y based on exponent 0
1279
1280   my $shiftx = 0;                       # correct _e of $x by this
1281   if ($x->{_e}->{sign} eq '-')          # has digits after dot
1282     {
1283     # 123.4 % 20 => 1234 % 200
1284     $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
1285     $ym->blsft($shiftx,10);
1286     }
1287   # 123e1 % 20 => 1230 % 20
1288   if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
1289     {
1290     $x->{_m}->blsft($x->{_e},10);
1291     }
1292   $x->{_e} = $MBI->bzero() unless $x->{_e}->is_zero();
1293
1294   $x->{_e}->bsub($shiftx) if $shiftx != 0;
1295   $x->{_e}->bsub($shifty) if $shifty != 0;
1296
1297   # now mantissas are equalized, exponent of $x is adjusted, so calc result
1298
1299   $x->{_m}->bmod($ym);
1300
1301   $x->{sign} = '+' if $x->{_m}->is_zero();              # fix sign for -0
1302   $x->bnorm();
1303
1304   if ($neg != 0)        # one of them negative => correct in place
1305     {
1306     my $r = $y - $x;
1307     $x->{_m} = $r->{_m};
1308     $x->{_e} = $r->{_e};
1309     $x->{sign} = '+' if $x->{_m}->is_zero();            # fix sign for -0
1310     $x->bnorm();
1311     }
1312
1313   $x->round($a,$p,$r,$y);       # round and return
1314   }
1315
1316 sub broot
1317   {
1318   # calculate $y'th root of $x
1319   my ($self,$x,$y,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(2,@_);
1320
1321   # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1322   return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1323          $y->{sign} !~ /^\+$/;
1324
1325   return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1326
1327   # we need to limit the accuracy to protect against overflow
1328   my $fallback = 0;
1329   my (@params,$scale);
1330   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1331
1332   return $x if $x->is_nan();            # error in _find_round_parameters?
1333
1334   # no rounding at all, so must use fallback
1335   if (scalar @params == 0)
1336     {
1337     # simulate old behaviour
1338     $params[0] = $self->div_scale();    # and round to it as accuracy
1339     $scale = $params[0]+4;              # at least four more for proper round
1340     $params[2] = $r;                    # round mode by caller or undef
1341     $fallback = 1;                      # to clear a/p afterwards
1342     }
1343   else
1344     {
1345     # the 4 below is empirical, and there might be cases where it is not
1346     # enough...
1347     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1348     }
1349
1350   # when user set globals, they would interfere with our calculation, so
1351   # disable them and later re-enable them
1352   no strict 'refs';
1353   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1354   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1355   # we also need to disable any set A or P on $x (_find_round_parameters took
1356   # them already into account), since these would interfere, too
1357   delete $x->{_a}; delete $x->{_p};
1358   # need to disable $upgrade in BigInt, to avoid deep recursion
1359   local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1360
1361   # remember sign and make $x positive, since -4 ** (1/2) => -2
1362   my $sign = 0; $sign = 1 if $x->is_negative(); $x->babs();
1363
1364   if ($y->bcmp(2) == 0)         # normal square root
1365     {
1366     $x->bsqrt($scale+4);
1367     }
1368   elsif ($y->is_one('-'))
1369     {
1370     # $x ** -1 => 1/$x
1371     my $u = $self->bone()->bdiv($x,$scale);
1372     # copy private parts over
1373     $x->{_m} = $u->{_m};
1374     $x->{_e} = $u->{_e};
1375     }
1376   else
1377     {
1378     my $u = $self->bone()->bdiv($y,$scale+4);
1379     delete $u->{_a}; delete $u->{_p};           # otherwise it conflicts
1380     $x->bpow($u,$scale+4);                      # el cheapo
1381     }
1382   $x->bneg() if $sign == 1;
1383
1384   # shortcut to not run trough _find_round_parameters again
1385   if (defined $params[0])
1386     {
1387     $x->bround($params[0],$params[2]);          # then round accordingly
1388     }
1389   else
1390     {
1391     $x->bfround($params[1],$params[2]);         # then round accordingly
1392     }
1393   if ($fallback)
1394     {
1395     # clear a/p after round, since user did not request it
1396     $x->{_a} = undef; $x->{_p} = undef;
1397     }
1398   # restore globals
1399   $$abr = $ab; $$pbr = $pb;
1400   $x;
1401   }
1402
1403 sub bsqrt
1404   {
1405   # calculate square root
1406   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1407
1408   return $x->bnan() if $x->{sign} !~ /^[+]/;    # NaN, -inf or < 0
1409   return $x if $x->{sign} eq '+inf';            # sqrt(inf) == inf
1410   return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
1411
1412   # we need to limit the accuracy to protect against overflow
1413   my $fallback = 0;
1414   my (@params,$scale);
1415   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1416
1417   return $x if $x->is_nan();            # error in _find_round_parameters?
1418
1419   # no rounding at all, so must use fallback
1420   if (scalar @params == 0)
1421     {
1422     # simulate old behaviour
1423     $params[0] = $self->div_scale();    # and round to it as accuracy
1424     $scale = $params[0]+4;              # at least four more for proper round
1425     $params[2] = $r;                    # round mode by caller or undef
1426     $fallback = 1;                      # to clear a/p afterwards
1427     }
1428   else
1429     {
1430     # the 4 below is empirical, and there might be cases where it is not
1431     # enough...
1432     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1433     }
1434
1435   # when user set globals, they would interfere with our calculation, so
1436   # disable them and later re-enable them
1437   no strict 'refs';
1438   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1439   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1440   # we also need to disable any set A or P on $x (_find_round_parameters took
1441   # them already into account), since these would interfere, too
1442   delete $x->{_a}; delete $x->{_p};
1443   # need to disable $upgrade in BigInt, to avoid deep recursion
1444   local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1445
1446   my $xas = $x->as_number();
1447   my $gs = $xas->copy()->bsqrt();       # some guess
1448
1449   if (($x->{_e}->{sign} ne '-')         # guess can't be accurate if there are
1450                                         # digits after the dot
1451    && ($xas->bacmp($gs * $gs) == 0))    # guess hit the nail on the head?
1452     {
1453     # exact result
1454     $x->{_m} = $gs; $x->{_e} = $MBI->bzero(); $x->bnorm();
1455     # shortcut to not run trough _find_round_parameters again
1456     if (defined $params[0])
1457       {
1458       $x->bround($params[0],$params[2]);        # then round accordingly
1459       }
1460     else
1461       {
1462       $x->bfround($params[1],$params[2]);       # then round accordingly
1463       }
1464     if ($fallback)
1465       {
1466       # clear a/p after round, since user did not request it
1467       $x->{_a} = undef; $x->{_p} = undef;
1468       }
1469     # re-enable A and P, upgrade is taken care of by "local"
1470     ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1471     return $x;
1472     }
1473
1474   # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1475   # of the result by multipyling the input by 100 and then divide the integer
1476   # result of sqrt(input) by 10. Rounding afterwards returns the real result.
1477   # this will transform 123.456 (in $x) into 123456 (in $y1)
1478   my $y1 = $x->{_m}->copy();
1479   # We now make sure that $y1 has the same odd or even number of digits than
1480   # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1481   # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1482   # steps of 10. The length of $x does not count, since an even or odd number
1483   # of digits before the dot is not changed by adding an even number of digits
1484   # after the dot (the result is still odd or even digits long).
1485   my $length = $y1->length();
1486   $y1->bmul(10) if $x->{_e}->is_odd();
1487   # now calculate how many digits the result of sqrt(y1) would have
1488   my $digits = int($length / 2);
1489   # but we need at least $scale digits, so calculate how many are missing
1490   my $shift = $scale - $digits;
1491   # that should never happen (we take care of integer guesses above)
1492   # $shift = 0 if $shift < 0;
1493   # multiply in steps of 100, by shifting left two times the "missing" digits
1494   $y1->blsft($shift*2,10);
1495   # now take the square root and truncate to integer
1496   $y1->bsqrt();
1497   # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1498   # result, which is than later rounded to the desired scale.
1499
1500   # calculate how many zeros $x had after the '.' (or before it, depending
1501   #  on sign of $dat, the result should have half as many:
1502   my $dat = $length + $x->{_e}->numify();
1503
1504   if ($dat > 0)
1505     {
1506     # no zeros after the dot (e.g. 1.23, 0.49 etc)
1507     # preserve half as many digits before the dot than the input had
1508     # (but round this "up")
1509     $dat = int(($dat+1)/2);
1510     }
1511   else
1512     {
1513     $dat = int(($dat)/2);
1514     }
1515   $x->{_e}= $MBI->new( $dat - $y1->length() );
1516
1517   $x->{_m} = $y1;
1518
1519   # shortcut to not run trough _find_round_parameters again
1520   if (defined $params[0])
1521     {
1522     $x->bround($params[0],$params[2]);          # then round accordingly
1523     }
1524   else
1525     {
1526     $x->bfround($params[1],$params[2]);         # then round accordingly
1527     }
1528   if ($fallback)
1529     {
1530     # clear a/p after round, since user did not request it
1531     $x->{_a} = undef; $x->{_p} = undef;
1532     }
1533   # restore globals
1534   $$abr = $ab; $$pbr = $pb;
1535   $x;
1536   }
1537
1538 sub bfac
1539   {
1540   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1541   # compute factorial numbers
1542   # modifies first argument
1543   my ($self,$x,@r) = objectify(1,@_);
1544
1545   return $x->bnan()
1546     if (($x->{sign} ne '+') ||          # inf, NaN, <0 etc => NaN
1547      ($x->{_e}->{sign} ne '+'));        # digits after dot?
1548
1549   return $x->bone('+',@r) if $x->is_zero() || $x->is_one();     # 0 or 1 => 1
1550
1551   # use BigInt's bfac() for faster calc
1552   $x->{_m}->blsft($x->{_e},10);         # un-norm m
1553   $x->{_e}->bzero();                    # norm $x again
1554   $x->{_m}->bfac();                     # factorial
1555   $x->bnorm()->round(@r);
1556   }
1557
1558 sub _pow
1559   {
1560   # Calculate a power where $y is a non-integer, like 2 ** 0.5
1561   my ($x,$y,$a,$p,$r) = @_;
1562   my $self = ref($x);
1563
1564   # if $y == 0.5, it is sqrt($x)
1565   return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1566
1567   # Using:
1568   # a ** x == e ** (x * ln a)
1569
1570   # u = y * ln x
1571   #                _                         _
1572   # Taylor:       |   u    u^2    u^3         |
1573   # x ** y  = 1 + |  --- + --- + ----- + ...  |
1574   #               |_  1    1*2   1*2*3       _|
1575
1576   # we need to limit the accuracy to protect against overflow
1577   my $fallback = 0;
1578   my ($scale,@params);
1579   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1580
1581   return $x if $x->is_nan();            # error in _find_round_parameters?
1582
1583   # no rounding at all, so must use fallback
1584   if (scalar @params == 0)
1585     {
1586     # simulate old behaviour
1587     $params[0] = $self->div_scale();    # and round to it as accuracy
1588     $params[1] = undef;                 # disable P
1589     $scale = $params[0]+4;              # at least four more for proper round
1590     $params[2] = $r;                    # round mode by caller or undef
1591     $fallback = 1;                      # to clear a/p afterwards
1592     }
1593   else
1594     {
1595     # the 4 below is empirical, and there might be cases where it is not
1596     # enough...
1597     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1598     }
1599
1600   # when user set globals, they would interfere with our calculation, so
1601   # disable them and later re-enable them
1602   no strict 'refs';
1603   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1604   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1605   # we also need to disable any set A or P on $x (_find_round_parameters took
1606   # them already into account), since these would interfere, too
1607   delete $x->{_a}; delete $x->{_p};
1608   # need to disable $upgrade in BigInt, to avoid deep recursion
1609   local $Math::BigInt::upgrade = undef;
1610
1611   my ($limit,$v,$u,$below,$factor,$next,$over);
1612
1613   $u = $x->copy()->blog(undef,$scale)->bmul($y);
1614   $v = $self->bone();                           # 1
1615   $factor = $self->new(2);                      # 2
1616   $x->bone();                                   # first term: 1
1617
1618   $below = $v->copy();
1619   $over = $u->copy();
1620
1621   $limit = $self->new("1E-". ($scale-1));
1622   #my $steps = 0;
1623   while (3 < 5)
1624     {
1625     # we calculate the next term, and add it to the last
1626     # when the next term is below our limit, it won't affect the outcome
1627     # anymore, so we stop
1628     $next = $over->copy()->bdiv($below,$scale);
1629     last if $next->bacmp($limit) <= 0;
1630     $x->badd($next);
1631     # calculate things for the next term
1632     $over *= $u; $below *= $factor; $factor->binc();
1633     #$steps++;
1634     }
1635
1636   # shortcut to not run trough _find_round_parameters again
1637   if (defined $params[0])
1638     {
1639     $x->bround($params[0],$params[2]);          # then round accordingly
1640     }
1641   else
1642     {
1643     $x->bfround($params[1],$params[2]);         # then round accordingly
1644     }
1645   if ($fallback)
1646     {
1647     # clear a/p after round, since user did not request it
1648     $x->{_a} = undef; $x->{_p} = undef;
1649     }
1650   # restore globals
1651   $$abr = $ab; $$pbr = $pb;
1652   $x;
1653   }
1654
1655 sub bpow
1656   {
1657   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1658   # compute power of two numbers, second arg is used as integer
1659   # modifies first argument
1660
1661   # set up parameters
1662   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1663   # objectify is costly, so avoid it
1664   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1665     {
1666     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1667     }
1668
1669   return $x if $x->{sign} =~ /^[+-]inf$/;
1670   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1671   return $x->bone() if $y->is_zero();
1672   return $x         if $x->is_one() || $y->is_one();
1673
1674   return $x->_pow($y,$a,$p,$r) if !$y->is_int();        # non-integer power
1675
1676   my $y1 = $y->as_number();             # make bigint
1677   # if ($x == -1)
1678   if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1679     {
1680     # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
1681     return $y1->is_odd() ? $x : $x->babs(1);
1682     }
1683   if ($x->is_zero())
1684     {
1685     return $x if $y->{sign} eq '+';     # 0**y => 0 (if not y <= 0)
1686     # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1687     $x->binf();
1688     }
1689
1690   # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1691   $y1->babs();
1692   $x->{_m}->bpow($y1);
1693   $x->{_e}->bmul($y1);
1694   $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1695   $x->bnorm();
1696   if ($y->{sign} eq '-')
1697     {
1698     # modify $x in place!
1699     my $z = $x->copy(); $x->bzero()->binc();
1700     return $x->bdiv($z,$a,$p,$r);       # round in one go (might ignore y's A!)
1701     }
1702   $x->round($a,$p,$r,$y);
1703   }
1704
1705 ###############################################################################
1706 # rounding functions
1707
1708 sub bfround
1709   {
1710   # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1711   # $n == 0 means round to integer
1712   # expects and returns normalized numbers!
1713   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1714
1715   return $x if $x->modify('bfround');
1716
1717   my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1718   return $x if !defined $scale;                 # no-op
1719
1720   # never round a 0, +-inf, NaN
1721   if ($x->is_zero())
1722     {
1723     $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1724     return $x;
1725     }
1726   return $x if $x->{sign} !~ /^[+-]$/;
1727
1728   # don't round if x already has lower precision
1729   return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1730
1731   $x->{_p} = $scale;                    # remember round in any case
1732   $x->{_a} = undef;                     # and clear A
1733   if ($scale < 0)
1734     {
1735     # round right from the '.'
1736
1737     return $x if $x->{_e}->{sign} eq '+';       # e >= 0 => nothing to round
1738
1739     $scale = -$scale;                           # positive for simplicity
1740     my $len = $x->{_m}->length();               # length of mantissa
1741
1742     # the following poses a restriction on _e, but if _e is bigger than a
1743     # scalar, you got other problems (memory etc) anyway
1744     my $dad = -($x->{_e}->numify());            # digits after dot
1745     my $zad = 0;                                # zeros after dot
1746     $zad = $dad - $len if (-$dad < -$len);      # for 0.00..00xxx style
1747
1748     #print "scale $scale dad $dad zad $zad len $len\n";
1749     # number  bsstr   len zad dad
1750     # 0.123   123e-3    3   0 3
1751     # 0.0123  123e-4    3   1 4
1752     # 0.001   1e-3      1   2 3
1753     # 1.23    123e-2    3   0 2
1754     # 1.2345  12345e-4  5   0 4
1755
1756     # do not round after/right of the $dad
1757     return $x if $scale > $dad;                 # 0.123, scale >= 3 => exit
1758
1759     # round to zero if rounding inside the $zad, but not for last zero like:
1760     # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1761     return $x->bzero() if $scale < $zad;
1762     if ($scale == $zad)                 # for 0.006, scale -3 and trunc
1763       {
1764       $scale = -$len;
1765       }
1766     else
1767       {
1768       # adjust round-point to be inside mantissa
1769       if ($zad != 0)
1770         {
1771         $scale = $scale-$zad;
1772         }
1773       else
1774         {
1775         my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;    # digits before dot
1776         $scale = $dbd+$scale;
1777         }
1778       }
1779     }
1780   else
1781     {
1782     # round left from the '.'
1783
1784     # 123 => 100 means length(123) = 3 - $scale (2) => 1
1785
1786     my $dbt = $x->{_m}->length();
1787     # digits before dot
1788     my $dbd = $dbt + $x->{_e}->numify();
1789     # should be the same, so treat it as this
1790     $scale = 1 if $scale == 0;
1791     # shortcut if already integer
1792     return $x if $scale == 1 && $dbt <= $dbd;
1793     # maximum digits before dot
1794     ++$dbd;
1795
1796     if ($scale > $dbd)
1797        {
1798        # not enough digits before dot, so round to zero
1799        return $x->bzero;
1800        }
1801     elsif ( $scale == $dbd )
1802        {
1803        # maximum
1804        $scale = -$dbt;
1805        }
1806     else
1807        {
1808        $scale = $dbd - $scale;
1809        }
1810     }
1811   # pass sign to bround for rounding modes '+inf' and '-inf'
1812   $x->{_m}->{sign} = $x->{sign};
1813   $x->{_m}->bround($scale,$mode);
1814   $x->{_m}->{sign} = '+';               # fix sign back
1815   $x->bnorm();
1816   }
1817
1818 sub bround
1819   {
1820   # accuracy: preserve $N digits, and overwrite the rest with 0's
1821   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1822
1823   if (($_[0] || 0) < 0)
1824     {
1825     require Carp; Carp::croak ('bround() needs positive accuracy');
1826     }
1827
1828   my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1829   return $x if !defined $scale;                         # no-op
1830
1831   return $x if $x->modify('bround');
1832
1833   # scale is now either $x->{_a}, $accuracy, or the user parameter
1834   # test whether $x already has lower accuracy, do nothing in this case
1835   # but do round if the accuracy is the same, since a math operation might
1836   # want to round a number with A=5 to 5 digits afterwards again
1837   return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1838
1839   # scale < 0 makes no sense
1840   # never round a +-inf, NaN
1841   return $x if ($scale < 0) ||  $x->{sign} !~ /^[+-]$/;
1842
1843   # 1: $scale == 0 => keep all digits
1844   # 2: never round a 0
1845   # 3: if we should keep more digits than the mantissa has, do nothing
1846   if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1847     {
1848     $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1849     return $x;
1850     }
1851
1852   # pass sign to bround for '+inf' and '-inf' rounding modes
1853   $x->{_m}->{sign} = $x->{sign};
1854   $x->{_m}->bround($scale,$mode);       # round mantissa
1855   $x->{_m}->{sign} = '+';               # fix sign back
1856   # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1857   $x->{_a} = $scale;                    # remember rounding
1858   $x->{_p} = undef;                     # and clear P
1859   $x->bnorm();                          # del trailing zeros gen. by bround()
1860   }
1861
1862 sub bfloor
1863   {
1864   # return integer less or equal then $x
1865   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1866
1867   return $x if $x->modify('bfloor');
1868
1869   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1870
1871   # if $x has digits after dot
1872   if ($x->{_e}->{sign} eq '-')
1873     {
1874     $x->{_e}->{sign} = '+';                     # negate e
1875     $x->{_m}->brsft($x->{_e},10);               # cut off digits after dot
1876     $x->{_e}->bzero();                          # trunc/norm
1877     $x->{_m}->binc() if $x->{sign} eq '-';      # decrement if negative
1878     }
1879   $x->round($a,$p,$r);
1880   }
1881
1882 sub bceil
1883   {
1884   # return integer greater or equal then $x
1885   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1886
1887   return $x if $x->modify('bceil');
1888   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1889
1890   # if $x has digits after dot
1891   if ($x->{_e}->{sign} eq '-')
1892     {
1893     #$x->{_m}->brsft(-$x->{_e},10);
1894     #$x->{_e}->bzero();
1895     #$x++ if $x->{sign} eq '+';
1896
1897     $x->{_e}->{sign} = '+';                     # negate e
1898     $x->{_m}->brsft($x->{_e},10);               # cut off digits after dot
1899     $x->{_e}->bzero();                          # trunc/norm
1900     $x->{_m}->binc() if $x->{sign} eq '+';      # decrement if negative
1901     }
1902   $x->round($a,$p,$r);
1903   }
1904
1905 sub brsft
1906   {
1907   # shift right by $y (divide by power of $n)
1908
1909   # set up parameters
1910   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1911   # objectify is costly, so avoid it
1912   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1913     {
1914     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1915     }
1916
1917   return $x if $x->modify('brsft');
1918   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1919
1920   $n = 2 if !defined $n; $n = $self->new($n);
1921   $x->bdiv($n->bpow($y),$a,$p,$r,$y);
1922   }
1923
1924 sub blsft
1925   {
1926   # shift left by $y (multiply by power of $n)
1927
1928   # set up parameters
1929   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1930   # objectify is costly, so avoid it
1931   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1932     {
1933     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1934     }
1935
1936   return $x if $x->modify('blsft');
1937   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1938
1939   $n = 2 if !defined $n; $n = $self->new($n);
1940   $x->bmul($n->bpow($y),$a,$p,$r,$y);
1941   }
1942
1943 ###############################################################################
1944
1945 sub DESTROY
1946   {
1947   # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1948   }
1949
1950 sub AUTOLOAD
1951   {
1952   # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
1953   # or falling back to MBI::bxxx()
1954   my $name = $AUTOLOAD;
1955
1956   $name =~ s/.*:://;    # split package
1957   no strict 'refs';
1958   $class->import() if $IMPORT == 0;
1959   if (!method_alias($name))
1960     {
1961     if (!defined $name)
1962       {
1963       # delayed load of Carp and avoid recursion
1964       require Carp;
1965       Carp::croak ("Can't call a method without name");
1966       }
1967     if (!method_hand_up($name))
1968       {
1969       # delayed load of Carp and avoid recursion
1970       require Carp;
1971       Carp::croak ("Can't call $class\-\>$name, not a valid method");
1972       }
1973     # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1974     $name =~ s/^f/b/;
1975     return &{"$MBI"."::$name"}(@_);
1976     }
1977   my $bname = $name; $bname =~ s/^f/b/;
1978   *{$class."::$name"} = \&$bname;
1979   &$bname;      # uses @_
1980   }
1981
1982 sub exponent
1983   {
1984   # return a copy of the exponent
1985   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1986
1987   if ($x->{sign} !~ /^[+-]$/)
1988     {
1989     my $s = $x->{sign}; $s =~ s/^[+-]//;
1990     return $self->new($s);                      # -inf, +inf => +inf
1991     }
1992   return $x->{_e}->copy();
1993   }
1994
1995 sub mantissa
1996   {
1997   # return a copy of the mantissa
1998   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1999
2000   if ($x->{sign} !~ /^[+-]$/)
2001     {
2002     my $s = $x->{sign}; $s =~ s/^[+]//;
2003     return $self->new($s);                      # -inf, +inf => +inf
2004     }
2005   my $m = $x->{_m}->copy();             # faster than going via bstr()
2006   $m->bneg() if $x->{sign} eq '-';
2007
2008   $m;
2009   }
2010
2011 sub parts
2012   {
2013   # return a copy of both the exponent and the mantissa
2014   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2015
2016   if ($x->{sign} !~ /^[+-]$/)
2017     {
2018     my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2019     return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2020     }
2021   my $m = $x->{_m}->copy();     # faster than going via bstr()
2022   $m->bneg() if $x->{sign} eq '-';
2023   return ($m,$x->{_e}->copy());
2024   }
2025
2026 ##############################################################################
2027 # private stuff (internal use only)
2028
2029 sub import
2030   {
2031   my $self = shift;
2032   my $l = scalar @_;
2033   my $lib = ''; my @a;
2034   $IMPORT=1;
2035   for ( my $i = 0; $i < $l ; $i++)
2036     {
2037     if ( $_[$i] eq ':constant' )
2038       {
2039       # this rest causes overlord er load to step in
2040       overload::constant float => sub { $self->new(shift); };
2041       }
2042     elsif ($_[$i] eq 'upgrade')
2043       {
2044       # this causes upgrading
2045       $upgrade = $_[$i+1];              # or undef to disable
2046       $i++;
2047       }
2048     elsif ($_[$i] eq 'downgrade')
2049       {
2050       # this causes downgrading
2051       $downgrade = $_[$i+1];            # or undef to disable
2052       $i++;
2053       }
2054     elsif ($_[$i] eq 'lib')
2055       {
2056       # alternative library
2057       $lib = $_[$i+1] || '';            # default Calc
2058       $i++;
2059       }
2060     elsif ($_[$i] eq 'with')
2061       {
2062       # alternative class for our private parts()
2063       $MBI = $_[$i+1] || 'Math::BigInt';        # default Math::BigInt
2064       $i++;
2065       }
2066     else
2067       {
2068       push @a, $_[$i];
2069       }
2070     }
2071
2072   # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2073   my $mbilib = eval { Math::BigInt->config()->{lib} };
2074   if ((defined $mbilib) && ($MBI eq 'Math::BigInt'))
2075     {
2076     # MBI already loaded
2077     $MBI->import('lib',"$lib,$mbilib", 'objectify');
2078     }
2079   else
2080     {
2081     # MBI not loaded, or with ne "Math::BigInt"
2082     $lib .= ",$mbilib" if defined $mbilib;
2083     $lib =~ s/^,//;                             # don't leave empty
2084     # replacement library can handle lib statement, but also could ignore it
2085     if ($] < 5.006)
2086       {
2087       # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2088       # used in the same script, or eval inside import().
2089       my @parts = split /::/, $MBI;             # Math::BigInt => Math BigInt
2090       my $file = pop @parts; $file .= '.pm';    # BigInt => BigInt.pm
2091       require File::Spec;
2092       $file = File::Spec->catfile (@parts, $file);
2093       eval { require "$file"; };
2094       $MBI->import( lib => $lib, 'objectify' );
2095       }
2096     else
2097       {
2098       my $rc = "use $MBI lib => '$lib', 'objectify';";
2099       eval $rc;
2100       }
2101     }
2102   if ($@)
2103     {
2104     require Carp; Carp::croak ("Couldn't load $MBI: $! $@");
2105     }
2106
2107   # any non :constant stuff is handled by our parent, Exporter
2108   # even if @_ is empty, to give it a chance
2109   $self->SUPER::import(@a);             # for subclasses
2110   $self->export_to_level(1,$self,@a);   # need this, too
2111   }
2112
2113 sub bnorm
2114   {
2115   # adjust m and e so that m is smallest possible
2116   # round number according to accuracy and precision settings
2117   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2118
2119   return $x if $x->{sign} !~ /^[+-]$/;          # inf, nan etc
2120
2121 #  if (!$x->{_m}->is_odd())
2122 #    {
2123     my $zeros = $x->{_m}->_trailing_zeros();    # correct for trailing zeros
2124     if ($zeros != 0)
2125       {
2126       $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros);
2127       }
2128     # for something like 0Ey, set y to 1, and -0 => +0
2129     $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
2130 #    }
2131   # this is to prevent automatically rounding when MBI's globals are set
2132   $x->{_m}->{_f} = MB_NEVER_ROUND;
2133   $x->{_e}->{_f} = MB_NEVER_ROUND;
2134   # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
2135   $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
2136   $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
2137   $x;                                   # MBI bnorm is no-op, so dont call it
2138   }
2139
2140 ##############################################################################
2141
2142 sub as_hex
2143   {
2144   # return number as hexadecimal string (only for integers defined)
2145   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2146
2147   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
2148   return '0x0' if $x->is_zero();
2149
2150   return $nan if $x->{_e}->{sign} ne '+';       # how to do 1e-1 in hex!?
2151
2152   my $z = $x->{_m}->copy();
2153   if (!$x->{_e}->is_zero())             # > 0
2154     {
2155     $z->blsft($x->{_e},10);
2156     }
2157   $z->{sign} = $x->{sign};
2158   $z->as_hex();
2159   }
2160
2161 sub as_bin
2162   {
2163   # return number as binary digit string (only for integers defined)
2164   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2165
2166   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
2167   return '0b0' if $x->is_zero();
2168
2169   return $nan if $x->{_e}->{sign} ne '+';       # how to do 1e-1 in hex!?
2170
2171   my $z = $x->{_m}->copy();
2172   if (!$x->{_e}->is_zero())             # > 0
2173     {
2174     $z->blsft($x->{_e},10);
2175     }
2176   $z->{sign} = $x->{sign};
2177   $z->as_bin();
2178   }
2179
2180 sub as_number
2181   {
2182   # return copy as a bigint representation of this BigFloat number
2183   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2184
2185   my $z = $x->{_m}->copy();
2186   if ($x->{_e}->{sign} eq '-')          # < 0
2187     {
2188     $x->{_e}->{sign} = '+';             # flip
2189     $z->brsft($x->{_e},10);
2190     $x->{_e}->{sign} = '-';             # flip back
2191     }
2192   elsif (!$x->{_e}->is_zero())          # > 0
2193     {
2194     $z->blsft($x->{_e},10);
2195     }
2196   $z->{sign} = $x->{sign};
2197   $z;
2198   }
2199
2200 sub length
2201   {
2202   my $x = shift;
2203   my $class = ref($x) || $x;
2204   $x = $class->new(shift) unless ref($x);
2205
2206   return 1 if $x->{_m}->is_zero();
2207   my $len = $x->{_m}->length();
2208   $len += $x->{_e} if $x->{_e}->sign() eq '+';
2209   if (wantarray())
2210     {
2211     my $t = $MBI->bzero();
2212     $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
2213     return ($len,$t);
2214     }
2215   $len;
2216   }
2217
2218 1;
2219 __END__
2220
2221 =head1 NAME
2222
2223 Math::BigFloat - Arbitrary size floating point math package
2224
2225 =head1 SYNOPSIS
2226
2227   use Math::BigFloat;
2228
2229   # Number creation
2230   $x = Math::BigFloat->new($str);       # defaults to 0
2231   $nan  = Math::BigFloat->bnan();       # create a NotANumber
2232   $zero = Math::BigFloat->bzero();      # create a +0
2233   $inf = Math::BigFloat->binf();        # create a +inf
2234   $inf = Math::BigFloat->binf('-');     # create a -inf
2235   $one = Math::BigFloat->bone();        # create a +1
2236   $one = Math::BigFloat->bone('-');     # create a -1
2237
2238   # Testing
2239   $x->is_zero();                # true if arg is +0
2240   $x->is_nan();                 # true if arg is NaN
2241   $x->is_one();                 # true if arg is +1
2242   $x->is_one('-');              # true if arg is -1
2243   $x->is_odd();                 # true if odd, false for even
2244   $x->is_even();                # true if even, false for odd
2245   $x->is_positive();            # true if >= 0
2246   $x->is_negative();            # true if <  0
2247   $x->is_inf(sign);             # true if +inf, or -inf (default is '+')
2248
2249   $x->bcmp($y);                 # compare numbers (undef,<0,=0,>0)
2250   $x->bacmp($y);                # compare absolutely (undef,<0,=0,>0)
2251   $x->sign();                   # return the sign, either +,- or NaN
2252   $x->digit($n);                # return the nth digit, counting from right
2253   $x->digit(-$n);               # return the nth digit, counting from left
2254
2255   # The following all modify their first argument. If you want to preserve
2256   # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2257   # neccessary when mixing $a = $b assigments with non-overloaded math.
2258
2259   # set
2260   $x->bzero();                  # set $i to 0
2261   $x->bnan();                   # set $i to NaN
2262   $x->bone();                   # set $x to +1
2263   $x->bone('-');                # set $x to -1
2264   $x->binf();                   # set $x to inf
2265   $x->binf('-');                # set $x to -inf
2266
2267   $x->bneg();                   # negation
2268   $x->babs();                   # absolute value
2269   $x->bnorm();                  # normalize (no-op)
2270   $x->bnot();                   # two's complement (bit wise not)
2271   $x->binc();                   # increment x by 1
2272   $x->bdec();                   # decrement x by 1
2273
2274   $x->badd($y);                 # addition (add $y to $x)
2275   $x->bsub($y);                 # subtraction (subtract $y from $x)
2276   $x->bmul($y);                 # multiplication (multiply $x by $y)
2277   $x->bdiv($y);                 # divide, set $x to quotient
2278                                 # return (quo,rem) or quo if scalar
2279
2280   $x->bmod($y);                 # modulus ($x % $y)
2281   $x->bpow($y);                 # power of arguments ($x ** $y)
2282   $x->blsft($y);                # left shift
2283   $x->brsft($y);                # right shift
2284                                 # return (quo,rem) or quo if scalar
2285
2286   $x->blog();                   # logarithm of $x to base e (Euler's number)
2287   $x->blog($base);              # logarithm of $x to base $base (f.i. 2)
2288
2289   $x->band($y);                 # bit-wise and
2290   $x->bior($y);                 # bit-wise inclusive or
2291   $x->bxor($y);                 # bit-wise exclusive or
2292   $x->bnot();                   # bit-wise not (two's complement)
2293
2294   $x->bsqrt();                  # calculate square-root
2295   $x->broot($y);                # $y'th root of $x (e.g. $y == 3 => cubic root)
2296   $x->bfac();                   # factorial of $x (1*2*3*4*..$x)
2297
2298   $x->bround($N);               # accuracy: preserve $N digits
2299   $x->bfround($N);              # precision: round to the $Nth digit
2300
2301   $x->bfloor();                 # return integer less or equal than $x
2302   $x->bceil();                  # return integer greater or equal than $x
2303
2304   # The following do not modify their arguments:
2305
2306   bgcd(@values);                # greatest common divisor
2307   blcm(@values);                # lowest common multiplicator
2308
2309   $x->bstr();                   # return string
2310   $x->bsstr();                  # return string in scientific notation
2311
2312   $x->exponent();               # return exponent as BigInt
2313   $x->mantissa();               # return mantissa as BigInt
2314   $x->parts();                  # return (mantissa,exponent) as BigInt
2315
2316   $x->length();                 # number of digits (w/o sign and '.')
2317   ($l,$f) = $x->length();       # number of digits, and length of fraction
2318
2319   $x->precision();              # return P of $x (or global, if P of $x undef)
2320   $x->precision($n);            # set P of $x to $n
2321   $x->accuracy();               # return A of $x (or global, if A of $x undef)
2322   $x->accuracy($n);             # set A $x to $n
2323
2324   # these get/set the appropriate global value for all BigFloat objects
2325   Math::BigFloat->precision();  # Precision
2326   Math::BigFloat->accuracy();   # Accuracy
2327   Math::BigFloat->round_mode(); # rounding mode
2328
2329 =head1 DESCRIPTION
2330
2331 All operators (inlcuding basic math operations) are overloaded if you
2332 declare your big floating point numbers as
2333
2334   $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2335
2336 Operations with overloaded operators preserve the arguments, which is
2337 exactly what you expect.
2338
2339 =head2 Canonical notation
2340
2341 Input to these routines are either BigFloat objects, or strings of the
2342 following four forms:
2343
2344 =over 2
2345
2346 =item *
2347
2348 C</^[+-]\d+$/>
2349
2350 =item *
2351
2352 C</^[+-]\d+\.\d*$/>
2353
2354 =item *
2355
2356 C</^[+-]\d+E[+-]?\d+$/>
2357
2358 =item *
2359
2360 C</^[+-]\d*\.\d+E[+-]?\d+$/>
2361
2362 =back
2363
2364 all with optional leading and trailing zeros and/or spaces. Additonally,
2365 numbers are allowed to have an underscore between any two digits.
2366
2367 Empty strings as well as other illegal numbers results in 'NaN'.
2368
2369 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2370 are always stored in normalized form. On a string, it creates a BigFloat
2371 object.
2372
2373 =head2 Output
2374
2375 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2376
2377 The string output will always have leading and trailing zeros stripped and drop
2378 a plus sign. C<bstr()> will give you always the form with a decimal point,
2379 while C<bsstr()> (s for scientific) gives you the scientific notation.
2380
2381         Input                   bstr()          bsstr()
2382         '-0'                    '0'             '0E1'
2383         '  -123 123 123'        '-123123123'    '-123123123E0'
2384         '00.0123'               '0.0123'        '123E-4'
2385         '123.45E-2'             '1.2345'        '12345E-4'
2386         '10E+3'                 '10000'         '1E4'
2387
2388 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2389 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2390 return either undef, <0, 0 or >0 and are suited for sort.
2391
2392 Actual math is done by using the class defined with C<with => Class;> (which
2393 defaults to BigInts) to represent the mantissa and exponent.
2394
2395 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2396 represent the result when input arguments are not numbers, as well as
2397 the result of dividing by zero.
2398
2399 =head2 C<mantissa()>, C<exponent()> and C<parts()>
2400
2401 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2402 as BigInts such that:
2403
2404         $m = $x->mantissa();
2405         $e = $x->exponent();
2406         $y = $m * ( 10 ** $e );
2407         print "ok\n" if $x == $y;
2408
2409 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2410
2411 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2412
2413 Currently the mantissa is reduced as much as possible, favouring higher
2414 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2415 This might change in the future, so do not depend on it.
2416
2417 =head2 Accuracy vs. Precision
2418
2419 See also: L<Rounding|Rounding>.
2420
2421 Math::BigFloat supports both precision and accuracy. For a full documentation,
2422 examples and tips on these topics please see the large section in
2423 L<Math::BigInt>.
2424
2425 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2426 a operation consumes all resources, each operation produces no more than
2427 the requested number of digits.
2428
2429 Please refer to BigInt's documentation for the precedence rules of which
2430 accuracy/precision setting will be used.
2431
2432 If there is no gloabl precision set, B<and> the operation inquestion was not
2433 called with a requested precision or accuracy, B<and> the input $x has no
2434 accuracy or precision set, then a fallback parameter will be used. For
2435 historical reasons, it is called C<div_scale> and can be accessed via:
2436
2437         $d = Math::BigFloat->div_scale();               # query
2438         Math::BigFloat->div_scale($n);                  # set to $n digits
2439
2440 The default value is 40 digits.
2441
2442 In case the result of one operation has more precision than specified,
2443 it is rounded. The rounding mode taken is either the default mode, or the one
2444 supplied to the operation after the I<scale>:
2445
2446         $x = Math::BigFloat->new(2);
2447         Math::BigFloat->precision(5);           # 5 digits max
2448         $y = $x->copy()->bdiv(3);               # will give 0.66666
2449         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
2450         $y = $x->copy()->bdiv(3,6,'odd');       # will give 0.666667
2451         Math::BigFloat->round_mode('zero');
2452         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
2453
2454 =head2 Rounding
2455
2456 =over 2
2457
2458 =item ffround ( +$scale )
2459
2460 Rounds to the $scale'th place left from the '.', counting from the dot.
2461 The first digit is numbered 1.
2462
2463 =item ffround ( -$scale )
2464
2465 Rounds to the $scale'th place right from the '.', counting from the dot.
2466
2467 =item ffround ( 0 )
2468
2469 Rounds to an integer.
2470
2471 =item fround  ( +$scale )
2472
2473 Preserves accuracy to $scale digits from the left (aka significant digits)
2474 and pads the rest with zeros. If the number is between 1 and -1, the
2475 significant digits count from the first non-zero after the '.'
2476
2477 =item fround  ( -$scale ) and fround ( 0 )
2478
2479 These are effectively no-ops.
2480
2481 =back
2482
2483 All rounding functions take as a second parameter a rounding mode from one of
2484 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2485
2486 The default rounding mode is 'even'. By using
2487 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
2488 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2489 no longer supported.
2490 The second parameter to the round functions then overrides the default
2491 temporarily.
2492
2493 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
2494 'trunc' as rounding mode to make it equivalent to:
2495
2496         $x = 2.5;
2497         $y = int($x) + 2;
2498
2499 You can override this by passing the desired rounding mode as parameter to
2500 C<as_number()>:
2501
2502         $x = Math::BigFloat->new(2.5);
2503         $y = $x->as_number('odd');      # $y = 3
2504
2505 =head1 EXAMPLES
2506
2507   # not ready yet
2508
2509 =head1 Autocreating constants
2510
2511 After C<use Math::BigFloat ':constant'> all the floating point constants
2512 in the given scope are converted to C<Math::BigFloat>. This conversion
2513 happens at compile time.
2514
2515 In particular
2516
2517   perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2518
2519 prints the value of C<2E-100>. Note that without conversion of
2520 constants the expression 2E-100 will be calculated as normal floating point
2521 number.
2522
2523 Please note that ':constant' does not affect integer constants, nor binary
2524 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2525 work.
2526
2527 =head2 Math library
2528
2529 Math with the numbers is done (by default) by a module called
2530 Math::BigInt::Calc. This is equivalent to saying:
2531
2532         use Math::BigFloat lib => 'Calc';
2533
2534 You can change this by using:
2535
2536         use Math::BigFloat lib => 'BitVect';
2537
2538 The following would first try to find Math::BigInt::Foo, then
2539 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2540
2541         use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2542
2543 Calc.pm uses as internal format an array of elements of some decimal base
2544 (usually 1e7, but this might be differen for some systems) with the least
2545 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2546 significant bit first. Other modules might use even different means of
2547 representing the numbers. See the respective module documentation for further
2548 details.
2549
2550 Please note that Math::BigFloat does B<not> use the denoted library itself,
2551 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2552 to do:
2553
2554         use Math::BigInt lib => 'GMP';
2555         use Math::BigFloat;
2556
2557 you can roll it all into one line:
2558
2559         use Math::BigFloat lib => 'GMP';
2560
2561 It is also possible to just require Math::BigFloat:
2562
2563         require Math::BigFloat;
2564
2565 This will load the neccessary things (like BigInt) when they are needed, and
2566 automatically.
2567
2568 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2569 you ever wanted to know about loading a different library.
2570
2571 =head2 Using Math::BigInt::Lite
2572
2573 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2574
2575         # 1
2576         use Math::BigFloat with => 'Math::BigInt::Lite';
2577
2578 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2579 can combine these if you want. For instance, you may want to use
2580 Math::BigInt objects in your main script, too.
2581
2582         # 2
2583         use Math::BigInt;
2584         use Math::BigFloat with => 'Math::BigInt::Lite';
2585
2586 Of course, you can combine this with the C<lib> parameter.
2587
2588         # 3
2589         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2590
2591 There is no need for a "use Math::BigInt;" statement, even if you want to
2592 use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2593 always loads it. But if you add it, add it B<before>:
2594
2595         # 4
2596         use Math::BigInt;
2597         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2598
2599 Notice that the module with the last C<lib> will "win" and thus
2600 it's lib will be used if the lib is available:
2601
2602         # 5
2603         use Math::BigInt lib => 'Bar,Baz';
2604         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2605
2606 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2607 words, Math::BigFloat will try to retain previously loaded libs when you
2608 don't specify it onem but if you specify one, it will try to load them.
2609
2610 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2611 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2612 same as trying the latter load alone, except for the fact that one of Bar or
2613 Baz might be loaded needlessly in an intermidiate step (and thus hang around
2614 and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2615 will still be tried to be loaded, but this is not as time/memory consuming as
2616 actually loading one of them. Still, this type of usage is not recommended due
2617 to these issues.
2618
2619 The old way (loading the lib only in BigInt) still works though:
2620
2621         # 6
2622         use Math::BigInt lib => 'Bar,Baz';
2623         use Math::BigFloat;
2624
2625 You can even load Math::BigInt afterwards:
2626
2627         # 7
2628         use Math::BigFloat;
2629         use Math::BigInt lib => 'Bar,Baz';
2630
2631 But this has the same problems like #5, it will first load Calc
2632 (Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2633 Baz, depending on which of them works and is usable/loadable. Since this
2634 loads Calc unnecc., it is not recommended.
2635
2636 Since it also possible to just require Math::BigFloat, this poses the question
2637 about what libary this will use:
2638
2639         require Math::BigFloat;
2640         my $x = Math::BigFloat->new(123); $x += 123;
2641
2642 It will use Calc. Please note that the call to import() is still done, but
2643 only when you use for the first time some Math::BigFloat math (it is triggered
2644 via any constructor, so the first time you create a Math::BigFloat, the load
2645 will happen in the background). This means:
2646
2647         require Math::BigFloat;
2648         Math::BigFloat->import ( lib => 'Foo,Bar' );
2649
2650 would be the same as:
2651
2652         use Math::BigFloat lib => 'Foo, Bar';
2653
2654 But don't try to be clever to insert some operations in between:
2655
2656         require Math::BigFloat;
2657         my $x = Math::BigFloat->bone() + 4;             # load BigInt and Calc
2658         Math::BigFloat->import( lib => 'Pari' );        # load Pari, too
2659         $x = Math::BigFloat->bone()+4;                  # now use Pari
2660
2661 While this works, it loads Calc needlessly. But maybe you just wanted that?
2662
2663 B<Examples #3 is highly recommended> for daily usage.
2664
2665 =head1 BUGS
2666
2667 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
2668
2669 =head1 CAVEATS
2670
2671 =over 1
2672
2673 =item stringify, bstr()
2674
2675 Both stringify and bstr() now drop the leading '+'. The old code would return
2676 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2677 reasoning and details.
2678
2679 =item bdiv
2680
2681 The following will probably not do what you expect:
2682
2683         print $c->bdiv(123.456),"\n";
2684
2685 It prints both quotient and reminder since print works in list context. Also,
2686 bdiv() will modify $c, so be carefull. You probably want to use
2687
2688         print $c / 123.456,"\n";
2689         print scalar $c->bdiv(123.456),"\n";  # or if you want to modify $c
2690
2691 instead.
2692
2693 =item Modifying and =
2694
2695 Beware of:
2696
2697         $x = Math::BigFloat->new(5);
2698         $y = $x;
2699
2700 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2701 a second reference to the B<same> object and stores it in $y. Thus anything
2702 that modifies $x will modify $y (except overloaded math operators), and vice
2703 versa. See L<Math::BigInt> for details and how to avoid that.
2704
2705 =item bpow
2706
2707 C<bpow()> now modifies the first argument, unlike the old code which left
2708 it alone and only returned the result. This is to be consistent with
2709 C<badd()> etc. The first will modify $x, the second one won't:
2710
2711         print bpow($x,$i),"\n";         # modify $x
2712         print $x->bpow($i),"\n";        # ditto
2713         print $x ** $i,"\n";            # leave $x alone
2714
2715 =back
2716
2717 =head1 SEE ALSO
2718
2719 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
2720 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and  L<Math::BigInt::GMP>.
2721
2722 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
2723 because they solve the autoupgrading/downgrading issue, at least partly.
2724
2725 The package at
2726 L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
2727 more documentation including a full version history, testcases, empty
2728 subclass files and benchmarks.
2729
2730 =head1 LICENSE
2731
2732 This program is free software; you may redistribute it and/or modify it under
2733 the same terms as Perl itself.
2734
2735 =head1 AUTHORS
2736
2737 Mark Biggar, overloaded interface by Ilya Zakharevich.
2738 Completely rewritten by Tels http://bloodgate.com in 2001, 2002, and still
2739 at it in 2003.
2740
2741 =cut