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