lib/Math/BigFloat.pm

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