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