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