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