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