1 package Math::BigFloat;
4 # Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
7 # The following hash values are internally used:
8 # _e: exponent (BigInt)
9 # _m: mantissa (absolute BigInt)
10 # sign: +,-,+inf,-inf, or "NaN" if not a number
13 # _f: flags, used to signal MBI not to touch our private parts
19 @ISA = qw(Exporter Math::BigInt);
22 # $_trap_inf and $_trap_nan are internal and should never be accessed from the outside
23 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode
24 $upgrade $downgrade $_trap_nan $_trap_inf/;
25 my $class = "Math::BigFloat";
28 '<=>' => sub { $_[2] ?
29 ref($_[0])->bcmp($_[1],$_[0]) :
30 ref($_[0])->bcmp($_[0],$_[1])},
31 'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
34 ##############################################################################
35 # global constants, flags and assorted stuff
37 # the following are public, but their usage is not recommended. Use the
38 # accessor methods instead.
40 # class constants, use Class->constant_name() to access
41 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
48 my $MBI = 'Math::BigInt'; # the package we are using for our private parts
49 # changable by use Math::BigFloat with => 'package'
51 # the following are private and not to be used from the outside:
53 sub MB_NEVER_ROUND () { 0x0001; }
55 # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
60 # constant for easier life
63 my $IMPORT = 0; # was import() called yet?
64 # used to make require work
66 # some digits of accuracy for blog(undef,10); which we use in blog() for speed
68 '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
69 my $LOG_10_A = length($LOG_10)-1;
72 '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
73 my $LOG_2_A = length($LOG_2)-1;
75 ##############################################################################
76 # the old code had $rnd_mode, so we need to support it, too
78 sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
79 sub FETCH { return $round_mode; }
80 sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
84 # when someone set's $rnd_mode, we catch this and check the value to see
85 # whether it is valid or not.
86 $rnd_mode = 'even'; tie $rnd_mode, 'Math::BigFloat';
89 ##############################################################################
91 # in case we call SUPER::->foo() and this wants to call modify()
92 # sub modify () { 0; }
95 # valid method aliases for AUTOLOAD
96 my %methods = map { $_ => 1 }
97 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
98 fint facmp fcmp fzero fnan finf finc fdec flog ffac
99 fceil ffloor frsft flsft fone flog froot
101 # valid method's that can be hand-ed up (for AUTOLOAD)
102 my %hand_ups = map { $_ => 1 }
103 qw / is_nan is_inf is_negative is_positive
104 accuracy precision div_scale round_mode fneg fabs fnot
105 objectify upgrade downgrade
109 sub method_alias { return exists $methods{$_[0]||''}; }
110 sub method_hand_up { return exists $hand_ups{$_[0]||''}; }
113 ##############################################################################
118 # create a new BigFloat object from a string or another bigfloat object.
121 # sign => sign (+/-), or "NaN"
123 my ($class,$wanted,@r) = @_;
125 # avoid numify-calls by not using || on $wanted!
126 return $class->bzero() if !defined $wanted; # default to 0
127 return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
129 $class->import() if $IMPORT == 0; # make require work
131 my $self = {}; bless $self, $class;
132 # shortcut for bigints and its subclasses
133 if ((ref($wanted)) && (ref($wanted) ne $class))
135 $self->{_m} = $wanted->as_number(); # get us a bigint copy
136 $self->{_e} = $MBI->bzero();
138 $self->{sign} = $wanted->sign();
139 return $self->bnorm();
142 # handle '+inf', '-inf' first
143 if ($wanted =~ /^[+-]?inf$/)
145 return $downgrade->new($wanted) if $downgrade;
147 $self->{_e} = $MBI->bzero();
148 $self->{_m} = $MBI->bzero();
149 $self->{sign} = $wanted;
150 $self->{sign} = '+inf' if $self->{sign} eq 'inf';
151 return $self->bnorm();
153 #print "new string '$wanted'\n";
155 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
161 Carp::croak ("$wanted is not a number initialized to $class");
164 return $downgrade->bnan() if $downgrade;
166 $self->{_e} = $MBI->bzero();
167 $self->{_m} = $MBI->bzero();
168 $self->{sign} = $nan;
172 # make integer from mantissa by adjusting exp, then convert to bigint
173 # undef,undef to signal MBI that we don't need no bloody rounding
174 $self->{_e} = $MBI->new("$$es$$ev",undef,undef); # exponent
175 $self->{_m} = $MBI->new("$$miv$$mfv",undef,undef); # create mant.
177 # this is to prevent automatically rounding when MBI's globals are set
178 $self->{_m}->{_f} = MB_NEVER_ROUND;
179 $self->{_e}->{_f} = MB_NEVER_ROUND;
181 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
182 $self->{_e}->bsub( $MBI->new(CORE::length($$mfv),undef,undef))
183 if CORE::length($$mfv) != 0;
184 $self->{sign} = $$mis;
186 #print "$$miv$$mfv $$es$$ev\n";
188 # we can only have trailing zeros on the mantissa of $$mfv eq ''
189 if (CORE::length($$mfv) == 0)
191 my $zeros = $self->{_m}->_trailing_zeros(); # correct for trailing zeros
194 $self->{_m}->brsft($zeros,10); $self->{_e}->badd($MBI->new($zeros));
199 # for something like 0Ey, set y to 1, and -0 => +0
200 $self->{sign} = '+', $self->{_e}->bone() if $self->{_m}->is_zero();
202 return $self->round(@r) if !$downgrade;
204 # if downgrade, inf, NaN or integers go down
206 if ($downgrade && $self->{_e}->{sign} eq '+')
208 #print "downgrading $$miv$$mfv"."E$$es$$ev";
209 if ($self->{_e}->is_zero())
211 $self->{_m}->{sign} = $$mis; # negative if wanted
212 return $downgrade->new($self->{_m});
214 return $downgrade->new($self->bsstr());
216 #print "mbf new $self->{sign} $self->{_m} e $self->{_e} ",ref($self),"\n";
217 $self->bnorm()->round(@r); # first normalize, then round
222 # used by parent class bone() to initialize number to NaN
228 my $class = ref($self);
229 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
232 $IMPORT=1; # call our import only once
233 $self->{_m} = $MBI->bzero();
234 $self->{_e} = $MBI->bzero();
239 # used by parent class bone() to initialize number to +-inf
245 my $class = ref($self);
246 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
249 $IMPORT=1; # call our import only once
250 $self->{_m} = $MBI->bzero();
251 $self->{_e} = $MBI->bzero();
256 # used by parent class bone() to initialize number to 1
258 $IMPORT=1; # call our import only once
259 $self->{_m} = $MBI->bone();
260 $self->{_e} = $MBI->bzero();
265 # used by parent class bone() to initialize number to 0
267 $IMPORT=1; # call our import only once
268 $self->{_m} = $MBI->bzero();
269 $self->{_e} = $MBI->bone();
274 my ($self,$class) = @_;
275 return if $class =~ /^Math::BigInt/; # we aren't one of these
276 UNIVERSAL::isa($self,$class);
281 # return (later set?) configuration data as hash ref
282 my $class = shift || 'Math::BigFloat';
284 my $cfg = $class->SUPER::config(@_);
286 # now we need only to override the ones that are different from our parent
287 $cfg->{class} = $class;
292 ##############################################################################
293 # string conversation
297 # (ref to BFLOAT or num_str ) return num_str
298 # Convert number from internal format to (non-scientific) string format.
299 # internal format is always normalized (no leading zeros, "-0" => "+0")
300 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
302 if ($x->{sign} !~ /^[+-]$/)
304 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
308 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
311 my $not_zero = !($x->{sign} eq '+' && $x->{_m}->is_zero());
314 $es = $x->{_m}->bstr();
315 $len = CORE::length($es);
316 my $e = $x->{_e}->numify();
320 # if _e is bigger than a scalar, the following will blow your memory
323 #print "style: 0.xxxx\n";
324 my $r = abs($e) - $len;
325 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
329 #print "insert '.' at $e in '$es'\n";
330 substr($es,$e,0) = '.'; $cad = $x->{_e};
336 $es .= '0' x $e; $len += $e; $cad = 0;
339 $es = '-'.$es if $x->{sign} eq '-';
340 # if set accuracy or precision, pad with zeros on the right side
341 if ((defined $x->{_a}) && ($not_zero))
343 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
344 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
345 $zeros = $x->{_a} - $len if $cad != $len;
346 $es .= $dot.'0' x $zeros if $zeros > 0;
348 elsif ((($x->{_p} || 0) < 0))
350 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
351 my $zeros = -$x->{_p} + $cad;
352 $es .= $dot.'0' x $zeros if $zeros > 0;
359 # (ref to BFLOAT or num_str ) return num_str
360 # Convert number from internal format to scientific string format.
361 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
362 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
364 if ($x->{sign} !~ /^[+-]$/)
366 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
369 # do $esign, because we need '1e+1', since $x->{_e}->bstr() misses the +
370 my $esign = $x->{_e}->{sign}; $esign = '' if $esign eq '-';
371 my $sep = 'e'.$esign;
372 my $sign = $x->{sign}; $sign = '' if $sign eq '+';
373 $sign . $x->{_m}->bstr() . $sep . $x->{_e}->bstr();
378 # Make a number from a BigFloat object
379 # simple return a string and let Perl's atoi()/atof() handle the rest
380 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
384 ##############################################################################
385 # public stuff (usually prefixed with "b")
388 # XXX TODO this must be overwritten and return NaN for non-integer values
389 # band(), bior(), bxor(), too
392 # $class->SUPER::bnot($class,@_);
397 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
400 my ($self,$x,$y) = (ref($_[0]),@_);
401 # objectify is costly, so avoid it
402 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
404 ($self,$x,$y) = objectify(2,@_);
407 return $upgrade->bcmp($x,$y) if defined $upgrade &&
408 ((!$x->isa($self)) || (!$y->isa($self)));
410 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
412 # handle +-inf and NaN
413 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
414 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
415 return +1 if $x->{sign} eq '+inf';
416 return -1 if $x->{sign} eq '-inf';
417 return -1 if $y->{sign} eq '+inf';
421 # check sign for speed first
422 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
423 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
426 my $xz = $x->is_zero();
427 my $yz = $y->is_zero();
428 return 0 if $xz && $yz; # 0 <=> 0
429 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
430 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
432 # adjust so that exponents are equal
433 my $lxm = $x->{_m}->length();
434 my $lym = $y->{_m}->length();
435 # the numify somewhat limits our length, but makes it much faster
436 my $lx = $lxm + $x->{_e}->numify();
437 my $ly = $lym + $y->{_e}->numify();
438 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
439 return $l <=> 0 if $l != 0;
441 # lengths (corrected by exponent) are equal
442 # so make mantissa equal length by padding with zero (shift left)
443 my $diff = $lxm - $lym;
444 my $xm = $x->{_m}; # not yet copy it
448 $ym = $y->{_m}->copy()->blsft($diff,10);
452 $xm = $x->{_m}->copy()->blsft(-$diff,10);
454 my $rc = $xm->bacmp($ym);
455 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
461 # Compares 2 values, ignoring their signs.
462 # Returns one of undef, <0, =0, >0. (suitable for sort)
465 my ($self,$x,$y) = (ref($_[0]),@_);
466 # objectify is costly, so avoid it
467 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
469 ($self,$x,$y) = objectify(2,@_);
472 return $upgrade->bacmp($x,$y) if defined $upgrade &&
473 ((!$x->isa($self)) || (!$y->isa($self)));
475 # handle +-inf and NaN's
476 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
478 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
479 return 0 if ($x->is_inf() && $y->is_inf());
480 return 1 if ($x->is_inf() && !$y->is_inf());
485 my $xz = $x->is_zero();
486 my $yz = $y->is_zero();
487 return 0 if $xz && $yz; # 0 <=> 0
488 return -1 if $xz && !$yz; # 0 <=> +y
489 return 1 if $yz && !$xz; # +x <=> 0
491 # adjust so that exponents are equal
492 my $lxm = $x->{_m}->length();
493 my $lym = $y->{_m}->length();
494 # the numify somewhat limits our length, but makes it much faster
495 my $lx = $lxm + $x->{_e}->numify();
496 my $ly = $lym + $y->{_e}->numify();
498 return $l <=> 0 if $l != 0;
500 # lengths (corrected by exponent) are equal
501 # so make mantissa equal-length by padding with zero (shift left)
502 my $diff = $lxm - $lym;
503 my $xm = $x->{_m}; # not yet copy it
507 $ym = $y->{_m}->copy()->blsft($diff,10);
511 $xm = $x->{_m}->copy()->blsft(-$diff,10);
513 $xm->bacmp($ym) <=> 0;
518 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
519 # return result as BFLOAT
522 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
523 # objectify is costly, so avoid it
524 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
526 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
529 # inf and NaN handling
530 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
533 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
535 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
537 # +inf++inf or -inf+-inf => same, rest is NaN
538 return $x if $x->{sign} eq $y->{sign};
541 # +-inf + something => +inf; something +-inf => +-inf
542 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
546 return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
547 ((!$x->isa($self)) || (!$y->isa($self)));
549 # speed: no add for 0+y or x+0
550 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
551 if ($x->is_zero()) # 0+y
553 # make copy, clobbering up x (modify in place!)
554 $x->{_e} = $y->{_e}->copy();
555 $x->{_m} = $y->{_m}->copy();
556 $x->{sign} = $y->{sign} || $nan;
557 return $x->round($a,$p,$r,$y);
560 # take lower of the two e's and adapt m1 to it to match m2
562 $e = $MBI->bzero() if !defined $e; # if no BFLOAT ?
563 $e = $e->copy(); # make copy (didn't do it yet)
564 $e->bsub($x->{_e}); # Ye - Xe
565 my $add = $y->{_m}->copy();
566 if ($e->{sign} eq '-') # < 0
568 $x->{_e} += $e; # need the sign of e
569 $x->{_m}->blsft($e->babs(),10); # destroys copy of _e
571 elsif (!$e->is_zero()) # > 0
575 # else: both e are the same, so just leave them
576 $x->{_m}->{sign} = $x->{sign}; # fiddle with signs
577 $add->{sign} = $y->{sign};
578 $x->{_m} += $add; # finally do add/sub
579 $x->{sign} = $x->{_m}->{sign}; # re-adjust signs
580 $x->{_m}->{sign} = '+'; # mantissa always positiv
581 # delete trailing zeros, then round
582 $x->bnorm()->round($a,$p,$r,$y);
587 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
588 # subtract second arg from first, modify first
591 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
592 # objectify is costly, so avoid it
593 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
595 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
598 if ($y->is_zero()) # still round for not adding zero
600 return $x->round($a,$p,$r);
604 $y->{sign} =~ tr/+-/-+/; # does nothing for NaN
605 $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
606 $y->{sign} =~ tr/+-/-+/; # refix $y (does nothing for NaN)
607 $x; # already rounded by badd()
612 # increment arg by one
613 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
615 if ($x->{_e}->sign() eq '-')
617 return $x->badd($self->bone(),@r); # digits after dot
620 if (!$x->{_e}->is_zero()) # _e == 0 for NaN, inf, -inf
622 # 1e2 => 100, so after the shift below _m has a '0' as last digit
623 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
624 $x->{_e}->bzero(); # normalize
625 # we know that the last digit of $x will be '1' or '9', depending on the
629 if ($x->{sign} eq '+')
632 return $x->bnorm()->bround(@r);
634 elsif ($x->{sign} eq '-')
637 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
638 return $x->bnorm()->bround(@r);
640 # inf, nan handling etc
641 $x->badd($self->bone(),@r); # badd() does round
646 # decrement arg by one
647 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
649 if ($x->{_e}->sign() eq '-')
651 return $x->badd($self->bone('-'),@r); # digits after dot
654 if (!$x->{_e}->is_zero())
656 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
660 my $zero = $x->is_zero();
662 if (($x->{sign} eq '-') || $zero)
665 $x->{sign} = '-' if $zero; # 0 => 1 => -1
666 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
667 return $x->bnorm()->round(@r);
670 elsif ($x->{sign} eq '+')
673 return $x->bnorm()->round(@r);
675 # inf, nan handling etc
676 $x->badd($self->bone('-'),@r); # does round
683 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
685 # $base > 0, $base != 1; if $base == undef default to $base == e
688 # we need to limit the accuracy to protect against overflow
691 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
693 # also takes care of the "error in _find_round_parameters?" case
694 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
696 # no rounding at all, so must use fallback
697 if (scalar @params == 0)
699 # simulate old behaviour
700 $params[0] = $self->div_scale(); # and round to it as accuracy
701 $params[1] = undef; # P = undef
702 $scale = $params[0]+4; # at least four more for proper round
703 $params[2] = $r; # round mode by caller or undef
704 $fallback = 1; # to clear a/p afterwards
708 # the 4 below is empirical, and there might be cases where it is not
710 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
713 return $x->bzero(@params) if $x->is_one();
714 # base not defined => base == Euler's constant e
717 # make object, since we don't feed it through objectify() to still get the
718 # case of $base == undef
719 $base = $self->new($base) unless ref($base);
720 # $base > 0; $base != 1
721 return $x->bnan() if $base->is_zero() || $base->is_one() ||
722 $base->{sign} ne '+';
723 # if $x == $base, we know the result must be 1.0
724 return $x->bone('+',@params) if $x->bcmp($base) == 0;
727 # when user set globals, they would interfere with our calculation, so
728 # disable them and later re-enable them
730 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
731 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
732 # we also need to disable any set A or P on $x (_find_round_parameters took
733 # them already into account), since these would interfere, too
734 delete $x->{_a}; delete $x->{_p};
735 # need to disable $upgrade in BigInt, to avoid deep recursion
736 local $Math::BigInt::upgrade = undef;
737 local $Math::BigFloat::downgrade = undef;
739 # upgrade $x if $x is not a BigFloat (handle BigInt input)
740 if (!$x->isa('Math::BigFloat'))
742 $x = Math::BigFloat->new($x);
748 # If the base is defined and an integer, try to calculate integer result
749 # first. This is very fast, and in case the real result was found, we can
751 if (defined $base && $base->is_int() && $x->is_int())
753 my $int = $x->{_m}->copy();
754 $int->blsft($x->{_e},10) unless $x->{_e}->is_zero();
755 $int->blog($base->as_number());
757 if ($base->copy()->bpow($int) == $x)
759 # found result, return it
761 $x->{_e} = $MBI->bzero();
769 # first calculate the log to base e (using reduction by 10 (and probably 2))
770 $self->_log_10($x,$scale);
772 # and if a different base was requested, convert it
775 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
776 # not ln, but some other base (don't modify $base)
777 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
781 # shortcut to not run through _find_round_parameters again
782 if (defined $params[0])
784 $x->bround($params[0],$params[2]); # then round accordingly
788 $x->bfround($params[1],$params[2]); # then round accordingly
792 # clear a/p after round, since user did not request it
793 $x->{_a} = undef; $x->{_p} = undef;
796 $$abr = $ab; $$pbr = $pb;
803 # internal log function to calculate ln() based on Taylor series.
804 # Modifies $x in place.
805 my ($self,$x,$scale) = @_;
807 # in case of $x == 1, result is 0
808 return $x->bzero() if $x->is_one();
810 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
814 # Taylor: | u 1 u^3 1 u^5 |
815 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
816 # |_ v 3 v^3 5 v^5 _|
818 # This takes much more steps to calculate the result and is thus not used
821 # Taylor: | u 1 u^2 1 u^3 |
822 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
823 # |_ x 2 x^2 3 x^3 _|
825 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
827 $v = $x->copy(); $v->binc(); # v = x+1
828 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
829 $x->bdiv($v,$scale); # first term: u/v
832 $u *= $u; $v *= $v; # u^2, v^2
833 $below->bmul($v); # u^3, v^3
835 $factor = $self->new(3); $f = $self->new(2);
837 my $steps = 0 if DEBUG;
838 $limit = $self->new("1E-". ($scale-1));
841 # we calculate the next term, and add it to the last
842 # when the next term is below our limit, it won't affect the outcome
843 # anymore, so we stop
845 # calculating the next term simple from over/below will result in quite
846 # a time hog if the input has many digits, since over and below will
847 # accumulate more and more digits, and the result will also have many
848 # digits, but in the end it is rounded to $scale digits anyway. So if we
849 # round $over and $below first, we save a lot of time for the division
850 # (not with log(1.2345), but try log (123**123) to see what I mean. This
851 # can introduce a rounding error if the division result would be f.i.
852 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
853 # if we truncated $over and $below we might get 0.12345. Does this matter
854 # for the end result? So we give $over and $below 4 more digits to be
855 # on the safe side (unscientific error handling as usual... :+D
857 $next = $over->copy->bround($scale+4)->bdiv(
858 $below->copy->bmul($factor)->bround($scale+4),
862 ## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
864 last if $next->bacmp($limit) <= 0;
866 delete $next->{_a}; delete $next->{_p};
868 #print "step $x\n ($next - $limit = ",$next - $limit,")\n";
869 # calculate things for the next term
870 $over *= $u; $below *= $v; $factor->badd($f);
873 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
876 $x->bmul($f); # $x *= 2
877 print "took $steps steps\n" if DEBUG;
882 # Internal log function based on reducing input to the range of 0.1 .. 9.99
883 # and then "correcting" the result to the proper one. Modifies $x in place.
884 my ($self,$x,$scale) = @_;
886 # taking blog() from numbers greater than 10 takes a *very long* time, so we
887 # break the computation down into parts based on the observation that:
888 # blog(x*y) = blog(x) + blog(y)
889 # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
890 # the faster it get's, especially because 2*$x takes about 10 times as long,
891 # so by dividing $x by 10 we make it at least factor 100 faster...)
893 # The same observation is valid for numbers smaller than 0.1 (e.g. computing
894 # log(1) is fastest, and the farther away we get from 1, the longer it takes)
895 # so we also 'break' this down by multiplying $x with 10 and subtract the
896 # log(10) afterwards to get the correct result.
898 # calculate nr of digits before dot
899 my $dbd = $x->{_m}->length() + $x->{_e}->numify();
901 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
904 my $calc = 1; # do some calculation?
906 # disable the shortcut for 10, since we need log(10) and this would recurse
908 if ($x->{_e}->is_one() && $x->{_m}->is_one())
910 $dbd = 0; # disable shortcut
911 # we can use the cached value in these cases
912 if ($scale <= $LOG_10_A)
914 $x->bzero(); $x->badd($LOG_10);
915 $calc = 0; # no need to calc, but round
920 # disable the shortcut for 2, since we maybe have it cached
921 if ($x->{_e}->is_zero() && $x->{_m}->bcmp(2) == 0)
923 $dbd = 0; # disable shortcut
924 # we can use the cached value in these cases
925 if ($scale <= $LOG_2_A)
927 $x->bzero(); $x->badd($LOG_2);
928 $calc = 0; # no need to calc, but round
933 # if $x = 0.1, we know the result must be 0-log(10)
934 if ($calc != 0 && $x->{_e}->is_one('-') && $x->{_m}->is_one())
936 $dbd = 0; # disable shortcut
937 # we can use the cached value in these cases
938 if ($scale <= $LOG_10_A)
940 $x->bzero(); $x->bsub($LOG_10);
941 $calc = 0; # no need to calc, but round
945 return if $calc == 0; # already have the result
947 # default: these correction factors are undef and thus not used
948 my $l_10; # value of ln(10) to A of $scale
949 my $l_2; # value of ln(2) to A of $scale
951 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
952 # so don't do this shortcut for 1 or 0
953 if (($dbd > 1) || ($dbd < 0))
955 # convert our cached value to an object if not already (avoid doing this
956 # at import() time, since not everybody needs this)
957 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
959 #print "x = $x, dbd = $dbd, calc = $calc\n";
960 # got more than one digit before the dot, or more than one zero after the
962 # log(123) == log(1.23) + log(10) * 2
963 # log(0.0123) == log(1.23) - log(10) * 2
965 if ($scale <= $LOG_10_A)
968 #print "using cached value for l_10\n";
969 $l_10 = $LOG_10->copy(); # copy for mul
973 # else: slower, compute it (but don't cache it, because it could be big)
974 # also disable downgrade for this code path
975 local $Math::BigFloat::downgrade = undef;
976 #print "l_10 = $l_10 (self = $self',
977 # ", ref(l_10) = ",ref($l_10)," scale $scale)\n";
978 #print "calculating value for l_10, scale $scale\n";
979 $l_10 = $self->new(10)->blog(undef,$scale); # scale+4, actually
981 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
983 $dbd = $self->new($dbd);
985 $l_10->bmul($dbd); # log(10) * (digits_before_dot-1)
986 #print "l_10 = $l_10\n";
988 $x->{_e}->bsub($dbd); # 123 => 1.23
990 #print "calculating log($x) with scale=$scale\n";
994 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
996 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
997 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
999 my $half = $self->new('0.5');
1000 my $twos = 0; # default: none (0 times)
1001 my $two = $self->new(2);
1002 while ($x->bacmp($half) <= 0)
1004 $twos--; $x->bmul($two);
1006 while ($x->bacmp($two) >= 0)
1008 $twos++; $x->bdiv($two,$scale+4); # keep all digits
1011 # $twos > 0 => did mul 2, < 0 => did div 2 (never both)
1012 # calculate correction factor based on ln(2)
1015 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1016 if ($scale <= $LOG_2_A)
1019 #print "using cached value for l_10\n";
1020 $l_2 = $LOG_2->copy(); # copy for mul
1024 # else: slower, compute it (but don't cache it, because it could be big)
1025 # also disable downgrade for this code path
1026 local $Math::BigFloat::downgrade = undef;
1027 #print "calculating value for l_2, scale $scale\n";
1028 $l_2 = $two->blog(undef,$scale); # scale+4, actually
1030 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
1033 $self->_log($x,$scale); # need to do the "normal" way
1034 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
1035 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
1036 # all done, $x contains now the result
1041 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1042 # does not modify arguments, but returns new object
1043 # Lowest Common Multiplicator
1045 my ($self,@arg) = objectify(0,@_);
1046 my $x = $self->new(shift @arg);
1047 while (@arg) { $x = _lcm($x,shift @arg); }
1053 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
1054 # does not modify arguments, but returns new object
1055 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
1057 my ($self,@arg) = objectify(0,@_);
1058 my $x = $self->new(shift @arg);
1059 while (@arg) { $x = _gcd($x,shift @arg); }
1063 ###############################################################################
1064 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1068 # internal, return true if BigInt arg is zero or one, saving the
1069 # two calls to is_zero() and is_one()
1072 $x->{sign} eq '+' && ($x->is_zero() || $x->is_one());
1077 # return true if arg (BFLOAT or num_str) is an integer
1078 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1080 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1081 $x->{_e}->{sign} eq '+'; # 1e-1 => no integer
1087 # return true if arg (BFLOAT or num_str) is zero
1088 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1090 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
1096 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1097 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1099 $sign = '+' if !defined $sign || $sign ne '-';
1101 if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
1107 # return true if arg (BFLOAT or num_str) is odd or false if even
1108 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1110 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1111 ($x->{_e}->is_zero() && $x->{_m}->is_odd());
1117 # return true if arg (BINT or num_str) is even or false if odd
1118 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1120 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
1121 return 1 if ($x->{_e}->{sign} eq '+' # 123.45 is never
1122 && $x->{_m}->is_even()); # but 1200 is
1128 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
1129 # (BINT or num_str, BINT or num_str) return BINT
1132 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1133 # objectify is costly, so avoid it
1134 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1136 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1139 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1142 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1144 return $x->bnan() if $x->is_zero() || $y->is_zero();
1145 # result will always be +-inf:
1146 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1147 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1148 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1149 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1150 return $x->binf('-');
1153 return $x->bzero() if $x->is_zero() || $y->is_zero();
1155 return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
1156 ((!$x->isa($self)) || (!$y->isa($self)));
1158 # aEb * cEd = (a*c)E(b+d)
1159 $x->{_m}->bmul($y->{_m});
1160 $x->{_e}->badd($y->{_e});
1162 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1163 return $x->bnorm()->round($a,$p,$r,$y);
1168 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1169 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1172 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1173 # objectify is costly, so avoid it
1174 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1176 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1179 return $self->_div_inf($x,$y)
1180 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1182 # x== 0 # also: or y == 1 or y == -1
1183 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1186 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1188 # we need to limit the accuracy to protect against overflow
1190 my (@params,$scale);
1191 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1193 return $x if $x->is_nan(); # error in _find_round_parameters?
1195 # no rounding at all, so must use fallback
1196 if (scalar @params == 0)
1198 # simulate old behaviour
1199 $params[0] = $self->div_scale(); # and round to it as accuracy
1200 $scale = $params[0]+4; # at least four more for proper round
1201 $params[2] = $r; # round mode by caller or undef
1202 $fallback = 1; # to clear a/p afterwards
1206 # the 4 below is empirical, and there might be cases where it is not
1208 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1210 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
1211 $scale = $lx if $lx > $scale;
1212 $scale = $ly if $ly > $scale;
1213 my $diff = $ly - $lx;
1214 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
1216 # make copy of $x in case of list context for later reminder calculation
1218 if (wantarray && !$y->is_one())
1223 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1225 # check for / +-1 ( +/- 1E0)
1228 # promote BigInts and it's subclasses (except when already a BigFloat)
1229 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1231 # need to disable $upgrade in BigInt, to avoid deep recursion
1232 local $Math::BigInt::upgrade = undef; # should be parent class vs MBI
1234 # calculate the result to $scale digits and then round it
1235 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1236 $x->{_m}->blsft($scale,10);
1237 $x->{_m}->bdiv( $y->{_m} ); # a/c
1238 $x->{_e}->bsub( $y->{_e} ); # b-d
1239 $x->{_e}->bsub($scale); # correct for 10**scale
1240 $x->bnorm(); # remove trailing 0's
1243 # shortcut to not run through _find_round_parameters again
1244 if (defined $params[0])
1246 $x->{_a} = undef; # clear before round
1247 $x->bround($params[0],$params[2]); # then round accordingly
1251 $x->{_p} = undef; # clear before round
1252 $x->bfround($params[1],$params[2]); # then round accordingly
1256 # clear a/p after round, since user did not request it
1257 $x->{_a} = undef; $x->{_p} = undef;
1264 $rem->bmod($y,@params); # copy already done
1268 $rem = $self->bzero();
1272 # clear a/p after round, since user did not request it
1273 $rem->{_a} = undef; $rem->{_p} = undef;
1282 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1285 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1286 # objectify is costly, so avoid it
1287 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1289 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1292 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1294 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1295 $x->{sign} = $re->{sign};
1296 $x->{_e} = $re->{_e};
1297 $x->{_m} = $re->{_m};
1298 return $x->round($a,$p,$r,$y);
1300 return $x->bnan() if $x->is_zero() && $y->is_zero();
1301 return $x if $y->is_zero();
1302 return $x->bnan() if $x->is_nan() || $y->is_nan();
1303 return $x->bzero() if $y->is_one() || $x->is_zero();
1305 # inf handling is missing here
1307 my $cmp = $x->bacmp($y); # equal or $x < $y?
1308 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1310 # only $y of the operands negative?
1311 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1313 $x->{sign} = $y->{sign}; # calc sign first
1314 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1316 my $ym = $y->{_m}->copy();
1319 $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
1321 # if $y has digits after dot
1322 my $shifty = 0; # correct _e of $x by this
1323 if ($y->{_e}->{sign} eq '-') # has digits after dot
1325 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1326 $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
1327 $x->blsft($shifty,10); # 123 => 1230, $y->{_m} is already 25
1329 # $ym is now mantissa of $y based on exponent 0
1331 my $shiftx = 0; # correct _e of $x by this
1332 if ($x->{_e}->{sign} eq '-') # has digits after dot
1334 # 123.4 % 20 => 1234 % 200
1335 $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
1336 $ym->blsft($shiftx,10);
1338 # 123e1 % 20 => 1230 % 20
1339 if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
1341 $x->{_m}->blsft($x->{_e},10);
1343 $x->{_e} = $MBI->bzero() unless $x->{_e}->is_zero();
1345 $x->{_e}->bsub($shiftx) if $shiftx != 0;
1346 $x->{_e}->bsub($shifty) if $shifty != 0;
1348 # now mantissas are equalized, exponent of $x is adjusted, so calc result
1350 $x->{_m}->bmod($ym);
1352 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
1355 if ($neg != 0) # one of them negative => correct in place
1358 $x->{_m} = $r->{_m};
1359 $x->{_e} = $r->{_e};
1360 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
1364 $x->round($a,$p,$r,$y); # round and return
1369 # calculate $y'th root of $x
1372 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1373 # objectify is costly, so avoid it
1374 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1376 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1379 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1380 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1381 $y->{sign} !~ /^\+$/;
1383 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1385 # we need to limit the accuracy to protect against overflow
1387 my (@params,$scale);
1388 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1390 return $x if $x->is_nan(); # error in _find_round_parameters?
1392 # no rounding at all, so must use fallback
1393 if (scalar @params == 0)
1395 # simulate old behaviour
1396 $params[0] = $self->div_scale(); # and round to it as accuracy
1397 $scale = $params[0]+4; # at least four more for proper round
1398 $params[2] = $r; # round mode by caller or undef
1399 $fallback = 1; # to clear a/p afterwards
1403 # the 4 below is empirical, and there might be cases where it is not
1405 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1408 # when user set globals, they would interfere with our calculation, so
1409 # disable them and later re-enable them
1411 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1412 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1413 # we also need to disable any set A or P on $x (_find_round_parameters took
1414 # them already into account), since these would interfere, too
1415 delete $x->{_a}; delete $x->{_p};
1416 # need to disable $upgrade in BigInt, to avoid deep recursion
1417 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1419 # remember sign and make $x positive, since -4 ** (1/2) => -2
1420 my $sign = 0; $sign = 1 if $x->is_negative(); $x->babs();
1422 if ($y->bcmp(2) == 0) # normal square root
1424 $x->bsqrt($scale+4);
1426 elsif ($y->is_one('-'))
1429 my $u = $self->bone()->bdiv($x,$scale);
1430 # copy private parts over
1431 $x->{_m} = $u->{_m};
1432 $x->{_e} = $u->{_e};
1436 # calculate the broot() as integer result first, and if it fits, return
1437 # it rightaway (but only if $x and $y are integer):
1439 my $done = 0; # not yet
1440 if ($y->is_int() && $x->is_int())
1442 my $int = $x->{_m}->copy();
1443 $int->blsft($x->{_e},10) unless $x->{_e}->is_zero();
1444 $int->broot($y->as_number());
1446 if ($int->copy()->bpow($y) == $x)
1448 # found result, return it
1450 $x->{_e} = $MBI->bzero();
1457 my $u = $self->bone()->bdiv($y,$scale+4);
1458 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
1459 $x->bpow($u,$scale+4); # el cheapo
1462 $x->bneg() if $sign == 1;
1464 # shortcut to not run through _find_round_parameters again
1465 if (defined $params[0])
1467 $x->bround($params[0],$params[2]); # then round accordingly
1471 $x->bfround($params[1],$params[2]); # then round accordingly
1475 # clear a/p after round, since user did not request it
1476 $x->{_a} = undef; $x->{_p} = undef;
1479 $$abr = $ab; $$pbr = $pb;
1485 # calculate square root
1486 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1488 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
1489 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
1490 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
1492 # we need to limit the accuracy to protect against overflow
1494 my (@params,$scale);
1495 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1497 return $x if $x->is_nan(); # error in _find_round_parameters?
1499 # no rounding at all, so must use fallback
1500 if (scalar @params == 0)
1502 # simulate old behaviour
1503 $params[0] = $self->div_scale(); # and round to it as accuracy
1504 $scale = $params[0]+4; # at least four more for proper round
1505 $params[2] = $r; # round mode by caller or undef
1506 $fallback = 1; # to clear a/p afterwards
1510 # the 4 below is empirical, and there might be cases where it is not
1512 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1515 # when user set globals, they would interfere with our calculation, so
1516 # disable them and later re-enable them
1518 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1519 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1520 # we also need to disable any set A or P on $x (_find_round_parameters took
1521 # them already into account), since these would interfere, too
1522 delete $x->{_a}; delete $x->{_p};
1523 # need to disable $upgrade in BigInt, to avoid deep recursion
1524 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1526 my $xas = $x->as_number();
1527 my $gs = $xas->copy()->bsqrt(); # some guess
1529 if (($x->{_e}->{sign} ne '-') # guess can't be accurate if there are
1530 # digits after the dot
1531 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
1534 $x->{_m} = $gs; $x->{_e} = $MBI->bzero(); $x->bnorm();
1535 # shortcut to not run through _find_round_parameters again
1536 if (defined $params[0])
1538 $x->bround($params[0],$params[2]); # then round accordingly
1542 $x->bfround($params[1],$params[2]); # then round accordingly
1546 # clear a/p after round, since user did not request it
1547 $x->{_a} = undef; $x->{_p} = undef;
1549 # re-enable A and P, upgrade is taken care of by "local"
1550 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1554 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1555 # of the result by multipyling the input by 100 and then divide the integer
1556 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
1557 # this will transform 123.456 (in $x) into 123456 (in $y1)
1558 my $y1 = $x->{_m}->copy();
1559 # We now make sure that $y1 has the same odd or even number of digits than
1560 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1561 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1562 # steps of 10. The length of $x does not count, since an even or odd number
1563 # of digits before the dot is not changed by adding an even number of digits
1564 # after the dot (the result is still odd or even digits long).
1565 my $length = $y1->length();
1566 $y1->bmul(10) if $x->{_e}->is_odd();
1567 # now calculate how many digits the result of sqrt(y1) would have
1568 my $digits = int($length / 2);
1569 # but we need at least $scale digits, so calculate how many are missing
1570 my $shift = $scale - $digits;
1571 # that should never happen (we take care of integer guesses above)
1572 # $shift = 0 if $shift < 0;
1573 # multiply in steps of 100, by shifting left two times the "missing" digits
1574 $y1->blsft($shift*2,10);
1575 # now take the square root and truncate to integer
1577 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1578 # result, which is than later rounded to the desired scale.
1580 # calculate how many zeros $x had after the '.' (or before it, depending
1581 # on sign of $dat, the result should have half as many:
1582 my $dat = $length + $x->{_e}->numify();
1586 # no zeros after the dot (e.g. 1.23, 0.49 etc)
1587 # preserve half as many digits before the dot than the input had
1588 # (but round this "up")
1589 $dat = int(($dat+1)/2);
1593 $dat = int(($dat)/2);
1595 $x->{_e}= $MBI->new( $dat - $y1->length() );
1599 # shortcut to not run through _find_round_parameters again
1600 if (defined $params[0])
1602 $x->bround($params[0],$params[2]); # then round accordingly
1606 $x->bfround($params[1],$params[2]); # then round accordingly
1610 # clear a/p after round, since user did not request it
1611 $x->{_a} = undef; $x->{_p} = undef;
1614 $$abr = $ab; $$pbr = $pb;
1620 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1621 # compute factorial number, modifies first argument
1624 my ($self,$x,@r) = (ref($_[0]),@_);
1625 # objectify is costly, so avoid it
1626 ($self,$x,@r) = objectify(1,@_) if !ref($x);
1628 return $x if $x->{sign} eq '+inf'; # inf => inf
1630 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1631 ($x->{_e}->{sign} ne '+')); # digits after dot?
1633 # use BigInt's bfac() for faster calc
1634 if (! $x->{_e}->is_zero())
1636 $x->{_m}->blsft($x->{_e},10); # change 12e1 to 120e0
1639 $x->{_m}->bfac(); # calculate factorial
1640 $x->bnorm()->round(@r); # norm again and round result
1645 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1646 my ($x,$y,$a,$p,$r) = @_;
1649 # if $y == 0.5, it is sqrt($x)
1650 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1653 # a ** x == e ** (x * ln a)
1657 # Taylor: | u u^2 u^3 |
1658 # x ** y = 1 + | --- + --- + ----- + ... |
1661 # we need to limit the accuracy to protect against overflow
1663 my ($scale,@params);
1664 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1666 return $x if $x->is_nan(); # error in _find_round_parameters?
1668 # no rounding at all, so must use fallback
1669 if (scalar @params == 0)
1671 # simulate old behaviour
1672 $params[0] = $self->div_scale(); # and round to it as accuracy
1673 $params[1] = undef; # disable P
1674 $scale = $params[0]+4; # at least four more for proper round
1675 $params[2] = $r; # round mode by caller or undef
1676 $fallback = 1; # to clear a/p afterwards
1680 # the 4 below is empirical, and there might be cases where it is not
1682 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1685 # when user set globals, they would interfere with our calculation, so
1686 # disable them and later re-enable them
1688 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1689 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1690 # we also need to disable any set A or P on $x (_find_round_parameters took
1691 # them already into account), since these would interfere, too
1692 delete $x->{_a}; delete $x->{_p};
1693 # need to disable $upgrade in BigInt, to avoid deep recursion
1694 local $Math::BigInt::upgrade = undef;
1696 my ($limit,$v,$u,$below,$factor,$next,$over);
1698 $u = $x->copy()->blog(undef,$scale)->bmul($y);
1699 $v = $self->bone(); # 1
1700 $factor = $self->new(2); # 2
1701 $x->bone(); # first term: 1
1703 $below = $v->copy();
1706 $limit = $self->new("1E-". ($scale-1));
1710 # we calculate the next term, and add it to the last
1711 # when the next term is below our limit, it won't affect the outcome
1712 # anymore, so we stop
1713 $next = $over->copy()->bdiv($below,$scale);
1714 last if $next->bacmp($limit) <= 0;
1716 # calculate things for the next term
1717 $over *= $u; $below *= $factor; $factor->binc();
1721 # shortcut to not run through _find_round_parameters again
1722 if (defined $params[0])
1724 $x->bround($params[0],$params[2]); # then round accordingly
1728 $x->bfround($params[1],$params[2]); # then round accordingly
1732 # clear a/p after round, since user did not request it
1733 $x->{_a} = undef; $x->{_p} = undef;
1736 $$abr = $ab; $$pbr = $pb;
1742 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1743 # compute power of two numbers, second arg is used as integer
1744 # modifies first argument
1747 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1748 # objectify is costly, so avoid it
1749 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1751 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1754 return $x if $x->{sign} =~ /^[+-]inf$/;
1755 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1756 return $x->bone() if $y->is_zero();
1757 return $x if $x->is_one() || $y->is_one();
1759 return $x->_pow($y,$a,$p,$r) if !$y->is_int(); # non-integer power
1761 my $y1 = $y->as_number(); # make bigint
1763 if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1765 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
1766 return $y1->is_odd() ? $x : $x->babs(1);
1770 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1771 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1775 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1777 $x->{_m}->bpow($y1);
1778 $x->{_e}->bmul($y1);
1779 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1781 if ($y->{sign} eq '-')
1783 # modify $x in place!
1784 my $z = $x->copy(); $x->bzero()->binc();
1785 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
1787 $x->round($a,$p,$r,$y);
1790 ###############################################################################
1791 # rounding functions
1795 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1796 # $n == 0 means round to integer
1797 # expects and returns normalized numbers!
1798 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1800 return $x if $x->modify('bfround');
1802 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1803 return $x if !defined $scale; # no-op
1805 # never round a 0, +-inf, NaN
1808 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1811 return $x if $x->{sign} !~ /^[+-]$/;
1813 # don't round if x already has lower precision
1814 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1816 $x->{_p} = $scale; # remember round in any case
1817 $x->{_a} = undef; # and clear A
1820 # round right from the '.'
1822 return $x if $x->{_e}->{sign} eq '+'; # e >= 0 => nothing to round
1824 $scale = -$scale; # positive for simplicity
1825 my $len = $x->{_m}->length(); # length of mantissa
1827 # the following poses a restriction on _e, but if _e is bigger than a
1828 # scalar, you got other problems (memory etc) anyway
1829 my $dad = -($x->{_e}->numify()); # digits after dot
1830 my $zad = 0; # zeros after dot
1831 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
1833 #print "scale $scale dad $dad zad $zad len $len\n";
1834 # number bsstr len zad dad
1835 # 0.123 123e-3 3 0 3
1836 # 0.0123 123e-4 3 1 4
1839 # 1.2345 12345e-4 5 0 4
1841 # do not round after/right of the $dad
1842 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1844 # round to zero if rounding inside the $zad, but not for last zero like:
1845 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1846 return $x->bzero() if $scale < $zad;
1847 if ($scale == $zad) # for 0.006, scale -3 and trunc
1853 # adjust round-point to be inside mantissa
1856 $scale = $scale-$zad;
1860 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
1861 $scale = $dbd+$scale;
1867 # round left from the '.'
1869 # 123 => 100 means length(123) = 3 - $scale (2) => 1
1871 my $dbt = $x->{_m}->length();
1873 my $dbd = $dbt + $x->{_e}->numify();
1874 # should be the same, so treat it as this
1875 $scale = 1 if $scale == 0;
1876 # shortcut if already integer
1877 return $x if $scale == 1 && $dbt <= $dbd;
1878 # maximum digits before dot
1883 # not enough digits before dot, so round to zero
1886 elsif ( $scale == $dbd )
1893 $scale = $dbd - $scale;
1896 # pass sign to bround for rounding modes '+inf' and '-inf'
1897 $x->{_m}->{sign} = $x->{sign};
1898 $x->{_m}->bround($scale,$mode);
1899 $x->{_m}->{sign} = '+'; # fix sign back
1905 # accuracy: preserve $N digits, and overwrite the rest with 0's
1906 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1908 if (($_[0] || 0) < 0)
1910 require Carp; Carp::croak ('bround() needs positive accuracy');
1913 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1914 return $x if !defined $scale; # no-op
1916 return $x if $x->modify('bround');
1918 # scale is now either $x->{_a}, $accuracy, or the user parameter
1919 # test whether $x already has lower accuracy, do nothing in this case
1920 # but do round if the accuracy is the same, since a math operation might
1921 # want to round a number with A=5 to 5 digits afterwards again
1922 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1924 # scale < 0 makes no sense
1925 # never round a +-inf, NaN
1926 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
1928 # 1: $scale == 0 => keep all digits
1929 # 2: never round a 0
1930 # 3: if we should keep more digits than the mantissa has, do nothing
1931 if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1933 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1937 # pass sign to bround for '+inf' and '-inf' rounding modes
1938 $x->{_m}->{sign} = $x->{sign};
1939 $x->{_m}->bround($scale,$mode); # round mantissa
1940 $x->{_m}->{sign} = '+'; # fix sign back
1941 # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1942 $x->{_a} = $scale; # remember rounding
1943 $x->{_p} = undef; # and clear P
1944 $x->bnorm(); # del trailing zeros gen. by bround()
1949 # return integer less or equal then $x
1950 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1952 return $x if $x->modify('bfloor');
1954 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1956 # if $x has digits after dot
1957 if ($x->{_e}->{sign} eq '-')
1959 $x->{_e}->{sign} = '+'; # negate e
1960 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1961 $x->{_e}->bzero(); # trunc/norm
1962 $x->{_m}->binc() if $x->{sign} eq '-'; # decrement if negative
1964 $x->round($a,$p,$r);
1969 # return integer greater or equal then $x
1970 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1972 return $x if $x->modify('bceil');
1973 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1975 # if $x has digits after dot
1976 if ($x->{_e}->{sign} eq '-')
1978 #$x->{_m}->brsft(-$x->{_e},10);
1980 #$x++ if $x->{sign} eq '+';
1982 $x->{_e}->{sign} = '+'; # negate e
1983 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1984 $x->{_e}->bzero(); # trunc/norm
1985 $x->{_m}->binc() if $x->{sign} eq '+'; # decrement if negative
1987 $x->round($a,$p,$r);
1992 # shift right by $y (divide by power of $n)
1995 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1996 # objectify is costly, so avoid it
1997 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1999 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2002 return $x if $x->modify('brsft');
2003 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2005 $n = 2 if !defined $n; $n = $self->new($n);
2006 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
2011 # shift left by $y (multiply by power of $n)
2014 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2015 # objectify is costly, so avoid it
2016 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2018 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2021 return $x if $x->modify('blsft');
2022 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2024 $n = 2 if !defined $n; $n = $self->new($n);
2025 $x->bmul($n->bpow($y),$a,$p,$r,$y);
2028 ###############################################################################
2032 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
2037 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
2038 # or falling back to MBI::bxxx()
2039 my $name = $AUTOLOAD;
2041 $name =~ s/.*:://; # split package
2043 $class->import() if $IMPORT == 0;
2044 if (!method_alias($name))
2048 # delayed load of Carp and avoid recursion
2050 Carp::croak ("Can't call a method without name");
2052 if (!method_hand_up($name))
2054 # delayed load of Carp and avoid recursion
2056 Carp::croak ("Can't call $class\-\>$name, not a valid method");
2058 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
2060 return &{"$MBI"."::$name"}(@_);
2062 my $bname = $name; $bname =~ s/^f/b/;
2063 *{$class."::$name"} = \&$bname;
2069 # return a copy of the exponent
2070 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2072 if ($x->{sign} !~ /^[+-]$/)
2074 my $s = $x->{sign}; $s =~ s/^[+-]//;
2075 return $self->new($s); # -inf, +inf => +inf
2077 return $x->{_e}->copy();
2082 # return a copy of the mantissa
2083 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2085 if ($x->{sign} !~ /^[+-]$/)
2087 my $s = $x->{sign}; $s =~ s/^[+]//;
2088 return $self->new($s); # -inf, +inf => +inf
2090 my $m = $x->{_m}->copy(); # faster than going via bstr()
2091 $m->bneg() if $x->{sign} eq '-';
2098 # return a copy of both the exponent and the mantissa
2099 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2101 if ($x->{sign} !~ /^[+-]$/)
2103 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2104 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2106 my $m = $x->{_m}->copy(); # faster than going via bstr()
2107 $m->bneg() if $x->{sign} eq '-';
2108 return ($m,$x->{_e}->copy());
2111 ##############################################################################
2112 # private stuff (internal use only)
2118 my $lib = ''; my @a;
2120 for ( my $i = 0; $i < $l ; $i++)
2122 if ( $_[$i] eq ':constant' )
2124 # This causes overlord er load to step in. 'binary' and 'integer'
2125 # are handled by BigInt.
2126 overload::constant float => sub { $self->new(shift); };
2128 elsif ($_[$i] eq 'upgrade')
2130 # this causes upgrading
2131 $upgrade = $_[$i+1]; # or undef to disable
2134 elsif ($_[$i] eq 'downgrade')
2136 # this causes downgrading
2137 $downgrade = $_[$i+1]; # or undef to disable
2140 elsif ($_[$i] eq 'lib')
2142 # alternative library
2143 $lib = $_[$i+1] || ''; # default Calc
2146 elsif ($_[$i] eq 'with')
2148 # alternative class for our private parts()
2149 $MBI = $_[$i+1] || 'Math::BigInt'; # default Math::BigInt
2158 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2159 my $mbilib = eval { Math::BigInt->config()->{lib} };
2160 if ((defined $mbilib) && ($MBI eq 'Math::BigInt'))
2162 # MBI already loaded
2163 $MBI->import('lib',"$lib,$mbilib", 'objectify');
2167 # MBI not loaded, or with ne "Math::BigInt"
2168 $lib .= ",$mbilib" if defined $mbilib;
2169 $lib =~ s/^,//; # don't leave empty
2170 # replacement library can handle lib statement, but also could ignore it
2173 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2174 # used in the same script, or eval inside import().
2175 my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
2176 my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
2178 $file = File::Spec->catfile (@parts, $file);
2179 eval { require "$file"; };
2180 $MBI->import( lib => $lib, 'objectify' );
2184 my $rc = "use $MBI lib => '$lib', 'objectify';";
2190 require Carp; Carp::croak ("Couldn't load $MBI: $! $@");
2193 # any non :constant stuff is handled by our parent, Exporter
2194 # even if @_ is empty, to give it a chance
2195 $self->SUPER::import(@a); # for subclasses
2196 $self->export_to_level(1,$self,@a); # need this, too
2201 # adjust m and e so that m is smallest possible
2202 # round number according to accuracy and precision settings
2203 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2205 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2207 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
2210 my $z = $MBI->new($zeros,undef,undef);
2211 $x->{_m}->brsft($z,10); $x->{_e}->badd($z);
2215 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
2216 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
2217 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
2220 # this is to prevent automatically rounding when MBI's globals are set
2221 $x->{_m}->{_f} = MB_NEVER_ROUND;
2222 $x->{_e}->{_f} = MB_NEVER_ROUND;
2223 # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
2224 $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
2225 $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
2226 $x; # MBI bnorm is no-op, so dont call it
2229 ##############################################################################
2233 # return number as hexadecimal string (only for integers defined)
2234 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2236 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2237 return '0x0' if $x->is_zero();
2239 return $nan if $x->{_e}->{sign} ne '+'; # how to do 1e-1 in hex!?
2241 my $z = $x->{_m}->copy();
2242 if (!$x->{_e}->is_zero()) # > 0
2244 $z->blsft($x->{_e},10);
2246 $z->{sign} = $x->{sign};
2252 # return number as binary digit string (only for integers defined)
2253 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2255 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2256 return '0b0' if $x->is_zero();
2258 return $nan if $x->{_e}->{sign} ne '+'; # how to do 1e-1 in hex!?
2260 my $z = $x->{_m}->copy();
2261 if (!$x->{_e}->is_zero()) # > 0
2263 $z->blsft($x->{_e},10);
2265 $z->{sign} = $x->{sign};
2271 # return copy as a bigint representation of this BigFloat number
2272 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2274 my $z = $x->{_m}->copy();
2275 if ($x->{_e}->{sign} eq '-') # < 0
2277 $x->{_e}->{sign} = '+'; # flip
2278 $z->brsft($x->{_e},10);
2279 $x->{_e}->{sign} = '-'; # flip back
2281 elsif (!$x->{_e}->is_zero()) # > 0
2283 $z->blsft($x->{_e},10);
2285 $z->{sign} = $x->{sign};
2292 my $class = ref($x) || $x;
2293 $x = $class->new(shift) unless ref($x);
2295 return 1 if $x->{_m}->is_zero();
2296 my $len = $x->{_m}->length();
2297 $len += $x->{_e} if $x->{_e}->sign() eq '+';
2300 my $t = $MBI->bzero();
2301 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
2312 Math::BigFloat - Arbitrary size floating point math package
2319 $x = Math::BigFloat->new($str); # defaults to 0
2320 $nan = Math::BigFloat->bnan(); # create a NotANumber
2321 $zero = Math::BigFloat->bzero(); # create a +0
2322 $inf = Math::BigFloat->binf(); # create a +inf
2323 $inf = Math::BigFloat->binf('-'); # create a -inf
2324 $one = Math::BigFloat->bone(); # create a +1
2325 $one = Math::BigFloat->bone('-'); # create a -1
2328 $x->is_zero(); # true if arg is +0
2329 $x->is_nan(); # true if arg is NaN
2330 $x->is_one(); # true if arg is +1
2331 $x->is_one('-'); # true if arg is -1
2332 $x->is_odd(); # true if odd, false for even
2333 $x->is_even(); # true if even, false for odd
2334 $x->is_positive(); # true if >= 0
2335 $x->is_negative(); # true if < 0
2336 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2338 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2339 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2340 $x->sign(); # return the sign, either +,- or NaN
2341 $x->digit($n); # return the nth digit, counting from right
2342 $x->digit(-$n); # return the nth digit, counting from left
2344 # The following all modify their first argument. If you want to preserve
2345 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2346 # neccessary when mixing $a = $b assigments with non-overloaded math.
2349 $x->bzero(); # set $i to 0
2350 $x->bnan(); # set $i to NaN
2351 $x->bone(); # set $x to +1
2352 $x->bone('-'); # set $x to -1
2353 $x->binf(); # set $x to inf
2354 $x->binf('-'); # set $x to -inf
2356 $x->bneg(); # negation
2357 $x->babs(); # absolute value
2358 $x->bnorm(); # normalize (no-op)
2359 $x->bnot(); # two's complement (bit wise not)
2360 $x->binc(); # increment x by 1
2361 $x->bdec(); # decrement x by 1
2363 $x->badd($y); # addition (add $y to $x)
2364 $x->bsub($y); # subtraction (subtract $y from $x)
2365 $x->bmul($y); # multiplication (multiply $x by $y)
2366 $x->bdiv($y); # divide, set $x to quotient
2367 # return (quo,rem) or quo if scalar
2369 $x->bmod($y); # modulus ($x % $y)
2370 $x->bpow($y); # power of arguments ($x ** $y)
2371 $x->blsft($y); # left shift
2372 $x->brsft($y); # right shift
2373 # return (quo,rem) or quo if scalar
2375 $x->blog(); # logarithm of $x to base e (Euler's number)
2376 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
2378 $x->band($y); # bit-wise and
2379 $x->bior($y); # bit-wise inclusive or
2380 $x->bxor($y); # bit-wise exclusive or
2381 $x->bnot(); # bit-wise not (two's complement)
2383 $x->bsqrt(); # calculate square-root
2384 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
2385 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2387 $x->bround($N); # accuracy: preserve $N digits
2388 $x->bfround($N); # precision: round to the $Nth digit
2390 $x->bfloor(); # return integer less or equal than $x
2391 $x->bceil(); # return integer greater or equal than $x
2393 # The following do not modify their arguments:
2395 bgcd(@values); # greatest common divisor
2396 blcm(@values); # lowest common multiplicator
2398 $x->bstr(); # return string
2399 $x->bsstr(); # return string in scientific notation
2401 $x->exponent(); # return exponent as BigInt
2402 $x->mantissa(); # return mantissa as BigInt
2403 $x->parts(); # return (mantissa,exponent) as BigInt
2405 $x->length(); # number of digits (w/o sign and '.')
2406 ($l,$f) = $x->length(); # number of digits, and length of fraction
2408 $x->precision(); # return P of $x (or global, if P of $x undef)
2409 $x->precision($n); # set P of $x to $n
2410 $x->accuracy(); # return A of $x (or global, if A of $x undef)
2411 $x->accuracy($n); # set A $x to $n
2413 # these get/set the appropriate global value for all BigFloat objects
2414 Math::BigFloat->precision(); # Precision
2415 Math::BigFloat->accuracy(); # Accuracy
2416 Math::BigFloat->round_mode(); # rounding mode
2420 All operators (inlcuding basic math operations) are overloaded if you
2421 declare your big floating point numbers as
2423 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2425 Operations with overloaded operators preserve the arguments, which is
2426 exactly what you expect.
2428 =head2 Canonical notation
2430 Input to these routines are either BigFloat objects, or strings of the
2431 following four forms:
2445 C</^[+-]\d+E[+-]?\d+$/>
2449 C</^[+-]\d*\.\d+E[+-]?\d+$/>
2453 all with optional leading and trailing zeros and/or spaces. Additonally,
2454 numbers are allowed to have an underscore between any two digits.
2456 Empty strings as well as other illegal numbers results in 'NaN'.
2458 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2459 are always stored in normalized form. On a string, it creates a BigFloat
2464 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2466 The string output will always have leading and trailing zeros stripped and drop
2467 a plus sign. C<bstr()> will give you always the form with a decimal point,
2468 while C<bsstr()> (s for scientific) gives you the scientific notation.
2470 Input bstr() bsstr()
2472 ' -123 123 123' '-123123123' '-123123123E0'
2473 '00.0123' '0.0123' '123E-4'
2474 '123.45E-2' '1.2345' '12345E-4'
2475 '10E+3' '10000' '1E4'
2477 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2478 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2479 return either undef, <0, 0 or >0 and are suited for sort.
2481 Actual math is done by using the class defined with C<with => Class;> (which
2482 defaults to BigInts) to represent the mantissa and exponent.
2484 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2485 represent the result when input arguments are not numbers, as well as
2486 the result of dividing by zero.
2488 =head2 C<mantissa()>, C<exponent()> and C<parts()>
2490 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2491 as BigInts such that:
2493 $m = $x->mantissa();
2494 $e = $x->exponent();
2495 $y = $m * ( 10 ** $e );
2496 print "ok\n" if $x == $y;
2498 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2500 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2502 Currently the mantissa is reduced as much as possible, favouring higher
2503 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2504 This might change in the future, so do not depend on it.
2506 =head2 Accuracy vs. Precision
2508 See also: L<Rounding|Rounding>.
2510 Math::BigFloat supports both precision and accuracy. For a full documentation,
2511 examples and tips on these topics please see the large section in
2514 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2515 a operation consumes all resources, each operation produces no more than
2516 the requested number of digits.
2518 Please refer to BigInt's documentation for the precedence rules of which
2519 accuracy/precision setting will be used.
2521 If there is no gloabl precision set, B<and> the operation inquestion was not
2522 called with a requested precision or accuracy, B<and> the input $x has no
2523 accuracy or precision set, then a fallback parameter will be used. For
2524 historical reasons, it is called C<div_scale> and can be accessed via:
2526 $d = Math::BigFloat->div_scale(); # query
2527 Math::BigFloat->div_scale($n); # set to $n digits
2529 The default value is 40 digits.
2531 In case the result of one operation has more precision than specified,
2532 it is rounded. The rounding mode taken is either the default mode, or the one
2533 supplied to the operation after the I<scale>:
2535 $x = Math::BigFloat->new(2);
2536 Math::BigFloat->precision(5); # 5 digits max
2537 $y = $x->copy()->bdiv(3); # will give 0.66666
2538 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2539 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
2540 Math::BigFloat->round_mode('zero');
2541 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2547 =item ffround ( +$scale )
2549 Rounds to the $scale'th place left from the '.', counting from the dot.
2550 The first digit is numbered 1.
2552 =item ffround ( -$scale )
2554 Rounds to the $scale'th place right from the '.', counting from the dot.
2558 Rounds to an integer.
2560 =item fround ( +$scale )
2562 Preserves accuracy to $scale digits from the left (aka significant digits)
2563 and pads the rest with zeros. If the number is between 1 and -1, the
2564 significant digits count from the first non-zero after the '.'
2566 =item fround ( -$scale ) and fround ( 0 )
2568 These are effectively no-ops.
2572 All rounding functions take as a second parameter a rounding mode from one of
2573 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2575 The default rounding mode is 'even'. By using
2576 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
2577 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2578 no longer supported.
2579 The second parameter to the round functions then overrides the default
2582 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
2583 'trunc' as rounding mode to make it equivalent to:
2588 You can override this by passing the desired rounding mode as parameter to
2591 $x = Math::BigFloat->new(2.5);
2592 $y = $x->as_number('odd'); # $y = 3
2598 =head1 Autocreating constants
2600 After C<use Math::BigFloat ':constant'> all the floating point constants
2601 in the given scope are converted to C<Math::BigFloat>. This conversion
2602 happens at compile time.
2606 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2608 prints the value of C<2E-100>. Note that without conversion of
2609 constants the expression 2E-100 will be calculated as normal floating point
2612 Please note that ':constant' does not affect integer constants, nor binary
2613 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2618 Math with the numbers is done (by default) by a module called
2619 Math::BigInt::Calc. This is equivalent to saying:
2621 use Math::BigFloat lib => 'Calc';
2623 You can change this by using:
2625 use Math::BigFloat lib => 'BitVect';
2627 The following would first try to find Math::BigInt::Foo, then
2628 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2630 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2632 Calc.pm uses as internal format an array of elements of some decimal base
2633 (usually 1e7, but this might be differen for some systems) with the least
2634 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2635 significant bit first. Other modules might use even different means of
2636 representing the numbers. See the respective module documentation for further
2639 Please note that Math::BigFloat does B<not> use the denoted library itself,
2640 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2643 use Math::BigInt lib => 'GMP';
2646 you can roll it all into one line:
2648 use Math::BigFloat lib => 'GMP';
2650 It is also possible to just require Math::BigFloat:
2652 require Math::BigFloat;
2654 This will load the neccessary things (like BigInt) when they are needed, and
2657 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2658 you ever wanted to know about loading a different library.
2660 =head2 Using Math::BigInt::Lite
2662 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2665 use Math::BigFloat with => 'Math::BigInt::Lite';
2667 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2668 can combine these if you want. For instance, you may want to use
2669 Math::BigInt objects in your main script, too.
2673 use Math::BigFloat with => 'Math::BigInt::Lite';
2675 Of course, you can combine this with the C<lib> parameter.
2678 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2680 There is no need for a "use Math::BigInt;" statement, even if you want to
2681 use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2682 always loads it. But if you add it, add it B<before>:
2686 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2688 Notice that the module with the last C<lib> will "win" and thus
2689 it's lib will be used if the lib is available:
2692 use Math::BigInt lib => 'Bar,Baz';
2693 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2695 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2696 words, Math::BigFloat will try to retain previously loaded libs when you
2697 don't specify it onem but if you specify one, it will try to load them.
2699 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2700 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2701 same as trying the latter load alone, except for the fact that one of Bar or
2702 Baz might be loaded needlessly in an intermidiate step (and thus hang around
2703 and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2704 will still be tried to be loaded, but this is not as time/memory consuming as
2705 actually loading one of them. Still, this type of usage is not recommended due
2708 The old way (loading the lib only in BigInt) still works though:
2711 use Math::BigInt lib => 'Bar,Baz';
2714 You can even load Math::BigInt afterwards:
2718 use Math::BigInt lib => 'Bar,Baz';
2720 But this has the same problems like #5, it will first load Calc
2721 (Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2722 Baz, depending on which of them works and is usable/loadable. Since this
2723 loads Calc unnecc., it is not recommended.
2725 Since it also possible to just require Math::BigFloat, this poses the question
2726 about what libary this will use:
2728 require Math::BigFloat;
2729 my $x = Math::BigFloat->new(123); $x += 123;
2731 It will use Calc. Please note that the call to import() is still done, but
2732 only when you use for the first time some Math::BigFloat math (it is triggered
2733 via any constructor, so the first time you create a Math::BigFloat, the load
2734 will happen in the background). This means:
2736 require Math::BigFloat;
2737 Math::BigFloat->import ( lib => 'Foo,Bar' );
2739 would be the same as:
2741 use Math::BigFloat lib => 'Foo, Bar';
2743 But don't try to be clever to insert some operations in between:
2745 require Math::BigFloat;
2746 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
2747 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
2748 $x = Math::BigFloat->bone()+4; # now use Pari
2750 While this works, it loads Calc needlessly. But maybe you just wanted that?
2752 B<Examples #3 is highly recommended> for daily usage.
2756 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
2762 =item stringify, bstr()
2764 Both stringify and bstr() now drop the leading '+'. The old code would return
2765 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2766 reasoning and details.
2770 The following will probably not do what you expect:
2772 print $c->bdiv(123.456),"\n";
2774 It prints both quotient and reminder since print works in list context. Also,
2775 bdiv() will modify $c, so be carefull. You probably want to use
2777 print $c / 123.456,"\n";
2778 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2782 =item Modifying and =
2786 $x = Math::BigFloat->new(5);
2789 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2790 a second reference to the B<same> object and stores it in $y. Thus anything
2791 that modifies $x will modify $y (except overloaded math operators), and vice
2792 versa. See L<Math::BigInt> for details and how to avoid that.
2796 C<bpow()> now modifies the first argument, unlike the old code which left
2797 it alone and only returned the result. This is to be consistent with
2798 C<badd()> etc. The first will modify $x, the second one won't:
2800 print bpow($x,$i),"\n"; # modify $x
2801 print $x->bpow($i),"\n"; # ditto
2802 print $x ** $i,"\n"; # leave $x alone
2808 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
2809 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
2811 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
2812 because they solve the autoupgrading/downgrading issue, at least partly.
2815 L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
2816 more documentation including a full version history, testcases, empty
2817 subclass files and benchmarks.
2821 This program is free software; you may redistribute it and/or modify it under
2822 the same terms as Perl itself.
2826 Mark Biggar, overloaded interface by Ilya Zakharevich.
2827 Completely rewritten by Tels http://bloodgate.com in 2001, 2002, and still