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