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
18 @ISA = qw(Exporter Math::BigInt);
21 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode/;
22 use vars qw/$upgrade $downgrade/;
23 # the following are internal and should never be accessed from the outside
24 use vars qw/$_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 use constant 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";
154 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
160 Carp::croak ("$wanted is not a number initialized to $class");
163 return $downgrade->bnan() if $downgrade;
165 $self->{_e} = $MBI->bzero();
166 $self->{_m} = $MBI->bzero();
167 $self->{sign} = $nan;
171 # make integer from mantissa by adjusting exp, then convert to bigint
172 # undef,undef to signal MBI that we don't need no bloody rounding
173 $self->{_e} = $MBI->new("$$es$$ev",undef,undef); # exponent
174 $self->{_m} = $MBI->new("$$miv$$mfv",undef,undef); # create mant.
175 # print $self->{_e}, " ", $self->{_m},"\n";
176 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
177 $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;
178 $self->{sign} = $$mis;
180 # if downgrade, inf, NaN or integers go down
182 if ($downgrade && $self->{_e}->{sign} eq '+')
184 #print "downgrading $$miv$$mfv"."E$$es$$ev";
185 if ($self->{_e}->is_zero())
187 $self->{_m}->{sign} = $$mis; # negative if wanted
188 return $downgrade->new($self->{_m});
190 return $downgrade->new($self->bsstr());
192 #print "mbf new $self->{sign} $self->{_m} e $self->{_e} ",ref($self),"\n";
193 $self->bnorm()->round(@r); # first normalize, then round
198 # used by parent class bone() to initialize number to NaN
204 my $class = ref($self);
205 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
208 $IMPORT=1; # call our import only once
209 $self->{_m} = $MBI->bzero();
210 $self->{_e} = $MBI->bzero();
215 # used by parent class bone() to initialize number to +-inf
221 my $class = ref($self);
222 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
225 $IMPORT=1; # call our import only once
226 $self->{_m} = $MBI->bzero();
227 $self->{_e} = $MBI->bzero();
232 # used by parent class bone() to initialize number to 1
234 $IMPORT=1; # call our import only once
235 $self->{_m} = $MBI->bone();
236 $self->{_e} = $MBI->bzero();
241 # used by parent class bone() to initialize number to 0
243 $IMPORT=1; # call our import only once
244 $self->{_m} = $MBI->bzero();
245 $self->{_e} = $MBI->bone();
250 my ($self,$class) = @_;
251 return if $class =~ /^Math::BigInt/; # we aren't one of these
252 UNIVERSAL::isa($self,$class);
257 # return (later set?) configuration data as hash ref
258 my $class = shift || 'Math::BigFloat';
260 my $cfg = $class->SUPER::config(@_);
262 # now we need only to override the ones that are different from our parent
263 $cfg->{class} = $class;
268 ##############################################################################
269 # string conversation
273 # (ref to BFLOAT or num_str ) return num_str
274 # Convert number from internal format to (non-scientific) string format.
275 # internal format is always normalized (no leading zeros, "-0" => "+0")
276 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
278 if ($x->{sign} !~ /^[+-]$/)
280 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
284 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
287 my $not_zero = !($x->{sign} eq '+' && $x->{_m}->is_zero());
290 $es = $x->{_m}->bstr();
291 $len = CORE::length($es);
292 my $e = $x->{_e}->numify();
296 # if _e is bigger than a scalar, the following will blow your memory
299 #print "style: 0.xxxx\n";
300 my $r = abs($e) - $len;
301 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
305 #print "insert '.' at $e in '$es'\n";
306 substr($es,$e,0) = '.'; $cad = $x->{_e};
312 $es .= '0' x $e; $len += $e; $cad = 0;
315 $es = '-'.$es if $x->{sign} eq '-';
316 # if set accuracy or precision, pad with zeros on the right side
317 if ((defined $x->{_a}) && ($not_zero))
319 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
320 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
321 $zeros = $x->{_a} - $len if $cad != $len;
322 $es .= $dot.'0' x $zeros if $zeros > 0;
324 elsif ((($x->{_p} || 0) < 0))
326 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
327 my $zeros = -$x->{_p} + $cad;
328 $es .= $dot.'0' x $zeros if $zeros > 0;
335 # (ref to BFLOAT or num_str ) return num_str
336 # Convert number from internal format to scientific string format.
337 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
338 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
340 if ($x->{sign} !~ /^[+-]$/)
342 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
345 # do $esign, because we need '1e+1', since $x->{_e}->bstr() misses the +
346 my $esign = $x->{_e}->{sign}; $esign = '' if $esign eq '-';
347 my $sep = 'e'.$esign;
348 my $sign = $x->{sign}; $sign = '' if $sign eq '+';
349 $sign . $x->{_m}->bstr() . $sep . $x->{_e}->bstr();
354 # Make a number from a BigFloat object
355 # simple return string and let Perl's atoi()/atof() handle the rest
356 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
360 ##############################################################################
361 # public stuff (usually prefixed with "b")
364 # todo: this must be overwritten and return NaN for non-integer values
365 # band(), bior(), bxor(), too
368 # $class->SUPER::bnot($class,@_);
373 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
374 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
377 my ($self,$x,$y) = (ref($_[0]),@_);
378 # objectify is costly, so avoid it
379 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
381 ($self,$x,$y) = objectify(2,@_);
384 return $upgrade->bcmp($x,$y) if defined $upgrade &&
385 ((!$x->isa($self)) || (!$y->isa($self)));
387 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
389 # handle +-inf and NaN
390 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
391 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
392 return +1 if $x->{sign} eq '+inf';
393 return -1 if $x->{sign} eq '-inf';
394 return -1 if $y->{sign} eq '+inf';
398 # check sign for speed first
399 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
400 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
403 my $xz = $x->is_zero();
404 my $yz = $y->is_zero();
405 return 0 if $xz && $yz; # 0 <=> 0
406 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
407 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
409 # adjust so that exponents are equal
410 my $lxm = $x->{_m}->length();
411 my $lym = $y->{_m}->length();
412 # the numify somewhat limits our length, but makes it much faster
413 my $lx = $lxm + $x->{_e}->numify();
414 my $ly = $lym + $y->{_e}->numify();
415 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
416 return $l <=> 0 if $l != 0;
418 # lengths (corrected by exponent) are equal
419 # so make mantissa equal length by padding with zero (shift left)
420 my $diff = $lxm - $lym;
421 my $xm = $x->{_m}; # not yet copy it
425 $ym = $y->{_m}->copy()->blsft($diff,10);
429 $xm = $x->{_m}->copy()->blsft(-$diff,10);
431 my $rc = $xm->bacmp($ym);
432 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
438 # Compares 2 values, ignoring their signs.
439 # Returns one of undef, <0, =0, >0. (suitable for sort)
440 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
443 my ($self,$x,$y) = (ref($_[0]),@_);
444 # objectify is costly, so avoid it
445 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
447 ($self,$x,$y) = objectify(2,@_);
450 return $upgrade->bacmp($x,$y) if defined $upgrade &&
451 ((!$x->isa($self)) || (!$y->isa($self)));
453 # handle +-inf and NaN's
454 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
456 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
457 return 0 if ($x->is_inf() && $y->is_inf());
458 return 1 if ($x->is_inf() && !$y->is_inf());
463 my $xz = $x->is_zero();
464 my $yz = $y->is_zero();
465 return 0 if $xz && $yz; # 0 <=> 0
466 return -1 if $xz && !$yz; # 0 <=> +y
467 return 1 if $yz && !$xz; # +x <=> 0
469 # adjust so that exponents are equal
470 my $lxm = $x->{_m}->length();
471 my $lym = $y->{_m}->length();
472 # the numify somewhat limits our length, but makes it much faster
473 my $lx = $lxm + $x->{_e}->numify();
474 my $ly = $lym + $y->{_e}->numify();
476 return $l <=> 0 if $l != 0;
478 # lengths (corrected by exponent) are equal
479 # so make mantissa equal-length by padding with zero (shift left)
480 my $diff = $lxm - $lym;
481 my $xm = $x->{_m}; # not yet copy it
485 $ym = $y->{_m}->copy()->blsft($diff,10);
489 $xm = $x->{_m}->copy()->blsft(-$diff,10);
491 $xm->bacmp($ym) <=> 0;
496 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
497 # return result as BFLOAT
500 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
501 # objectify is costly, so avoid it
502 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
504 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
507 # inf and NaN handling
508 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
511 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
513 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
515 # +inf++inf or -inf+-inf => same, rest is NaN
516 return $x if $x->{sign} eq $y->{sign};
519 # +-inf + something => +inf; something +-inf => +-inf
520 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
524 return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
525 ((!$x->isa($self)) || (!$y->isa($self)));
527 # speed: no add for 0+y or x+0
528 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
529 if ($x->is_zero()) # 0+y
531 # make copy, clobbering up x (modify in place!)
532 $x->{_e} = $y->{_e}->copy();
533 $x->{_m} = $y->{_m}->copy();
534 $x->{sign} = $y->{sign} || $nan;
535 return $x->round($a,$p,$r,$y);
538 # take lower of the two e's and adapt m1 to it to match m2
540 $e = $MBI->bzero() if !defined $e; # if no BFLOAT ?
541 $e = $e->copy(); # make copy (didn't do it yet)
542 $e->bsub($x->{_e}); # Ye - Xe
543 my $add = $y->{_m}->copy();
544 if ($e->{sign} eq '-') # < 0
546 $x->{_e} += $e; # need the sign of e
547 $x->{_m}->blsft($e->babs(),10); # destroys copy of _e
549 elsif (!$e->is_zero()) # > 0
553 # else: both e are the same, so just leave them
554 $x->{_m}->{sign} = $x->{sign}; # fiddle with signs
555 $add->{sign} = $y->{sign};
556 $x->{_m} += $add; # finally do add/sub
557 $x->{sign} = $x->{_m}->{sign}; # re-adjust signs
558 $x->{_m}->{sign} = '+'; # mantissa always positiv
559 # delete trailing zeros, then round
560 $x->bnorm()->round($a,$p,$r,$y);
565 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
566 # subtract second arg from first, modify first
569 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
570 # objectify is costly, so avoid it
571 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
573 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
577 if ($y->is_zero()) # still round for not adding zero
579 return $x->round($a,$p,$r);
583 $y->{sign} =~ tr/+-/-+/; # does nothing for NaN
584 $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
585 $y->{sign} =~ tr/+-/-+/; # refix $y (does nothing for NaN)
586 $x; # already rounded by badd()
591 # increment arg by one
592 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
594 if ($x->{_e}->sign() eq '-')
596 return $x->badd($self->bone(),$a,$p,$r); # digits after dot
599 if (!$x->{_e}->is_zero())
601 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
605 if ($x->{sign} eq '+')
608 return $x->bnorm()->bround($a,$p,$r);
610 elsif ($x->{sign} eq '-')
613 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
614 return $x->bnorm()->bround($a,$p,$r);
616 # inf, nan handling etc
617 $x->badd($self->bone(),$a,$p,$r); # does round
622 # decrement arg by one
623 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
625 if ($x->{_e}->sign() eq '-')
627 return $x->badd($self->bone('-'),$a,$p,$r); # digits after dot
630 if (!$x->{_e}->is_zero())
632 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
636 my $zero = $x->is_zero();
638 if (($x->{sign} eq '-') || $zero)
641 $x->{sign} = '-' if $zero; # 0 => 1 => -1
642 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
643 return $x->bnorm()->round($a,$p,$r);
646 elsif ($x->{sign} eq '+')
649 return $x->bnorm()->round($a,$p,$r);
651 # inf, nan handling etc
652 $x->badd($self->bone('-'),$a,$p,$r); # does round
659 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
661 # $base > 0, $base != 1; if $base == undef default to $base == e
664 # we need to limit the accuracy to protect against overflow
667 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
669 # also takes care of the "error in _find_round_parameters?" case
670 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
672 # no rounding at all, so must use fallback
673 if (scalar @params == 0)
675 # simulate old behaviour
676 $params[0] = $self->div_scale(); # and round to it as accuracy
677 $params[1] = undef; # P = undef
678 $scale = $params[0]+4; # at least four more for proper round
679 $params[2] = $r; # round mode by caller or undef
680 $fallback = 1; # to clear a/p afterwards
684 # the 4 below is empirical, and there might be cases where it is not
686 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
689 return $x->bzero(@params) if $x->is_one();
690 # base not defined => base == Euler's constant e
693 # make object, since we don't feed it through objectify() to still get the
694 # case of $base == undef
695 $base = $self->new($base) unless ref($base);
696 # $base > 0; $base != 1
697 return $x->bnan() if $base->is_zero() || $base->is_one() ||
698 $base->{sign} ne '+';
699 # if $x == $base, we know the result must be 1.0
700 return $x->bone('+',@params) if $x->bcmp($base) == 0;
703 # when user set globals, they would interfere with our calculation, so
704 # disable them and later re-enable them
706 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
707 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
708 # we also need to disable any set A or P on $x (_find_round_parameters took
709 # them already into account), since these would interfere, too
710 delete $x->{_a}; delete $x->{_p};
711 # need to disable $upgrade in BigInt, to avoid deep recursion
712 local $Math::BigInt::upgrade = undef;
713 local $Math::BigFloat::downgrade = undef;
715 # upgrade $x if $x is not a BigFloat (handle BigInt input)
716 if (!$x->isa('Math::BigFloat'))
718 $x = Math::BigFloat->new($x);
721 # first calculate the log to base e (using reduction by 10 (and probably 2))
722 $self->_log_10($x,$scale);
724 # and if a different base was requested, convert it
727 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
728 # not ln, but some other base (don't modify $base)
729 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
732 # shortcut to not run through _find_round_parameters again
733 if (defined $params[0])
735 $x->bround($params[0],$params[2]); # then round accordingly
739 $x->bfround($params[1],$params[2]); # then round accordingly
743 # clear a/p after round, since user did not request it
744 $x->{_a} = undef; $x->{_p} = undef;
747 $$abr = $ab; $$pbr = $pb;
754 # internal log function to calculate ln() based on Taylor series.
755 # Modifies $x in place.
756 my ($self,$x,$scale) = @_;
758 # in case of $x == 1, result is 0
759 return $x->bzero() if $x->is_one();
761 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
765 # Taylor: | u 1 u^3 1 u^5 |
766 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
767 # |_ v 3 v^3 5 v^5 _|
769 # This takes much more steps to calculate the result and is thus not used
772 # Taylor: | u 1 u^2 1 u^3 |
773 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
774 # |_ x 2 x^2 3 x^3 _|
776 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
778 $v = $x->copy(); $v->binc(); # v = x+1
779 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
780 $x->bdiv($v,$scale); # first term: u/v
783 $u *= $u; $v *= $v; # u^2, v^2
784 $below->bmul($v); # u^3, v^3
786 $factor = $self->new(3); $f = $self->new(2);
788 my $steps = 0 if DEBUG;
789 $limit = $self->new("1E-". ($scale-1));
792 # we calculate the next term, and add it to the last
793 # when the next term is below our limit, it won't affect the outcome
794 # anymore, so we stop
796 # calculating the next term simple from over/below will result in quite
797 # a time hog if the input has many digits, since over and below will
798 # accumulate more and more digits, and the result will also have many
799 # digits, but in the end it is rounded to $scale digits anyway. So if we
800 # round $over and $below first, we save a lot of time for the division
801 # (not with log(1.2345), but try log (123**123) to see what I mean. This
802 # can introduce a rounding error if the division result would be f.i.
803 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
804 # if we truncated $over and $below we might get 0.12345. Does this matter
805 # for the end result? So we give $over and $below 4 more digits to be
806 # on the safe side (unscientific error handling as usual... :+D
808 $next = $over->copy->bround($scale+4)->bdiv(
809 $below->copy->bmul($factor)->bround($scale+4),
813 ## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
815 last if $next->bacmp($limit) <= 0;
817 delete $next->{_a}; delete $next->{_p};
819 #print "step $x\n ($next - $limit = ",$next - $limit,")\n";
820 # calculate things for the next term
821 $over *= $u; $below *= $v; $factor->badd($f);
824 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
827 $x->bmul($f); # $x *= 2
828 print "took $steps steps\n" if DEBUG;
833 # Internal log function based on reducing input to the range of 0.1 .. 9.99
834 # and then "correcting" the result to the proper one. Modifies $x in place.
835 my ($self,$x,$scale) = @_;
837 # taking blog() from numbers greater than 10 takes a *very long* time, so we
838 # break the computation down into parts based on the observation that:
839 # blog(x*y) = blog(x) + blog(y)
840 # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
841 # the faster it get's, especially because 2*$x takes about 10 times as long,
842 # so by dividing $x by 10 we make it at least factor 100 faster...)
844 # The same observation is valid for numbers smaller than 0.1 (e.g. computing
845 # log(1) is fastest, and the farther away we get from 1, the longer it takes)
846 # so we also 'break' this down by multiplying $x with 10 and subtract the
847 # log(10) afterwards to get the correct result.
849 # calculate nr of digits before dot
850 my $dbd = $x->{_m}->length() + $x->{_e}->numify();
852 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
855 my $calc = 1; # do some calculation?
857 # disable the shortcut for 10, since we need log(10) and this would recurse
859 if ($x->{_e}->is_one() && $x->{_m}->is_one())
861 $dbd = 0; # disable shortcut
862 # we can use the cached value in these cases
863 if ($scale <= $LOG_10_A)
865 $x->bzero(); $x->badd($LOG_10);
866 $calc = 0; # no need to calc, but round
871 # disable the shortcut for 2, since we maybe have it cached
872 if ($x->{_e}->is_zero() && $x->{_m}->bcmp(2) == 0)
874 $dbd = 0; # disable shortcut
875 # we can use the cached value in these cases
876 if ($scale <= $LOG_2_A)
878 $x->bzero(); $x->badd($LOG_2);
879 $calc = 0; # no need to calc, but round
884 # if $x = 0.1, we know the result must be 0-log(10)
885 if ($calc != 0 && $x->{_e}->is_one('-') && $x->{_m}->is_one())
887 $dbd = 0; # disable shortcut
888 # we can use the cached value in these cases
889 if ($scale <= $LOG_10_A)
891 $x->bzero(); $x->bsub($LOG_10);
892 $calc = 0; # no need to calc, but round
896 return if $calc == 0; # already have the result
898 # default: these correction factors are undef and thus not used
899 my $l_10; # value of ln(10) to A of $scale
900 my $l_2; # value of ln(2) to A of $scale
902 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
903 # so don't do this shortcut for 1 or 0
904 if (($dbd > 1) || ($dbd < 0))
906 # convert our cached value to an object if not already (avoid doing this
907 # at import() time, since not everybody needs this)
908 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
910 #print "x = $x, dbd = $dbd, calc = $calc\n";
911 # got more than one digit before the dot, or more than one zero after the
913 # log(123) == log(1.23) + log(10) * 2
914 # log(0.0123) == log(1.23) - log(10) * 2
916 if ($scale <= $LOG_10_A)
919 #print "using cached value for l_10\n";
920 $l_10 = $LOG_10->copy(); # copy for mul
924 # else: slower, compute it (but don't cache it, because it could be big)
925 # also disable downgrade for this code path
926 local $Math::BigFloat::downgrade = undef;
927 #print "l_10 = $l_10 (self = $self',
928 # ", ref(l_10) = ",ref($l_10)," scale $scale)\n";
929 #print "calculating value for l_10, scale $scale\n";
930 $l_10 = $self->new(10)->blog(undef,$scale); # scale+4, actually
932 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
934 $dbd = $self->new($dbd);
936 $l_10->bmul($dbd); # log(10) * (digits_before_dot-1)
937 #print "l_10 = $l_10\n";
939 $x->{_e}->bsub($dbd); # 123 => 1.23
941 #print "calculating log($x) with scale=$scale\n";
945 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
947 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
948 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
950 my $half = $self->new('0.5');
951 my $twos = 0; # default: none (0 times)
952 my $two = $self->new(2);
953 while ($x->bacmp($half) <= 0)
955 $twos--; $x->bmul($two);
957 while ($x->bacmp($two) >= 0)
959 $twos++; $x->bdiv($two,$scale+4); # keep all digits
962 # $twos > 0 => did mul 2, < 0 => did div 2 (never both)
963 # calculate correction factor based on ln(2)
966 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
967 if ($scale <= $LOG_2_A)
970 #print "using cached value for l_10\n";
971 $l_2 = $LOG_2->copy(); # copy for mul
975 # else: slower, compute it (but don't cache it, because it could be big)
976 # also disable downgrade for this code path
977 local $Math::BigFloat::downgrade = undef;
978 #print "calculating value for l_2, scale $scale\n";
979 $l_2 = $two->blog(undef,$scale); # scale+4, actually
981 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
984 $self->_log($x,$scale); # need to do the "normal" way
985 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
986 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
987 # all done, $x contains now the result
992 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
993 # does not modify arguments, but returns new object
994 # Lowest Common Multiplicator
996 my ($self,@arg) = objectify(0,@_);
997 my $x = $self->new(shift @arg);
998 while (@arg) { $x = _lcm($x,shift @arg); }
1004 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
1005 # does not modify arguments, but returns new object
1006 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
1008 my ($self,@arg) = objectify(0,@_);
1009 my $x = $self->new(shift @arg);
1010 while (@arg) { $x = _gcd($x,shift @arg); }
1014 ###############################################################################
1015 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1019 # internal, return true if BigInt arg is zero or one, saving the
1020 # two calls to is_zero() and is_one()
1023 $x->{sign} eq '+' && ($x->is_zero() || $x->is_one());
1028 # return true if arg (BFLOAT or num_str) is an integer
1029 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1031 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1032 $x->{_e}->{sign} eq '+'; # 1e-1 => no integer
1038 # return true if arg (BFLOAT or num_str) is zero
1039 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1041 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
1047 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1048 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1050 $sign = '+' if !defined $sign || $sign ne '-';
1052 if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
1058 # return true if arg (BFLOAT or num_str) is odd or false if even
1059 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1061 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1062 ($x->{_e}->is_zero() && $x->{_m}->is_odd());
1068 # return true if arg (BINT or num_str) is even or false if odd
1069 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1071 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
1072 return 1 if ($x->{_e}->{sign} eq '+' # 123.45 is never
1073 && $x->{_m}->is_even()); # but 1200 is
1079 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
1080 # (BINT or num_str, BINT or num_str) return BINT
1083 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1084 # objectify is costly, so avoid it
1085 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1087 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1090 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1093 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1095 return $x->bnan() if $x->is_zero() || $y->is_zero();
1096 # result will always be +-inf:
1097 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1098 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1099 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1100 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1101 return $x->binf('-');
1104 return $x->bzero() if $x->is_zero() || $y->is_zero();
1106 return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
1107 ((!$x->isa($self)) || (!$y->isa($self)));
1109 # aEb * cEd = (a*c)E(b+d)
1110 $x->{_m}->bmul($y->{_m});
1111 $x->{_e}->badd($y->{_e});
1113 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1114 return $x->bnorm()->round($a,$p,$r,$y);
1119 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1120 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1123 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1124 # objectify is costly, so avoid it
1125 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1127 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1130 return $self->_div_inf($x,$y)
1131 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1133 # x== 0 # also: or y == 1 or y == -1
1134 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1137 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1139 # we need to limit the accuracy to protect against overflow
1141 my (@params,$scale);
1142 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1144 return $x if $x->is_nan(); # error in _find_round_parameters?
1146 # no rounding at all, so must use fallback
1147 if (scalar @params == 0)
1149 # simulate old behaviour
1150 $params[0] = $self->div_scale(); # and round to it as accuracy
1151 $scale = $params[0]+4; # at least four more for proper round
1152 $params[2] = $r; # round mode by caller or undef
1153 $fallback = 1; # to clear a/p afterwards
1157 # the 4 below is empirical, and there might be cases where it is not
1159 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1161 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
1162 $scale = $lx if $lx > $scale;
1163 $scale = $ly if $ly > $scale;
1164 my $diff = $ly - $lx;
1165 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
1167 # make copy of $x in case of list context for later reminder calculation
1169 if (wantarray && !$y->is_one())
1174 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1176 # check for / +-1 ( +/- 1E0)
1179 # promote BigInts and it's subclasses (except when already a BigFloat)
1180 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1182 # need to disable $upgrade in BigInt, to avoid deep recursion
1183 local $Math::BigInt::upgrade = undef; # should be parent class vs MBI
1185 # calculate the result to $scale digits and then round it
1186 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1187 $x->{_m}->blsft($scale,10);
1188 $x->{_m}->bdiv( $y->{_m} ); # a/c
1189 $x->{_e}->bsub( $y->{_e} ); # b-d
1190 $x->{_e}->bsub($scale); # correct for 10**scale
1191 $x->bnorm(); # remove trailing 0's
1194 # shortcut to not run through _find_round_parameters again
1195 if (defined $params[0])
1197 $x->{_a} = undef; # clear before round
1198 $x->bround($params[0],$params[2]); # then round accordingly
1202 $x->{_p} = undef; # clear before round
1203 $x->bfround($params[1],$params[2]); # then round accordingly
1207 # clear a/p after round, since user did not request it
1208 $x->{_a} = undef; $x->{_p} = undef;
1215 $rem->bmod($y,@params); # copy already done
1219 $rem = $self->bzero();
1223 # clear a/p after round, since user did not request it
1224 $rem->{_a} = undef; $rem->{_p} = undef;
1233 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1236 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1237 # objectify is costly, so avoid it
1238 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1240 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1243 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1245 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1246 $x->{sign} = $re->{sign};
1247 $x->{_e} = $re->{_e};
1248 $x->{_m} = $re->{_m};
1249 return $x->round($a,$p,$r,$y);
1251 return $x->bnan() if $x->is_zero() && $y->is_zero();
1252 return $x if $y->is_zero();
1253 return $x->bnan() if $x->is_nan() || $y->is_nan();
1254 return $x->bzero() if $y->is_one() || $x->is_zero();
1256 # inf handling is missing here
1258 my $cmp = $x->bacmp($y); # equal or $x < $y?
1259 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1261 # only $y of the operands negative?
1262 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1264 $x->{sign} = $y->{sign}; # calc sign first
1265 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1267 my $ym = $y->{_m}->copy();
1270 $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
1272 # if $y has digits after dot
1273 my $shifty = 0; # correct _e of $x by this
1274 if ($y->{_e}->{sign} eq '-') # has digits after dot
1276 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1277 $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
1278 $x->blsft($shifty,10); # 123 => 1230, $y->{_m} is already 25
1280 # $ym is now mantissa of $y based on exponent 0
1282 my $shiftx = 0; # correct _e of $x by this
1283 if ($x->{_e}->{sign} eq '-') # has digits after dot
1285 # 123.4 % 20 => 1234 % 200
1286 $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
1287 $ym->blsft($shiftx,10);
1289 # 123e1 % 20 => 1230 % 20
1290 if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
1292 $x->{_m}->blsft($x->{_e},10);
1294 $x->{_e} = $MBI->bzero() unless $x->{_e}->is_zero();
1296 $x->{_e}->bsub($shiftx) if $shiftx != 0;
1297 $x->{_e}->bsub($shifty) if $shifty != 0;
1299 # now mantissas are equalized, exponent of $x is adjusted, so calc result
1301 $x->{_m}->bmod($ym);
1303 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
1306 if ($neg != 0) # one of them negative => correct in place
1309 $x->{_m} = $r->{_m};
1310 $x->{_e} = $r->{_e};
1311 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
1315 $x->round($a,$p,$r,$y); # round and return
1320 # calculate $y'th root of $x
1321 my ($self,$x,$y,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(2,@_);
1323 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1324 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1325 $y->{sign} !~ /^\+$/;
1327 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1329 # we need to limit the accuracy to protect against overflow
1331 my (@params,$scale);
1332 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1334 return $x if $x->is_nan(); # error in _find_round_parameters?
1336 # no rounding at all, so must use fallback
1337 if (scalar @params == 0)
1339 # simulate old behaviour
1340 $params[0] = $self->div_scale(); # and round to it as accuracy
1341 $scale = $params[0]+4; # at least four more for proper round
1342 $params[2] = $r; # round mode by caller or undef
1343 $fallback = 1; # to clear a/p afterwards
1347 # the 4 below is empirical, and there might be cases where it is not
1349 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1352 # when user set globals, they would interfere with our calculation, so
1353 # disable them and later re-enable them
1355 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1356 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1357 # we also need to disable any set A or P on $x (_find_round_parameters took
1358 # them already into account), since these would interfere, too
1359 delete $x->{_a}; delete $x->{_p};
1360 # need to disable $upgrade in BigInt, to avoid deep recursion
1361 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1363 # remember sign and make $x positive, since -4 ** (1/2) => -2
1364 my $sign = 0; $sign = 1 if $x->is_negative(); $x->babs();
1366 if ($y->bcmp(2) == 0) # normal square root
1368 $x->bsqrt($scale+4);
1370 elsif ($y->is_one('-'))
1373 my $u = $self->bone()->bdiv($x,$scale);
1374 # copy private parts over
1375 $x->{_m} = $u->{_m};
1376 $x->{_e} = $u->{_e};
1380 my $u = $self->bone()->bdiv($y,$scale+4);
1381 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
1382 $x->bpow($u,$scale+4); # el cheapo
1384 $x->bneg() if $sign == 1;
1386 # shortcut to not run through _find_round_parameters again
1387 if (defined $params[0])
1389 $x->bround($params[0],$params[2]); # then round accordingly
1393 $x->bfround($params[1],$params[2]); # then round accordingly
1397 # clear a/p after round, since user did not request it
1398 $x->{_a} = undef; $x->{_p} = undef;
1401 $$abr = $ab; $$pbr = $pb;
1407 # calculate square root
1408 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1410 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
1411 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
1412 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
1414 # we need to limit the accuracy to protect against overflow
1416 my (@params,$scale);
1417 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1419 return $x if $x->is_nan(); # error in _find_round_parameters?
1421 # no rounding at all, so must use fallback
1422 if (scalar @params == 0)
1424 # simulate old behaviour
1425 $params[0] = $self->div_scale(); # and round to it as accuracy
1426 $scale = $params[0]+4; # at least four more for proper round
1427 $params[2] = $r; # round mode by caller or undef
1428 $fallback = 1; # to clear a/p afterwards
1432 # the 4 below is empirical, and there might be cases where it is not
1434 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1437 # when user set globals, they would interfere with our calculation, so
1438 # disable them and later re-enable them
1440 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1441 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1442 # we also need to disable any set A or P on $x (_find_round_parameters took
1443 # them already into account), since these would interfere, too
1444 delete $x->{_a}; delete $x->{_p};
1445 # need to disable $upgrade in BigInt, to avoid deep recursion
1446 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1448 my $xas = $x->as_number();
1449 my $gs = $xas->copy()->bsqrt(); # some guess
1451 if (($x->{_e}->{sign} ne '-') # guess can't be accurate if there are
1452 # digits after the dot
1453 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
1456 $x->{_m} = $gs; $x->{_e} = $MBI->bzero(); $x->bnorm();
1457 # shortcut to not run through _find_round_parameters again
1458 if (defined $params[0])
1460 $x->bround($params[0],$params[2]); # then round accordingly
1464 $x->bfround($params[1],$params[2]); # then round accordingly
1468 # clear a/p after round, since user did not request it
1469 $x->{_a} = undef; $x->{_p} = undef;
1471 # re-enable A and P, upgrade is taken care of by "local"
1472 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1476 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1477 # of the result by multipyling the input by 100 and then divide the integer
1478 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
1479 # this will transform 123.456 (in $x) into 123456 (in $y1)
1480 my $y1 = $x->{_m}->copy();
1481 # We now make sure that $y1 has the same odd or even number of digits than
1482 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1483 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1484 # steps of 10. The length of $x does not count, since an even or odd number
1485 # of digits before the dot is not changed by adding an even number of digits
1486 # after the dot (the result is still odd or even digits long).
1487 my $length = $y1->length();
1488 $y1->bmul(10) if $x->{_e}->is_odd();
1489 # now calculate how many digits the result of sqrt(y1) would have
1490 my $digits = int($length / 2);
1491 # but we need at least $scale digits, so calculate how many are missing
1492 my $shift = $scale - $digits;
1493 # that should never happen (we take care of integer guesses above)
1494 # $shift = 0 if $shift < 0;
1495 # multiply in steps of 100, by shifting left two times the "missing" digits
1496 $y1->blsft($shift*2,10);
1497 # now take the square root and truncate to integer
1499 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1500 # result, which is than later rounded to the desired scale.
1502 # calculate how many zeros $x had after the '.' (or before it, depending
1503 # on sign of $dat, the result should have half as many:
1504 my $dat = $length + $x->{_e}->numify();
1508 # no zeros after the dot (e.g. 1.23, 0.49 etc)
1509 # preserve half as many digits before the dot than the input had
1510 # (but round this "up")
1511 $dat = int(($dat+1)/2);
1515 $dat = int(($dat)/2);
1517 $x->{_e}= $MBI->new( $dat - $y1->length() );
1521 # shortcut to not run through _find_round_parameters again
1522 if (defined $params[0])
1524 $x->bround($params[0],$params[2]); # then round accordingly
1528 $x->bfround($params[1],$params[2]); # then round accordingly
1532 # clear a/p after round, since user did not request it
1533 $x->{_a} = undef; $x->{_p} = undef;
1536 $$abr = $ab; $$pbr = $pb;
1542 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1543 # compute factorial number, modifies first argument
1544 my ($self,$x,@r) = objectify(1,@_);
1547 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1548 ($x->{_e}->{sign} ne '+')); # digits after dot?
1550 # use BigInt's bfac() for faster calc
1551 if (! _is_zero_or_one($x->{_e}))
1553 $x->{_m}->blsft($x->{_e},10); # unnorm
1554 $x->{_e}->bzero(); # norm again
1556 $x->{_m}->blsft($x->{_e},10); # un-norm m
1557 $x->{_e}->bzero(); # norm again
1558 $x->{_m}->bfac(); # calculate factorial
1559 $x->bnorm()->round(@r); # norm again and round result
1564 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1565 my ($x,$y,$a,$p,$r) = @_;
1568 # if $y == 0.5, it is sqrt($x)
1569 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1572 # a ** x == e ** (x * ln a)
1576 # Taylor: | u u^2 u^3 |
1577 # x ** y = 1 + | --- + --- + ----- + ... |
1580 # we need to limit the accuracy to protect against overflow
1582 my ($scale,@params);
1583 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1585 return $x if $x->is_nan(); # error in _find_round_parameters?
1587 # no rounding at all, so must use fallback
1588 if (scalar @params == 0)
1590 # simulate old behaviour
1591 $params[0] = $self->div_scale(); # and round to it as accuracy
1592 $params[1] = undef; # disable P
1593 $scale = $params[0]+4; # at least four more for proper round
1594 $params[2] = $r; # round mode by caller or undef
1595 $fallback = 1; # to clear a/p afterwards
1599 # the 4 below is empirical, and there might be cases where it is not
1601 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1604 # when user set globals, they would interfere with our calculation, so
1605 # disable them and later re-enable them
1607 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1608 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1609 # we also need to disable any set A or P on $x (_find_round_parameters took
1610 # them already into account), since these would interfere, too
1611 delete $x->{_a}; delete $x->{_p};
1612 # need to disable $upgrade in BigInt, to avoid deep recursion
1613 local $Math::BigInt::upgrade = undef;
1615 my ($limit,$v,$u,$below,$factor,$next,$over);
1617 $u = $x->copy()->blog(undef,$scale)->bmul($y);
1618 $v = $self->bone(); # 1
1619 $factor = $self->new(2); # 2
1620 $x->bone(); # first term: 1
1622 $below = $v->copy();
1625 $limit = $self->new("1E-". ($scale-1));
1629 # we calculate the next term, and add it to the last
1630 # when the next term is below our limit, it won't affect the outcome
1631 # anymore, so we stop
1632 $next = $over->copy()->bdiv($below,$scale);
1633 last if $next->bacmp($limit) <= 0;
1635 # calculate things for the next term
1636 $over *= $u; $below *= $factor; $factor->binc();
1640 # shortcut to not run through _find_round_parameters again
1641 if (defined $params[0])
1643 $x->bround($params[0],$params[2]); # then round accordingly
1647 $x->bfround($params[1],$params[2]); # then round accordingly
1651 # clear a/p after round, since user did not request it
1652 $x->{_a} = undef; $x->{_p} = undef;
1655 $$abr = $ab; $$pbr = $pb;
1661 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1662 # compute power of two numbers, second arg is used as integer
1663 # modifies first argument
1666 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1667 # objectify is costly, so avoid it
1668 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1670 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1673 return $x if $x->{sign} =~ /^[+-]inf$/;
1674 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1675 return $x->bone() if $y->is_zero();
1676 return $x if $x->is_one() || $y->is_one();
1678 return $x->_pow($y,$a,$p,$r) if !$y->is_int(); # non-integer power
1680 my $y1 = $y->as_number(); # make bigint
1682 if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1684 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
1685 return $y1->is_odd() ? $x : $x->babs(1);
1689 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1690 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1694 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1696 $x->{_m}->bpow($y1);
1697 $x->{_e}->bmul($y1);
1698 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1700 if ($y->{sign} eq '-')
1702 # modify $x in place!
1703 my $z = $x->copy(); $x->bzero()->binc();
1704 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
1706 $x->round($a,$p,$r,$y);
1709 ###############################################################################
1710 # rounding functions
1714 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1715 # $n == 0 means round to integer
1716 # expects and returns normalized numbers!
1717 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1719 return $x if $x->modify('bfround');
1721 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1722 return $x if !defined $scale; # no-op
1724 # never round a 0, +-inf, NaN
1727 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1730 return $x if $x->{sign} !~ /^[+-]$/;
1732 # don't round if x already has lower precision
1733 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1735 $x->{_p} = $scale; # remember round in any case
1736 $x->{_a} = undef; # and clear A
1739 # round right from the '.'
1741 return $x if $x->{_e}->{sign} eq '+'; # e >= 0 => nothing to round
1743 $scale = -$scale; # positive for simplicity
1744 my $len = $x->{_m}->length(); # length of mantissa
1746 # the following poses a restriction on _e, but if _e is bigger than a
1747 # scalar, you got other problems (memory etc) anyway
1748 my $dad = -($x->{_e}->numify()); # digits after dot
1749 my $zad = 0; # zeros after dot
1750 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
1752 #print "scale $scale dad $dad zad $zad len $len\n";
1753 # number bsstr len zad dad
1754 # 0.123 123e-3 3 0 3
1755 # 0.0123 123e-4 3 1 4
1758 # 1.2345 12345e-4 5 0 4
1760 # do not round after/right of the $dad
1761 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1763 # round to zero if rounding inside the $zad, but not for last zero like:
1764 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1765 return $x->bzero() if $scale < $zad;
1766 if ($scale == $zad) # for 0.006, scale -3 and trunc
1772 # adjust round-point to be inside mantissa
1775 $scale = $scale-$zad;
1779 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
1780 $scale = $dbd+$scale;
1786 # round left from the '.'
1788 # 123 => 100 means length(123) = 3 - $scale (2) => 1
1790 my $dbt = $x->{_m}->length();
1792 my $dbd = $dbt + $x->{_e}->numify();
1793 # should be the same, so treat it as this
1794 $scale = 1 if $scale == 0;
1795 # shortcut if already integer
1796 return $x if $scale == 1 && $dbt <= $dbd;
1797 # maximum digits before dot
1802 # not enough digits before dot, so round to zero
1805 elsif ( $scale == $dbd )
1812 $scale = $dbd - $scale;
1815 # pass sign to bround for rounding modes '+inf' and '-inf'
1816 $x->{_m}->{sign} = $x->{sign};
1817 $x->{_m}->bround($scale,$mode);
1818 $x->{_m}->{sign} = '+'; # fix sign back
1824 # accuracy: preserve $N digits, and overwrite the rest with 0's
1825 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1827 if (($_[0] || 0) < 0)
1829 require Carp; Carp::croak ('bround() needs positive accuracy');
1832 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1833 return $x if !defined $scale; # no-op
1835 return $x if $x->modify('bround');
1837 # scale is now either $x->{_a}, $accuracy, or the user parameter
1838 # test whether $x already has lower accuracy, do nothing in this case
1839 # but do round if the accuracy is the same, since a math operation might
1840 # want to round a number with A=5 to 5 digits afterwards again
1841 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1843 # scale < 0 makes no sense
1844 # never round a +-inf, NaN
1845 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
1847 # 1: $scale == 0 => keep all digits
1848 # 2: never round a 0
1849 # 3: if we should keep more digits than the mantissa has, do nothing
1850 if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1852 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1856 # pass sign to bround for '+inf' and '-inf' rounding modes
1857 $x->{_m}->{sign} = $x->{sign};
1858 $x->{_m}->bround($scale,$mode); # round mantissa
1859 $x->{_m}->{sign} = '+'; # fix sign back
1860 # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1861 $x->{_a} = $scale; # remember rounding
1862 $x->{_p} = undef; # and clear P
1863 $x->bnorm(); # del trailing zeros gen. by bround()
1868 # return integer less or equal then $x
1869 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1871 return $x if $x->modify('bfloor');
1873 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1875 # if $x has digits after dot
1876 if ($x->{_e}->{sign} eq '-')
1878 $x->{_e}->{sign} = '+'; # negate e
1879 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1880 $x->{_e}->bzero(); # trunc/norm
1881 $x->{_m}->binc() if $x->{sign} eq '-'; # decrement if negative
1883 $x->round($a,$p,$r);
1888 # return integer greater or equal then $x
1889 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1891 return $x if $x->modify('bceil');
1892 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1894 # if $x has digits after dot
1895 if ($x->{_e}->{sign} eq '-')
1897 #$x->{_m}->brsft(-$x->{_e},10);
1899 #$x++ if $x->{sign} eq '+';
1901 $x->{_e}->{sign} = '+'; # negate e
1902 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1903 $x->{_e}->bzero(); # trunc/norm
1904 $x->{_m}->binc() if $x->{sign} eq '+'; # decrement if negative
1906 $x->round($a,$p,$r);
1911 # shift right by $y (divide by power of $n)
1914 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1915 # objectify is costly, so avoid it
1916 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1918 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1921 return $x if $x->modify('brsft');
1922 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1924 $n = 2 if !defined $n; $n = $self->new($n);
1925 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
1930 # shift left by $y (multiply by power of $n)
1933 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1934 # objectify is costly, so avoid it
1935 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1937 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1940 return $x if $x->modify('blsft');
1941 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1943 $n = 2 if !defined $n; $n = $self->new($n);
1944 $x->bmul($n->bpow($y),$a,$p,$r,$y);
1947 ###############################################################################
1951 # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1956 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
1957 # or falling back to MBI::bxxx()
1958 my $name = $AUTOLOAD;
1960 $name =~ s/.*:://; # split package
1962 $class->import() if $IMPORT == 0;
1963 if (!method_alias($name))
1967 # delayed load of Carp and avoid recursion
1969 Carp::croak ("Can't call a method without name");
1971 if (!method_hand_up($name))
1973 # delayed load of Carp and avoid recursion
1975 Carp::croak ("Can't call $class\-\>$name, not a valid method");
1977 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1979 return &{"$MBI"."::$name"}(@_);
1981 my $bname = $name; $bname =~ s/^f/b/;
1982 *{$class."::$name"} = \&$bname;
1988 # return a copy of the exponent
1989 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1991 if ($x->{sign} !~ /^[+-]$/)
1993 my $s = $x->{sign}; $s =~ s/^[+-]//;
1994 return $self->new($s); # -inf, +inf => +inf
1996 return $x->{_e}->copy();
2001 # return a copy of the mantissa
2002 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2004 if ($x->{sign} !~ /^[+-]$/)
2006 my $s = $x->{sign}; $s =~ s/^[+]//;
2007 return $self->new($s); # -inf, +inf => +inf
2009 my $m = $x->{_m}->copy(); # faster than going via bstr()
2010 $m->bneg() if $x->{sign} eq '-';
2017 # return a copy of both the exponent and the mantissa
2018 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2020 if ($x->{sign} !~ /^[+-]$/)
2022 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2023 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2025 my $m = $x->{_m}->copy(); # faster than going via bstr()
2026 $m->bneg() if $x->{sign} eq '-';
2027 return ($m,$x->{_e}->copy());
2030 ##############################################################################
2031 # private stuff (internal use only)
2037 my $lib = ''; my @a;
2039 for ( my $i = 0; $i < $l ; $i++)
2041 if ( $_[$i] eq ':constant' )
2043 # This causes overlord er load to step in. 'binary' and 'integer'
2044 # are handled by BigInt.
2045 overload::constant float => sub { $self->new(shift); };
2047 elsif ($_[$i] eq 'upgrade')
2049 # this causes upgrading
2050 $upgrade = $_[$i+1]; # or undef to disable
2053 elsif ($_[$i] eq 'downgrade')
2055 # this causes downgrading
2056 $downgrade = $_[$i+1]; # or undef to disable
2059 elsif ($_[$i] eq 'lib')
2061 # alternative library
2062 $lib = $_[$i+1] || ''; # default Calc
2065 elsif ($_[$i] eq 'with')
2067 # alternative class for our private parts()
2068 $MBI = $_[$i+1] || 'Math::BigInt'; # default Math::BigInt
2077 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2078 my $mbilib = eval { Math::BigInt->config()->{lib} };
2079 if ((defined $mbilib) && ($MBI eq 'Math::BigInt'))
2081 # MBI already loaded
2082 $MBI->import('lib',"$lib,$mbilib", 'objectify');
2086 # MBI not loaded, or with ne "Math::BigInt"
2087 $lib .= ",$mbilib" if defined $mbilib;
2088 $lib =~ s/^,//; # don't leave empty
2089 # replacement library can handle lib statement, but also could ignore it
2092 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2093 # used in the same script, or eval inside import().
2094 my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
2095 my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
2097 $file = File::Spec->catfile (@parts, $file);
2098 eval { require "$file"; };
2099 $MBI->import( lib => $lib, 'objectify' );
2103 my $rc = "use $MBI lib => '$lib', 'objectify';";
2109 require Carp; Carp::croak ("Couldn't load $MBI: $! $@");
2112 # any non :constant stuff is handled by our parent, Exporter
2113 # even if @_ is empty, to give it a chance
2114 $self->SUPER::import(@a); # for subclasses
2115 $self->export_to_level(1,$self,@a); # need this, too
2120 # adjust m and e so that m is smallest possible
2121 # round number according to accuracy and precision settings
2122 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2124 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2126 # if (!$x->{_m}->is_odd())
2128 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
2131 $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros);
2133 # for something like 0Ey, set y to 1, and -0 => +0
2134 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
2136 # this is to prevent automatically rounding when MBI's globals are set
2137 $x->{_m}->{_f} = MB_NEVER_ROUND;
2138 $x->{_e}->{_f} = MB_NEVER_ROUND;
2139 # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
2140 $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
2141 $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
2142 $x; # MBI bnorm is no-op, so dont call it
2145 ##############################################################################
2149 # return number as hexadecimal string (only for integers defined)
2150 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2152 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2153 return '0x0' if $x->is_zero();
2155 return $nan if $x->{_e}->{sign} ne '+'; # how to do 1e-1 in hex!?
2157 my $z = $x->{_m}->copy();
2158 if (!$x->{_e}->is_zero()) # > 0
2160 $z->blsft($x->{_e},10);
2162 $z->{sign} = $x->{sign};
2168 # return number as binary digit string (only for integers defined)
2169 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2171 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2172 return '0b0' if $x->is_zero();
2174 return $nan if $x->{_e}->{sign} ne '+'; # how to do 1e-1 in hex!?
2176 my $z = $x->{_m}->copy();
2177 if (!$x->{_e}->is_zero()) # > 0
2179 $z->blsft($x->{_e},10);
2181 $z->{sign} = $x->{sign};
2187 # return copy as a bigint representation of this BigFloat number
2188 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2190 my $z = $x->{_m}->copy();
2191 if ($x->{_e}->{sign} eq '-') # < 0
2193 $x->{_e}->{sign} = '+'; # flip
2194 $z->brsft($x->{_e},10);
2195 $x->{_e}->{sign} = '-'; # flip back
2197 elsif (!$x->{_e}->is_zero()) # > 0
2199 $z->blsft($x->{_e},10);
2201 $z->{sign} = $x->{sign};
2208 my $class = ref($x) || $x;
2209 $x = $class->new(shift) unless ref($x);
2211 return 1 if $x->{_m}->is_zero();
2212 my $len = $x->{_m}->length();
2213 $len += $x->{_e} if $x->{_e}->sign() eq '+';
2216 my $t = $MBI->bzero();
2217 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
2228 Math::BigFloat - Arbitrary size floating point math package
2235 $x = Math::BigFloat->new($str); # defaults to 0
2236 $nan = Math::BigFloat->bnan(); # create a NotANumber
2237 $zero = Math::BigFloat->bzero(); # create a +0
2238 $inf = Math::BigFloat->binf(); # create a +inf
2239 $inf = Math::BigFloat->binf('-'); # create a -inf
2240 $one = Math::BigFloat->bone(); # create a +1
2241 $one = Math::BigFloat->bone('-'); # create a -1
2244 $x->is_zero(); # true if arg is +0
2245 $x->is_nan(); # true if arg is NaN
2246 $x->is_one(); # true if arg is +1
2247 $x->is_one('-'); # true if arg is -1
2248 $x->is_odd(); # true if odd, false for even
2249 $x->is_even(); # true if even, false for odd
2250 $x->is_positive(); # true if >= 0
2251 $x->is_negative(); # true if < 0
2252 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2254 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2255 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2256 $x->sign(); # return the sign, either +,- or NaN
2257 $x->digit($n); # return the nth digit, counting from right
2258 $x->digit(-$n); # return the nth digit, counting from left
2260 # The following all modify their first argument. If you want to preserve
2261 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2262 # neccessary when mixing $a = $b assigments with non-overloaded math.
2265 $x->bzero(); # set $i to 0
2266 $x->bnan(); # set $i to NaN
2267 $x->bone(); # set $x to +1
2268 $x->bone('-'); # set $x to -1
2269 $x->binf(); # set $x to inf
2270 $x->binf('-'); # set $x to -inf
2272 $x->bneg(); # negation
2273 $x->babs(); # absolute value
2274 $x->bnorm(); # normalize (no-op)
2275 $x->bnot(); # two's complement (bit wise not)
2276 $x->binc(); # increment x by 1
2277 $x->bdec(); # decrement x by 1
2279 $x->badd($y); # addition (add $y to $x)
2280 $x->bsub($y); # subtraction (subtract $y from $x)
2281 $x->bmul($y); # multiplication (multiply $x by $y)
2282 $x->bdiv($y); # divide, set $x to quotient
2283 # return (quo,rem) or quo if scalar
2285 $x->bmod($y); # modulus ($x % $y)
2286 $x->bpow($y); # power of arguments ($x ** $y)
2287 $x->blsft($y); # left shift
2288 $x->brsft($y); # right shift
2289 # return (quo,rem) or quo if scalar
2291 $x->blog(); # logarithm of $x to base e (Euler's number)
2292 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
2294 $x->band($y); # bit-wise and
2295 $x->bior($y); # bit-wise inclusive or
2296 $x->bxor($y); # bit-wise exclusive or
2297 $x->bnot(); # bit-wise not (two's complement)
2299 $x->bsqrt(); # calculate square-root
2300 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
2301 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2303 $x->bround($N); # accuracy: preserve $N digits
2304 $x->bfround($N); # precision: round to the $Nth digit
2306 $x->bfloor(); # return integer less or equal than $x
2307 $x->bceil(); # return integer greater or equal than $x
2309 # The following do not modify their arguments:
2311 bgcd(@values); # greatest common divisor
2312 blcm(@values); # lowest common multiplicator
2314 $x->bstr(); # return string
2315 $x->bsstr(); # return string in scientific notation
2317 $x->exponent(); # return exponent as BigInt
2318 $x->mantissa(); # return mantissa as BigInt
2319 $x->parts(); # return (mantissa,exponent) as BigInt
2321 $x->length(); # number of digits (w/o sign and '.')
2322 ($l,$f) = $x->length(); # number of digits, and length of fraction
2324 $x->precision(); # return P of $x (or global, if P of $x undef)
2325 $x->precision($n); # set P of $x to $n
2326 $x->accuracy(); # return A of $x (or global, if A of $x undef)
2327 $x->accuracy($n); # set A $x to $n
2329 # these get/set the appropriate global value for all BigFloat objects
2330 Math::BigFloat->precision(); # Precision
2331 Math::BigFloat->accuracy(); # Accuracy
2332 Math::BigFloat->round_mode(); # rounding mode
2336 All operators (inlcuding basic math operations) are overloaded if you
2337 declare your big floating point numbers as
2339 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2341 Operations with overloaded operators preserve the arguments, which is
2342 exactly what you expect.
2344 =head2 Canonical notation
2346 Input to these routines are either BigFloat objects, or strings of the
2347 following four forms:
2361 C</^[+-]\d+E[+-]?\d+$/>
2365 C</^[+-]\d*\.\d+E[+-]?\d+$/>
2369 all with optional leading and trailing zeros and/or spaces. Additonally,
2370 numbers are allowed to have an underscore between any two digits.
2372 Empty strings as well as other illegal numbers results in 'NaN'.
2374 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2375 are always stored in normalized form. On a string, it creates a BigFloat
2380 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2382 The string output will always have leading and trailing zeros stripped and drop
2383 a plus sign. C<bstr()> will give you always the form with a decimal point,
2384 while C<bsstr()> (s for scientific) gives you the scientific notation.
2386 Input bstr() bsstr()
2388 ' -123 123 123' '-123123123' '-123123123E0'
2389 '00.0123' '0.0123' '123E-4'
2390 '123.45E-2' '1.2345' '12345E-4'
2391 '10E+3' '10000' '1E4'
2393 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2394 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2395 return either undef, <0, 0 or >0 and are suited for sort.
2397 Actual math is done by using the class defined with C<with => Class;> (which
2398 defaults to BigInts) to represent the mantissa and exponent.
2400 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2401 represent the result when input arguments are not numbers, as well as
2402 the result of dividing by zero.
2404 =head2 C<mantissa()>, C<exponent()> and C<parts()>
2406 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2407 as BigInts such that:
2409 $m = $x->mantissa();
2410 $e = $x->exponent();
2411 $y = $m * ( 10 ** $e );
2412 print "ok\n" if $x == $y;
2414 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2416 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2418 Currently the mantissa is reduced as much as possible, favouring higher
2419 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2420 This might change in the future, so do not depend on it.
2422 =head2 Accuracy vs. Precision
2424 See also: L<Rounding|Rounding>.
2426 Math::BigFloat supports both precision and accuracy. For a full documentation,
2427 examples and tips on these topics please see the large section in
2430 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2431 a operation consumes all resources, each operation produces no more than
2432 the requested number of digits.
2434 Please refer to BigInt's documentation for the precedence rules of which
2435 accuracy/precision setting will be used.
2437 If there is no gloabl precision set, B<and> the operation inquestion was not
2438 called with a requested precision or accuracy, B<and> the input $x has no
2439 accuracy or precision set, then a fallback parameter will be used. For
2440 historical reasons, it is called C<div_scale> and can be accessed via:
2442 $d = Math::BigFloat->div_scale(); # query
2443 Math::BigFloat->div_scale($n); # set to $n digits
2445 The default value is 40 digits.
2447 In case the result of one operation has more precision than specified,
2448 it is rounded. The rounding mode taken is either the default mode, or the one
2449 supplied to the operation after the I<scale>:
2451 $x = Math::BigFloat->new(2);
2452 Math::BigFloat->precision(5); # 5 digits max
2453 $y = $x->copy()->bdiv(3); # will give 0.66666
2454 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2455 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
2456 Math::BigFloat->round_mode('zero');
2457 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2463 =item ffround ( +$scale )
2465 Rounds to the $scale'th place left from the '.', counting from the dot.
2466 The first digit is numbered 1.
2468 =item ffround ( -$scale )
2470 Rounds to the $scale'th place right from the '.', counting from the dot.
2474 Rounds to an integer.
2476 =item fround ( +$scale )
2478 Preserves accuracy to $scale digits from the left (aka significant digits)
2479 and pads the rest with zeros. If the number is between 1 and -1, the
2480 significant digits count from the first non-zero after the '.'
2482 =item fround ( -$scale ) and fround ( 0 )
2484 These are effectively no-ops.
2488 All rounding functions take as a second parameter a rounding mode from one of
2489 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2491 The default rounding mode is 'even'. By using
2492 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
2493 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2494 no longer supported.
2495 The second parameter to the round functions then overrides the default
2498 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
2499 'trunc' as rounding mode to make it equivalent to:
2504 You can override this by passing the desired rounding mode as parameter to
2507 $x = Math::BigFloat->new(2.5);
2508 $y = $x->as_number('odd'); # $y = 3
2514 =head1 Autocreating constants
2516 After C<use Math::BigFloat ':constant'> all the floating point constants
2517 in the given scope are converted to C<Math::BigFloat>. This conversion
2518 happens at compile time.
2522 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2524 prints the value of C<2E-100>. Note that without conversion of
2525 constants the expression 2E-100 will be calculated as normal floating point
2528 Please note that ':constant' does not affect integer constants, nor binary
2529 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2534 Math with the numbers is done (by default) by a module called
2535 Math::BigInt::Calc. This is equivalent to saying:
2537 use Math::BigFloat lib => 'Calc';
2539 You can change this by using:
2541 use Math::BigFloat lib => 'BitVect';
2543 The following would first try to find Math::BigInt::Foo, then
2544 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2546 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2548 Calc.pm uses as internal format an array of elements of some decimal base
2549 (usually 1e7, but this might be differen for some systems) with the least
2550 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2551 significant bit first. Other modules might use even different means of
2552 representing the numbers. See the respective module documentation for further
2555 Please note that Math::BigFloat does B<not> use the denoted library itself,
2556 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2559 use Math::BigInt lib => 'GMP';
2562 you can roll it all into one line:
2564 use Math::BigFloat lib => 'GMP';
2566 It is also possible to just require Math::BigFloat:
2568 require Math::BigFloat;
2570 This will load the neccessary things (like BigInt) when they are needed, and
2573 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2574 you ever wanted to know about loading a different library.
2576 =head2 Using Math::BigInt::Lite
2578 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2581 use Math::BigFloat with => 'Math::BigInt::Lite';
2583 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2584 can combine these if you want. For instance, you may want to use
2585 Math::BigInt objects in your main script, too.
2589 use Math::BigFloat with => 'Math::BigInt::Lite';
2591 Of course, you can combine this with the C<lib> parameter.
2594 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2596 There is no need for a "use Math::BigInt;" statement, even if you want to
2597 use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2598 always loads it. But if you add it, add it B<before>:
2602 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2604 Notice that the module with the last C<lib> will "win" and thus
2605 it's lib will be used if the lib is available:
2608 use Math::BigInt lib => 'Bar,Baz';
2609 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2611 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2612 words, Math::BigFloat will try to retain previously loaded libs when you
2613 don't specify it onem but if you specify one, it will try to load them.
2615 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2616 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2617 same as trying the latter load alone, except for the fact that one of Bar or
2618 Baz might be loaded needlessly in an intermidiate step (and thus hang around
2619 and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2620 will still be tried to be loaded, but this is not as time/memory consuming as
2621 actually loading one of them. Still, this type of usage is not recommended due
2624 The old way (loading the lib only in BigInt) still works though:
2627 use Math::BigInt lib => 'Bar,Baz';
2630 You can even load Math::BigInt afterwards:
2634 use Math::BigInt lib => 'Bar,Baz';
2636 But this has the same problems like #5, it will first load Calc
2637 (Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2638 Baz, depending on which of them works and is usable/loadable. Since this
2639 loads Calc unnecc., it is not recommended.
2641 Since it also possible to just require Math::BigFloat, this poses the question
2642 about what libary this will use:
2644 require Math::BigFloat;
2645 my $x = Math::BigFloat->new(123); $x += 123;
2647 It will use Calc. Please note that the call to import() is still done, but
2648 only when you use for the first time some Math::BigFloat math (it is triggered
2649 via any constructor, so the first time you create a Math::BigFloat, the load
2650 will happen in the background). This means:
2652 require Math::BigFloat;
2653 Math::BigFloat->import ( lib => 'Foo,Bar' );
2655 would be the same as:
2657 use Math::BigFloat lib => 'Foo, Bar';
2659 But don't try to be clever to insert some operations in between:
2661 require Math::BigFloat;
2662 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
2663 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
2664 $x = Math::BigFloat->bone()+4; # now use Pari
2666 While this works, it loads Calc needlessly. But maybe you just wanted that?
2668 B<Examples #3 is highly recommended> for daily usage.
2672 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
2678 =item stringify, bstr()
2680 Both stringify and bstr() now drop the leading '+'. The old code would return
2681 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2682 reasoning and details.
2686 The following will probably not do what you expect:
2688 print $c->bdiv(123.456),"\n";
2690 It prints both quotient and reminder since print works in list context. Also,
2691 bdiv() will modify $c, so be carefull. You probably want to use
2693 print $c / 123.456,"\n";
2694 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2698 =item Modifying and =
2702 $x = Math::BigFloat->new(5);
2705 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2706 a second reference to the B<same> object and stores it in $y. Thus anything
2707 that modifies $x will modify $y (except overloaded math operators), and vice
2708 versa. See L<Math::BigInt> for details and how to avoid that.
2712 C<bpow()> now modifies the first argument, unlike the old code which left
2713 it alone and only returned the result. This is to be consistent with
2714 C<badd()> etc. The first will modify $x, the second one won't:
2716 print bpow($x,$i),"\n"; # modify $x
2717 print $x->bpow($i),"\n"; # ditto
2718 print $x ** $i,"\n"; # leave $x alone
2724 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
2725 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
2727 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
2728 because they solve the autoupgrading/downgrading issue, at least partly.
2731 L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
2732 more documentation including a full version history, testcases, empty
2733 subclass files and benchmarks.
2737 This program is free software; you may redistribute it and/or modify it under
2738 the same terms as Perl itself.
2742 Mark Biggar, overloaded interface by Ilya Zakharevich.
2743 Completely rewritten by Tels http://bloodgate.com in 2001, 2002, and still