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