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