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