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