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