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 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
176 $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;
177 $self->{sign} = $$mis;
179 # if downgrade, inf, NaN or integers go down
181 if ($downgrade && $self->{_e}->{sign} eq '+')
183 #print "downgrading $$miv$$mfv"."E$$es$$ev";
184 if ($self->{_e}->is_zero())
186 $self->{_m}->{sign} = $$mis; # negative if wanted
187 return $downgrade->new($self->{_m});
189 return $downgrade->new("$$mis$$miv$$mfv"."E$$es$$ev");
191 #print "mbf new $self->{sign} $self->{_m} e $self->{_e} ",ref($self),"\n";
192 $self->bnorm()->round(@r); # first normalize, then round
197 # used by parent class bone() to initialize number to NaN
203 my $class = ref($self);
204 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
207 $IMPORT=1; # call our import only once
208 $self->{_m} = $MBI->bzero();
209 $self->{_e} = $MBI->bzero();
214 # used by parent class bone() to initialize number to +-inf
220 my $class = ref($self);
221 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
224 $IMPORT=1; # call our import only once
225 $self->{_m} = $MBI->bzero();
226 $self->{_e} = $MBI->bzero();
231 # used by parent class bone() to initialize number to 1
233 $IMPORT=1; # call our import only once
234 $self->{_m} = $MBI->bone();
235 $self->{_e} = $MBI->bzero();
240 # used by parent class bone() to initialize number to 0
242 $IMPORT=1; # call our import only once
243 $self->{_m} = $MBI->bzero();
244 $self->{_e} = $MBI->bone();
249 my ($self,$class) = @_;
250 return if $class =~ /^Math::BigInt/; # we aren't one of these
251 UNIVERSAL::isa($self,$class);
256 # return (later set?) configuration data as hash ref
257 my $class = shift || 'Math::BigFloat';
259 my $cfg = $class->SUPER::config(@_);
261 # now we need only to override the ones that are different from our parent
262 $cfg->{class} = $class;
267 ##############################################################################
268 # string conversation
272 # (ref to BFLOAT or num_str ) return num_str
273 # Convert number from internal format to (non-scientific) string format.
274 # internal format is always normalized (no leading zeros, "-0" => "+0")
275 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
277 if ($x->{sign} !~ /^[+-]$/)
279 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
283 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
286 my $not_zero = !($x->{sign} eq '+' && $x->{_m}->is_zero());
289 $es = $x->{_m}->bstr();
290 $len = CORE::length($es);
291 my $e = $x->{_e}->numify();
295 # if _e is bigger than a scalar, the following will blow your memory
298 #print "style: 0.xxxx\n";
299 my $r = abs($e) - $len;
300 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
304 #print "insert '.' at $e in '$es'\n";
305 substr($es,$e,0) = '.'; $cad = $x->{_e};
311 $es .= '0' x $e; $len += $e; $cad = 0;
314 $es = '-'.$es if $x->{sign} eq '-';
315 # if set accuracy or precision, pad with zeros on the right side
316 if ((defined $x->{_a}) && ($not_zero))
318 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
319 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
320 $zeros = $x->{_a} - $len if $cad != $len;
321 $es .= $dot.'0' x $zeros if $zeros > 0;
323 elsif ((($x->{_p} || 0) < 0))
325 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
326 my $zeros = -$x->{_p} + $cad;
327 $es .= $dot.'0' x $zeros if $zeros > 0;
334 # (ref to BFLOAT or num_str ) return num_str
335 # Convert number from internal format to scientific string format.
336 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
337 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
338 #my $x = shift; my $class = ref($x) || $x;
339 #$x = $class->new(shift) unless ref($x);
341 if ($x->{sign} !~ /^[+-]$/)
343 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
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)
543 my $add = $y->{_m}->copy();
544 if ($e->{sign} eq '-') # < 0
546 my $e1 = $e->copy()->babs();
547 #$x->{_m} *= (10 ** $e1);
548 $x->{_m}->blsft($e1,10);
549 $x->{_e} += $e; # need the sign of e
551 elsif (!$e->is_zero()) # > 0
556 # else: both e are the same, so just leave them
557 $x->{_m}->{sign} = $x->{sign}; # fiddle with signs
558 $add->{sign} = $y->{sign};
559 $x->{_m} += $add; # finally do add/sub
560 $x->{sign} = $x->{_m}->{sign}; # re-adjust signs
561 $x->{_m}->{sign} = '+'; # mantissa always positiv
562 # delete trailing zeros, then round
563 return $x->bnorm()->round($a,$p,$r,$y);
568 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
569 # subtract second arg from first, modify first
572 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
573 # objectify is costly, so avoid it
574 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
576 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
579 if ($y->is_zero()) # still round for not adding zero
581 return $x->round($a,$p,$r);
584 $y->{sign} =~ tr/+\-/-+/; # does nothing for NaN
585 $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
586 $y->{sign} =~ tr/+\-/-+/; # refix $y (does nothing for NaN)
587 $x; # already rounded by badd()
592 # increment arg by one
593 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
595 if ($x->{_e}->sign() eq '-')
597 return $x->badd($self->bone(),$a,$p,$r); # digits after dot
600 if (!$x->{_e}->is_zero())
602 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
606 if ($x->{sign} eq '+')
609 return $x->bnorm()->bround($a,$p,$r);
611 elsif ($x->{sign} eq '-')
614 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
615 return $x->bnorm()->bround($a,$p,$r);
617 # inf, nan handling etc
618 $x->badd($self->__one(),$a,$p,$r); # does round
623 # decrement arg by one
624 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
626 if ($x->{_e}->sign() eq '-')
628 return $x->badd($self->bone('-'),$a,$p,$r); # digits after dot
631 if (!$x->{_e}->is_zero())
633 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
637 my $zero = $x->is_zero();
639 if (($x->{sign} eq '-') || $zero)
642 $x->{sign} = '-' if $zero; # 0 => 1 => -1
643 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
644 return $x->bnorm()->round($a,$p,$r);
647 elsif ($x->{sign} eq '+')
650 return $x->bnorm()->round($a,$p,$r);
652 # inf, nan handling etc
653 $x->badd($self->bone('-'),$a,$p,$r); # does round
660 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
662 # $base > 0, $base != 1; if $base == undef default to $base == e
665 # we need to limit the accuracy to protect against overflow
668 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
670 # also takes care of the "error in _find_round_parameters?" case
671 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
673 # no rounding at all, so must use fallback
674 if (scalar @params == 0)
676 # simulate old behaviour
677 $params[0] = $self->div_scale(); # and round to it as accuracy
678 $params[1] = undef; # P = undef
679 $scale = $params[0]+4; # at least four more for proper round
680 $params[2] = $r; # round mode by caller or undef
681 $fallback = 1; # to clear a/p afterwards
685 # the 4 below is empirical, and there might be cases where it is not
687 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
690 return $x->bzero(@params) if $x->is_one();
691 # base not defined => base == Euler's constant e
694 # make object, since we don't feed it trough objectify() to still get the
695 # case of $base == undef
696 $base = $self->new($base) unless ref($base);
697 # $base > 0; $base != 1
698 return $x->bnan() if $base->is_zero() || $base->is_one() ||
699 $base->{sign} ne '+';
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
729 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
732 # shortcut to not run trough _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 log based on Taylor.
755 # Modifies $x in place.
756 my ($self,$x,$scale) = @_;
758 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
762 # Taylor: | u 1 u^3 1 u^5 |
763 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
764 # |_ v 3 v^3 5 v^5 _|
766 # This takes much more steps to calculate the result and is thus not used
769 # Taylor: | u 1 u^2 1 u^3 |
770 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
771 # |_ x 2 x^2 3 x^3 _|
773 # "normal" log algorithmn
775 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
777 $v = $x->copy(); $v->binc(); # v = x+1
778 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
779 $x->bdiv($v,$scale); # first term: u/v
782 $u *= $u; $v *= $v; # u^2, v^2
783 $below->bmul($v); # u^3, v^3
785 $factor = $self->new(3); $f = $self->new(2);
787 my $steps = 0 if DEBUG;
788 $limit = $self->new("1E-". ($scale-1));
791 # we calculate the next term, and add it to the last
792 # when the next term is below our limit, it won't affect the outcome
793 # anymore, so we stop
795 # calculating the next term simple from over/below will result in quite
796 # a time hog if the input has many digits, since over and below will
797 # accumulate more and more digits, and the result will also have many
798 # digits, but in the end it is rounded to $scale digits anyway. So if we
799 # round $over and $below first, we save a lot of time for the division
800 # (not with log(1.2345), but try log (123**123) to see what I mean. This
801 # can introduce a rounding error if the division result would be f.i.
802 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
803 # if we truncated the $over and $below we might get 0.12345. Does this
804 # matter for the end result? So we give over and below 4 more digits to be
805 # on the safe side (unscientific error handling as usual...)
806 # Makes blog(1.23) *slightly* slower, but try blog(123*123) w/o it :o)
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 my ($self,$x,$scale) = @_;
836 # taking blog() from numbers greater than 10 takes a *very long* time, so we
837 # break the computation down into parts based on the observation that:
838 # blog(x*y) = blog(x) + blog(y)
839 # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
840 # the faster it get's, especially because 2*$x takes about 10 times as long,
841 # so by dividing $x by 10 we make it at least factor 100 faster...)
843 # The same observation is valid for numbers smaller than 0.1 (e.g. computing
844 # log(1) is fastest, and the farther away we get from 1, the longer it takes)
845 # so we also 'break' this down by multiplying $x with 10 and subtract the
846 # log(10) afterwards to get the correct result.
848 # calculate nr of digits before dot
849 my $dbd = $x->{_m}->length() + $x->{_e}->numify();
851 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
854 my $calc = 1; # do some calculation?
856 # disable the shortcut for 10, since we need log(10) and this would recurse
858 if ($x->{_e}->is_one() && $x->{_m}->is_one())
860 $dbd = 0; # disable shortcut
861 # we can use the cached value in these cases
862 if ($scale <= $LOG_10_A)
864 $x->bzero(); $x->badd($LOG_10);
865 $calc = 0; # no need to calc, but round
868 # disable the shortcut for 2, since we maybe have it cached
869 my $two = $self->new(2); # also used later on
870 if ($x->{_e}->is_zero() && $x->{_m}->bcmp($two) == 0)
872 $dbd = 0; # disable shortcut
873 # we can use the cached value in these cases
874 if ($scale <= $LOG_2_A)
876 $x->bzero(); $x->badd($LOG_2);
877 $calc = 0; # no need to calc, but round
881 # if $x = 0.1, we know the result must be 0-log(10)
882 if ($x->{_e}->is_one('-') && $x->{_m}->is_one())
884 $dbd = 0; # disable shortcut
885 # we can use the cached value in these cases
886 if ($scale <= $LOG_10_A)
888 $x->bzero(); $x->bsub($LOG_10);
889 $calc = 0; # no need to calc, but round
893 # default: these correction factors are undef and thus not used
894 my $l_10; # value of ln(10) to A of $scale
895 my $l_2; # value of ln(2) to A of $scale
897 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
898 # so don't do this shortcut for 1 or 0
899 if (($dbd > 1) || ($dbd < 0))
901 # convert our cached value to an object if not already (avoid doing this
902 # at import() time, since not everybody needs this)
903 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
905 #print "x = $x, dbd = $dbd, calc = $calc\n";
906 # got more than one digit before the dot, or more than one zero after the
908 # log(123) == log(1.23) + log(10) * 2
909 # log(0.0123) == log(1.23) - log(10) * 2
911 if ($scale <= $LOG_10_A)
914 #print "using cached value for l_10\n";
915 $l_10 = $LOG_10->copy(); # copy for mul
919 # else: slower, compute it (but don't cache it, because it could be big)
920 # also disable downgrade for this code path
921 local $Math::BigFloat::downgrade = undef;
922 #print "l_10 = $l_10 (self = $self',
923 # ", ref(l_10) = ",ref($l_10)," scale $scale)\n";
924 #print "calculating value for l_10, scale $scale\n";
925 $l_10 = $self->new(10)->blog(undef,$scale); # scale+4, actually
927 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
929 $dbd = $self->new($dbd);
931 $l_10->bmul($dbd); # log(10) * (digits_before_dot-1)
932 #print "l_10 = $l_10\n";
934 $x->{_e}->bsub($dbd); # 123 => 1.23
936 #print "calculating log($x) with scale=$scale\n";
940 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
942 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
943 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
947 my $half = $self->new('0.5');
948 my $twos = 0; # default: none (0 times)
949 while ($x->bacmp($half) < 0)
952 $twos--; $x->bmul($two);
954 while ($x->bacmp($two) > 0)
957 $twos++; $x->bdiv($two,$scale+4); # keep all digits
960 # $twos > 0 => did mul 2, < 0 => did div 2 (never both)
961 # calculate correction factor based on ln(2)
964 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
965 if ($scale <= $LOG_2_A)
968 #print "using cached value for l_10\n";
969 $l_2 = $LOG_2->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 "calculating value for l_2, scale $scale\n";
977 $l_2 = $two->blog(undef,$scale); # scale+4, actually
980 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
987 $self->_log($x,$scale); # need to do the "normal" way
988 #print "log(x) = $x\n";
989 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
990 #print "result = $x\n";
991 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
992 #print "result = $x\n";
994 # all done, $x contains now the result
999 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1000 # does not modify arguments, but returns new object
1001 # Lowest Common Multiplicator
1003 my ($self,@arg) = objectify(0,@_);
1004 my $x = $self->new(shift @arg);
1005 while (@arg) { $x = _lcm($x,shift @arg); }
1011 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
1012 # does not modify arguments, but returns new object
1013 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
1015 my ($self,@arg) = objectify(0,@_);
1016 my $x = $self->new(shift @arg);
1017 while (@arg) { $x = _gcd($x,shift @arg); }
1021 ###############################################################################
1022 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1026 # return true if arg (BFLOAT or num_str) is an integer
1027 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1029 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1030 $x->{_e}->{sign} eq '+'; # 1e-1 => no integer
1036 # return true if arg (BFLOAT or num_str) is zero
1037 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1039 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
1045 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1046 my ($self,$x,$sign) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1048 $sign = '+' if !defined $sign || $sign ne '-';
1050 if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
1056 # return true if arg (BFLOAT or num_str) is odd or false if even
1057 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1059 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1060 ($x->{_e}->is_zero() && $x->{_m}->is_odd());
1066 # return true if arg (BINT or num_str) is even or false if odd
1067 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1069 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
1070 return 1 if ($x->{_e}->{sign} eq '+' # 123.45 is never
1071 && $x->{_m}->is_even()); # but 1200 is
1077 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
1078 # (BINT or num_str, BINT or num_str) return BINT
1081 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1082 # objectify is costly, so avoid it
1083 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1085 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1088 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1091 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1093 return $x->bnan() if $x->is_zero() || $y->is_zero();
1094 # result will always be +-inf:
1095 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1096 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1097 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1098 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1099 return $x->binf('-');
1102 return $x->bzero() if $x->is_zero() || $y->is_zero();
1104 return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
1105 ((!$x->isa($self)) || (!$y->isa($self)));
1107 # aEb * cEd = (a*c)E(b+d)
1108 $x->{_m}->bmul($y->{_m});
1109 $x->{_e}->badd($y->{_e});
1111 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1112 return $x->bnorm()->round($a,$p,$r,$y);
1117 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1118 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1121 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1122 # objectify is costly, so avoid it
1123 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1125 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1128 return $self->_div_inf($x,$y)
1129 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1131 # x== 0 # also: or y == 1 or y == -1
1132 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1135 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1137 # we need to limit the accuracy to protect against overflow
1139 my (@params,$scale);
1140 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1142 return $x if $x->is_nan(); # error in _find_round_parameters?
1144 # no rounding at all, so must use fallback
1145 if (scalar @params == 0)
1147 # simulate old behaviour
1148 $params[0] = $self->div_scale(); # and round to it as accuracy
1149 $scale = $params[0]+4; # at least four more for proper round
1150 $params[2] = $r; # round mode by caller or undef
1151 $fallback = 1; # to clear a/p afterwards
1155 # the 4 below is empirical, and there might be cases where it is not
1157 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1159 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
1160 $scale = $lx if $lx > $scale;
1161 $scale = $ly if $ly > $scale;
1162 my $diff = $ly - $lx;
1163 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
1165 # make copy of $x in case of list context for later reminder calculation
1167 if (wantarray && !$y->is_one())
1172 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1174 # check for / +-1 ( +/- 1E0)
1177 # promote BigInts and it's subclasses (except when already a BigFloat)
1178 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1180 # need to disable $upgrade in BigInt, to avoid deep recursion
1181 local $Math::BigInt::upgrade = undef; # should be parent class vs MBI
1183 # calculate the result to $scale digits and then round it
1184 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1185 $x->{_m}->blsft($scale,10);
1186 $x->{_m}->bdiv( $y->{_m} ); # a/c
1187 $x->{_e}->bsub( $y->{_e} ); # b-d
1188 $x->{_e}->bsub($scale); # correct for 10**scale
1189 $x->bnorm(); # remove trailing 0's
1192 # shortcut to not run trough _find_round_parameters again
1193 if (defined $params[0])
1195 $x->{_a} = undef; # clear before round
1196 $x->bround($params[0],$params[2]); # then round accordingly
1200 $x->{_p} = undef; # clear before round
1201 $x->bfround($params[1],$params[2]); # then round accordingly
1205 # clear a/p after round, since user did not request it
1206 $x->{_a} = undef; $x->{_p} = undef;
1213 $rem->bmod($y,@params); # copy already done
1217 $rem = $self->bzero();
1221 # clear a/p after round, since user did not request it
1222 $rem->{_a} = undef; $rem->{_p} = undef;
1231 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1234 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1235 # objectify is costly, so avoid it
1236 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1238 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1241 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1243 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1244 $x->{sign} = $re->{sign};
1245 $x->{_e} = $re->{_e};
1246 $x->{_m} = $re->{_m};
1247 return $x->round($a,$p,$r,$y);
1249 return $x->bnan() if $x->is_zero() && $y->is_zero();
1250 return $x if $y->is_zero();
1251 return $x->bnan() if $x->is_nan() || $y->is_nan();
1252 return $x->bzero() if $y->is_one() || $x->is_zero();
1254 # inf handling is missing here
1256 my $cmp = $x->bacmp($y); # equal or $x < $y?
1257 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1259 # only $y of the operands negative?
1260 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1262 $x->{sign} = $y->{sign}; # calc sign first
1263 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1265 my $ym = $y->{_m}->copy();
1268 $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
1270 # if $y has digits after dot
1271 my $shifty = 0; # correct _e of $x by this
1272 if ($y->{_e}->{sign} eq '-') # has digits after dot
1274 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1275 $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
1276 $x->blsft($shifty,10); # 123 => 1230, $y->{_m} is already 25
1278 # $ym is now mantissa of $y based on exponent 0
1280 my $shiftx = 0; # correct _e of $x by this
1281 if ($x->{_e}->{sign} eq '-') # has digits after dot
1283 # 123.4 % 20 => 1234 % 200
1284 $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
1285 $ym->blsft($shiftx,10);
1287 # 123e1 % 20 => 1230 % 20
1288 if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
1290 $x->{_m}->blsft($x->{_e},10);
1292 $x->{_e} = $MBI->bzero() unless $x->{_e}->is_zero();
1294 $x->{_e}->bsub($shiftx) if $shiftx != 0;
1295 $x->{_e}->bsub($shifty) if $shifty != 0;
1297 # now mantissas are equalized, exponent of $x is adjusted, so calc result
1299 $x->{_m}->bmod($ym);
1301 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
1304 if ($neg != 0) # one of them negative => correct in place
1307 $x->{_m} = $r->{_m};
1308 $x->{_e} = $r->{_e};
1309 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
1313 $x->round($a,$p,$r,$y); # round and return
1318 # calculate $y'th root of $x
1319 my ($self,$x,$y,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(2,@_);
1321 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1322 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1323 $y->{sign} !~ /^\+$/;
1325 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1327 # we need to limit the accuracy to protect against overflow
1329 my (@params,$scale);
1330 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1332 return $x if $x->is_nan(); # error in _find_round_parameters?
1334 # no rounding at all, so must use fallback
1335 if (scalar @params == 0)
1337 # simulate old behaviour
1338 $params[0] = $self->div_scale(); # and round to it as accuracy
1339 $scale = $params[0]+4; # at least four more for proper round
1340 $params[2] = $r; # round mode by caller or undef
1341 $fallback = 1; # to clear a/p afterwards
1345 # the 4 below is empirical, and there might be cases where it is not
1347 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1350 # when user set globals, they would interfere with our calculation, so
1351 # disable them and later re-enable them
1353 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1354 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1355 # we also need to disable any set A or P on $x (_find_round_parameters took
1356 # them already into account), since these would interfere, too
1357 delete $x->{_a}; delete $x->{_p};
1358 # need to disable $upgrade in BigInt, to avoid deep recursion
1359 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1361 # remember sign and make $x positive, since -4 ** (1/2) => -2
1362 my $sign = 0; $sign = 1 if $x->is_negative(); $x->babs();
1364 if ($y->bcmp(2) == 0) # normal square root
1366 $x->bsqrt($scale+4);
1368 elsif ($y->is_one('-'))
1371 my $u = $self->bone()->bdiv($x,$scale);
1372 # copy private parts over
1373 $x->{_m} = $u->{_m};
1374 $x->{_e} = $u->{_e};
1378 my $u = $self->bone()->bdiv($y,$scale+4);
1379 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
1380 $x->bpow($u,$scale+4); # el cheapo
1382 $x->bneg() if $sign == 1;
1384 # shortcut to not run trough _find_round_parameters again
1385 if (defined $params[0])
1387 $x->bround($params[0],$params[2]); # then round accordingly
1391 $x->bfround($params[1],$params[2]); # then round accordingly
1395 # clear a/p after round, since user did not request it
1396 $x->{_a} = undef; $x->{_p} = undef;
1399 $$abr = $ab; $$pbr = $pb;
1405 # calculate square root
1406 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1408 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
1409 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
1410 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
1412 # we need to limit the accuracy to protect against overflow
1414 my (@params,$scale);
1415 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1417 return $x if $x->is_nan(); # error in _find_round_parameters?
1419 # no rounding at all, so must use fallback
1420 if (scalar @params == 0)
1422 # simulate old behaviour
1423 $params[0] = $self->div_scale(); # and round to it as accuracy
1424 $scale = $params[0]+4; # at least four more for proper round
1425 $params[2] = $r; # round mode by caller or undef
1426 $fallback = 1; # to clear a/p afterwards
1430 # the 4 below is empirical, and there might be cases where it is not
1432 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1435 # when user set globals, they would interfere with our calculation, so
1436 # disable them and later re-enable them
1438 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1439 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1440 # we also need to disable any set A or P on $x (_find_round_parameters took
1441 # them already into account), since these would interfere, too
1442 delete $x->{_a}; delete $x->{_p};
1443 # need to disable $upgrade in BigInt, to avoid deep recursion
1444 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1446 my $xas = $x->as_number();
1447 my $gs = $xas->copy()->bsqrt(); # some guess
1449 if (($x->{_e}->{sign} ne '-') # guess can't be accurate if there are
1450 # digits after the dot
1451 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
1454 $x->{_m} = $gs; $x->{_e} = $MBI->bzero(); $x->bnorm();
1455 # shortcut to not run trough _find_round_parameters again
1456 if (defined $params[0])
1458 $x->bround($params[0],$params[2]); # then round accordingly
1462 $x->bfround($params[1],$params[2]); # then round accordingly
1466 # clear a/p after round, since user did not request it
1467 $x->{_a} = undef; $x->{_p} = undef;
1469 # re-enable A and P, upgrade is taken care of by "local"
1470 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1474 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1475 # of the result by multipyling the input by 100 and then divide the integer
1476 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
1477 # this will transform 123.456 (in $x) into 123456 (in $y1)
1478 my $y1 = $x->{_m}->copy();
1479 # We now make sure that $y1 has the same odd or even number of digits than
1480 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1481 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1482 # steps of 10. The length of $x does not count, since an even or odd number
1483 # of digits before the dot is not changed by adding an even number of digits
1484 # after the dot (the result is still odd or even digits long).
1485 my $length = $y1->length();
1486 $y1->bmul(10) if $x->{_e}->is_odd();
1487 # now calculate how many digits the result of sqrt(y1) would have
1488 my $digits = int($length / 2);
1489 # but we need at least $scale digits, so calculate how many are missing
1490 my $shift = $scale - $digits;
1491 # that should never happen (we take care of integer guesses above)
1492 # $shift = 0 if $shift < 0;
1493 # multiply in steps of 100, by shifting left two times the "missing" digits
1494 $y1->blsft($shift*2,10);
1495 # now take the square root and truncate to integer
1497 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1498 # result, which is than later rounded to the desired scale.
1500 # calculate how many zeros $x had after the '.' (or before it, depending
1501 # on sign of $dat, the result should have half as many:
1502 my $dat = $length + $x->{_e}->numify();
1506 # no zeros after the dot (e.g. 1.23, 0.49 etc)
1507 # preserve half as many digits before the dot than the input had
1508 # (but round this "up")
1509 $dat = int(($dat+1)/2);
1513 $dat = int(($dat)/2);
1515 $x->{_e}= $MBI->new( $dat - $y1->length() );
1519 # shortcut to not run trough _find_round_parameters again
1520 if (defined $params[0])
1522 $x->bround($params[0],$params[2]); # then round accordingly
1526 $x->bfround($params[1],$params[2]); # then round accordingly
1530 # clear a/p after round, since user did not request it
1531 $x->{_a} = undef; $x->{_p} = undef;
1534 $$abr = $ab; $$pbr = $pb;
1540 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1541 # compute factorial numbers
1542 # modifies first argument
1543 my ($self,$x,@r) = objectify(1,@_);
1546 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1547 ($x->{_e}->{sign} ne '+')); # digits after dot?
1549 return $x->bone('+',@r) if $x->is_zero() || $x->is_one(); # 0 or 1 => 1
1551 # use BigInt's bfac() for faster calc
1552 $x->{_m}->blsft($x->{_e},10); # un-norm m
1553 $x->{_e}->bzero(); # norm $x again
1554 $x->{_m}->bfac(); # factorial
1555 $x->bnorm()->round(@r);
1560 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1561 my ($x,$y,$a,$p,$r) = @_;
1564 # if $y == 0.5, it is sqrt($x)
1565 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1568 # a ** x == e ** (x * ln a)
1572 # Taylor: | u u^2 u^3 |
1573 # x ** y = 1 + | --- + --- + ----- + ... |
1576 # we need to limit the accuracy to protect against overflow
1578 my ($scale,@params);
1579 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1581 return $x if $x->is_nan(); # error in _find_round_parameters?
1583 # no rounding at all, so must use fallback
1584 if (scalar @params == 0)
1586 # simulate old behaviour
1587 $params[0] = $self->div_scale(); # and round to it as accuracy
1588 $params[1] = undef; # disable P
1589 $scale = $params[0]+4; # at least four more for proper round
1590 $params[2] = $r; # round mode by caller or undef
1591 $fallback = 1; # to clear a/p afterwards
1595 # the 4 below is empirical, and there might be cases where it is not
1597 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1600 # when user set globals, they would interfere with our calculation, so
1601 # disable them and later re-enable them
1603 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1604 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1605 # we also need to disable any set A or P on $x (_find_round_parameters took
1606 # them already into account), since these would interfere, too
1607 delete $x->{_a}; delete $x->{_p};
1608 # need to disable $upgrade in BigInt, to avoid deep recursion
1609 local $Math::BigInt::upgrade = undef;
1611 my ($limit,$v,$u,$below,$factor,$next,$over);
1613 $u = $x->copy()->blog(undef,$scale)->bmul($y);
1614 $v = $self->bone(); # 1
1615 $factor = $self->new(2); # 2
1616 $x->bone(); # first term: 1
1618 $below = $v->copy();
1621 $limit = $self->new("1E-". ($scale-1));
1625 # we calculate the next term, and add it to the last
1626 # when the next term is below our limit, it won't affect the outcome
1627 # anymore, so we stop
1628 $next = $over->copy()->bdiv($below,$scale);
1629 last if $next->bacmp($limit) <= 0;
1631 # calculate things for the next term
1632 $over *= $u; $below *= $factor; $factor->binc();
1636 # shortcut to not run trough _find_round_parameters again
1637 if (defined $params[0])
1639 $x->bround($params[0],$params[2]); # then round accordingly
1643 $x->bfround($params[1],$params[2]); # then round accordingly
1647 # clear a/p after round, since user did not request it
1648 $x->{_a} = undef; $x->{_p} = undef;
1651 $$abr = $ab; $$pbr = $pb;
1657 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1658 # compute power of two numbers, second arg is used as integer
1659 # modifies first argument
1662 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1663 # objectify is costly, so avoid it
1664 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1666 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1669 return $x if $x->{sign} =~ /^[+-]inf$/;
1670 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1671 return $x->bone() if $y->is_zero();
1672 return $x if $x->is_one() || $y->is_one();
1674 return $x->_pow($y,$a,$p,$r) if !$y->is_int(); # non-integer power
1676 my $y1 = $y->as_number(); # make bigint
1678 if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1680 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
1681 return $y1->is_odd() ? $x : $x->babs(1);
1685 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1686 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1690 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1692 $x->{_m}->bpow($y1);
1693 $x->{_e}->bmul($y1);
1694 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1696 if ($y->{sign} eq '-')
1698 # modify $x in place!
1699 my $z = $x->copy(); $x->bzero()->binc();
1700 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
1702 $x->round($a,$p,$r,$y);
1705 ###############################################################################
1706 # rounding functions
1710 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1711 # $n == 0 means round to integer
1712 # expects and returns normalized numbers!
1713 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1715 return $x if $x->modify('bfround');
1717 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1718 return $x if !defined $scale; # no-op
1720 # never round a 0, +-inf, NaN
1723 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1726 return $x if $x->{sign} !~ /^[+-]$/;
1728 # don't round if x already has lower precision
1729 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1731 $x->{_p} = $scale; # remember round in any case
1732 $x->{_a} = undef; # and clear A
1735 # round right from the '.'
1737 return $x if $x->{_e}->{sign} eq '+'; # e >= 0 => nothing to round
1739 $scale = -$scale; # positive for simplicity
1740 my $len = $x->{_m}->length(); # length of mantissa
1742 # the following poses a restriction on _e, but if _e is bigger than a
1743 # scalar, you got other problems (memory etc) anyway
1744 my $dad = -($x->{_e}->numify()); # digits after dot
1745 my $zad = 0; # zeros after dot
1746 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
1748 #print "scale $scale dad $dad zad $zad len $len\n";
1749 # number bsstr len zad dad
1750 # 0.123 123e-3 3 0 3
1751 # 0.0123 123e-4 3 1 4
1754 # 1.2345 12345e-4 5 0 4
1756 # do not round after/right of the $dad
1757 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1759 # round to zero if rounding inside the $zad, but not for last zero like:
1760 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1761 return $x->bzero() if $scale < $zad;
1762 if ($scale == $zad) # for 0.006, scale -3 and trunc
1768 # adjust round-point to be inside mantissa
1771 $scale = $scale-$zad;
1775 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
1776 $scale = $dbd+$scale;
1782 # round left from the '.'
1784 # 123 => 100 means length(123) = 3 - $scale (2) => 1
1786 my $dbt = $x->{_m}->length();
1788 my $dbd = $dbt + $x->{_e}->numify();
1789 # should be the same, so treat it as this
1790 $scale = 1 if $scale == 0;
1791 # shortcut if already integer
1792 return $x if $scale == 1 && $dbt <= $dbd;
1793 # maximum digits before dot
1798 # not enough digits before dot, so round to zero
1801 elsif ( $scale == $dbd )
1808 $scale = $dbd - $scale;
1811 # pass sign to bround for rounding modes '+inf' and '-inf'
1812 $x->{_m}->{sign} = $x->{sign};
1813 $x->{_m}->bround($scale,$mode);
1814 $x->{_m}->{sign} = '+'; # fix sign back
1820 # accuracy: preserve $N digits, and overwrite the rest with 0's
1821 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1823 if (($_[0] || 0) < 0)
1825 require Carp; Carp::croak ('bround() needs positive accuracy');
1828 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1829 return $x if !defined $scale; # no-op
1831 return $x if $x->modify('bround');
1833 # scale is now either $x->{_a}, $accuracy, or the user parameter
1834 # test whether $x already has lower accuracy, do nothing in this case
1835 # but do round if the accuracy is the same, since a math operation might
1836 # want to round a number with A=5 to 5 digits afterwards again
1837 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1839 # scale < 0 makes no sense
1840 # never round a +-inf, NaN
1841 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
1843 # 1: $scale == 0 => keep all digits
1844 # 2: never round a 0
1845 # 3: if we should keep more digits than the mantissa has, do nothing
1846 if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1848 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1852 # pass sign to bround for '+inf' and '-inf' rounding modes
1853 $x->{_m}->{sign} = $x->{sign};
1854 $x->{_m}->bround($scale,$mode); # round mantissa
1855 $x->{_m}->{sign} = '+'; # fix sign back
1856 # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1857 $x->{_a} = $scale; # remember rounding
1858 $x->{_p} = undef; # and clear P
1859 $x->bnorm(); # del trailing zeros gen. by bround()
1864 # return integer less or equal then $x
1865 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1867 return $x if $x->modify('bfloor');
1869 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1871 # if $x has digits after dot
1872 if ($x->{_e}->{sign} eq '-')
1874 $x->{_e}->{sign} = '+'; # negate e
1875 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1876 $x->{_e}->bzero(); # trunc/norm
1877 $x->{_m}->binc() if $x->{sign} eq '-'; # decrement if negative
1879 $x->round($a,$p,$r);
1884 # return integer greater or equal then $x
1885 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1887 return $x if $x->modify('bceil');
1888 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1890 # if $x has digits after dot
1891 if ($x->{_e}->{sign} eq '-')
1893 #$x->{_m}->brsft(-$x->{_e},10);
1895 #$x++ if $x->{sign} eq '+';
1897 $x->{_e}->{sign} = '+'; # negate e
1898 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1899 $x->{_e}->bzero(); # trunc/norm
1900 $x->{_m}->binc() if $x->{sign} eq '+'; # decrement if negative
1902 $x->round($a,$p,$r);
1907 # shift right by $y (divide by power of $n)
1910 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1911 # objectify is costly, so avoid it
1912 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1914 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1917 return $x if $x->modify('brsft');
1918 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1920 $n = 2 if !defined $n; $n = $self->new($n);
1921 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
1926 # shift left by $y (multiply by power of $n)
1929 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1930 # objectify is costly, so avoid it
1931 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1933 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1936 return $x if $x->modify('blsft');
1937 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1939 $n = 2 if !defined $n; $n = $self->new($n);
1940 $x->bmul($n->bpow($y),$a,$p,$r,$y);
1943 ###############################################################################
1947 # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1952 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
1953 # or falling back to MBI::bxxx()
1954 my $name = $AUTOLOAD;
1956 $name =~ s/.*:://; # split package
1958 $class->import() if $IMPORT == 0;
1959 if (!method_alias($name))
1963 # delayed load of Carp and avoid recursion
1965 Carp::croak ("Can't call a method without name");
1967 if (!method_hand_up($name))
1969 # delayed load of Carp and avoid recursion
1971 Carp::croak ("Can't call $class\-\>$name, not a valid method");
1973 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1975 return &{"$MBI"."::$name"}(@_);
1977 my $bname = $name; $bname =~ s/^f/b/;
1978 *{$class."::$name"} = \&$bname;
1984 # return a copy of the exponent
1985 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1987 if ($x->{sign} !~ /^[+-]$/)
1989 my $s = $x->{sign}; $s =~ s/^[+-]//;
1990 return $self->new($s); # -inf, +inf => +inf
1992 return $x->{_e}->copy();
1997 # return a copy of the mantissa
1998 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2000 if ($x->{sign} !~ /^[+-]$/)
2002 my $s = $x->{sign}; $s =~ s/^[+]//;
2003 return $self->new($s); # -inf, +inf => +inf
2005 my $m = $x->{_m}->copy(); # faster than going via bstr()
2006 $m->bneg() if $x->{sign} eq '-';
2013 # return a copy of both the exponent and the mantissa
2014 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2016 if ($x->{sign} !~ /^[+-]$/)
2018 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2019 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2021 my $m = $x->{_m}->copy(); # faster than going via bstr()
2022 $m->bneg() if $x->{sign} eq '-';
2023 return ($m,$x->{_e}->copy());
2026 ##############################################################################
2027 # private stuff (internal use only)
2033 my $lib = ''; my @a;
2035 for ( my $i = 0; $i < $l ; $i++)
2037 if ( $_[$i] eq ':constant' )
2039 # this rest causes overlord er load to step in
2040 overload::constant float => sub { $self->new(shift); };
2042 elsif ($_[$i] eq 'upgrade')
2044 # this causes upgrading
2045 $upgrade = $_[$i+1]; # or undef to disable
2048 elsif ($_[$i] eq 'downgrade')
2050 # this causes downgrading
2051 $downgrade = $_[$i+1]; # or undef to disable
2054 elsif ($_[$i] eq 'lib')
2056 # alternative library
2057 $lib = $_[$i+1] || ''; # default Calc
2060 elsif ($_[$i] eq 'with')
2062 # alternative class for our private parts()
2063 $MBI = $_[$i+1] || 'Math::BigInt'; # default Math::BigInt
2072 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2073 my $mbilib = eval { Math::BigInt->config()->{lib} };
2074 if ((defined $mbilib) && ($MBI eq 'Math::BigInt'))
2076 # MBI already loaded
2077 $MBI->import('lib',"$lib,$mbilib", 'objectify');
2081 # MBI not loaded, or with ne "Math::BigInt"
2082 $lib .= ",$mbilib" if defined $mbilib;
2083 $lib =~ s/^,//; # don't leave empty
2084 # replacement library can handle lib statement, but also could ignore it
2087 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2088 # used in the same script, or eval inside import().
2089 my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
2090 my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
2092 $file = File::Spec->catfile (@parts, $file);
2093 eval { require "$file"; };
2094 $MBI->import( lib => $lib, 'objectify' );
2098 my $rc = "use $MBI lib => '$lib', 'objectify';";
2104 require Carp; Carp::croak ("Couldn't load $MBI: $! $@");
2107 # any non :constant stuff is handled by our parent, Exporter
2108 # even if @_ is empty, to give it a chance
2109 $self->SUPER::import(@a); # for subclasses
2110 $self->export_to_level(1,$self,@a); # need this, too
2115 # adjust m and e so that m is smallest possible
2116 # round number according to accuracy and precision settings
2117 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2119 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2121 # if (!$x->{_m}->is_odd())
2123 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
2126 $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros);
2128 # for something like 0Ey, set y to 1, and -0 => +0
2129 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
2131 # this is to prevent automatically rounding when MBI's globals are set
2132 $x->{_m}->{_f} = MB_NEVER_ROUND;
2133 $x->{_e}->{_f} = MB_NEVER_ROUND;
2134 # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
2135 $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
2136 $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
2137 $x; # MBI bnorm is no-op, so dont call it
2140 ##############################################################################
2144 # return number as hexadecimal string (only for integers defined)
2145 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2147 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2148 return '0x0' if $x->is_zero();
2150 return $nan if $x->{_e}->{sign} ne '+'; # how to do 1e-1 in hex!?
2152 my $z = $x->{_m}->copy();
2153 if (!$x->{_e}->is_zero()) # > 0
2155 $z->blsft($x->{_e},10);
2157 $z->{sign} = $x->{sign};
2163 # return number as binary digit string (only for integers defined)
2164 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2166 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2167 return '0b0' if $x->is_zero();
2169 return $nan if $x->{_e}->{sign} ne '+'; # how to do 1e-1 in hex!?
2171 my $z = $x->{_m}->copy();
2172 if (!$x->{_e}->is_zero()) # > 0
2174 $z->blsft($x->{_e},10);
2176 $z->{sign} = $x->{sign};
2182 # return copy as a bigint representation of this BigFloat number
2183 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2185 my $z = $x->{_m}->copy();
2186 if ($x->{_e}->{sign} eq '-') # < 0
2188 $x->{_e}->{sign} = '+'; # flip
2189 $z->brsft($x->{_e},10);
2190 $x->{_e}->{sign} = '-'; # flip back
2192 elsif (!$x->{_e}->is_zero()) # > 0
2194 $z->blsft($x->{_e},10);
2196 $z->{sign} = $x->{sign};
2203 my $class = ref($x) || $x;
2204 $x = $class->new(shift) unless ref($x);
2206 return 1 if $x->{_m}->is_zero();
2207 my $len = $x->{_m}->length();
2208 $len += $x->{_e} if $x->{_e}->sign() eq '+';
2211 my $t = $MBI->bzero();
2212 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
2223 Math::BigFloat - Arbitrary size floating point math package
2230 $x = Math::BigFloat->new($str); # defaults to 0
2231 $nan = Math::BigFloat->bnan(); # create a NotANumber
2232 $zero = Math::BigFloat->bzero(); # create a +0
2233 $inf = Math::BigFloat->binf(); # create a +inf
2234 $inf = Math::BigFloat->binf('-'); # create a -inf
2235 $one = Math::BigFloat->bone(); # create a +1
2236 $one = Math::BigFloat->bone('-'); # create a -1
2239 $x->is_zero(); # true if arg is +0
2240 $x->is_nan(); # true if arg is NaN
2241 $x->is_one(); # true if arg is +1
2242 $x->is_one('-'); # true if arg is -1
2243 $x->is_odd(); # true if odd, false for even
2244 $x->is_even(); # true if even, false for odd
2245 $x->is_positive(); # true if >= 0
2246 $x->is_negative(); # true if < 0
2247 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2249 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2250 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2251 $x->sign(); # return the sign, either +,- or NaN
2252 $x->digit($n); # return the nth digit, counting from right
2253 $x->digit(-$n); # return the nth digit, counting from left
2255 # The following all modify their first argument. If you want to preserve
2256 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2257 # neccessary when mixing $a = $b assigments with non-overloaded math.
2260 $x->bzero(); # set $i to 0
2261 $x->bnan(); # set $i to NaN
2262 $x->bone(); # set $x to +1
2263 $x->bone('-'); # set $x to -1
2264 $x->binf(); # set $x to inf
2265 $x->binf('-'); # set $x to -inf
2267 $x->bneg(); # negation
2268 $x->babs(); # absolute value
2269 $x->bnorm(); # normalize (no-op)
2270 $x->bnot(); # two's complement (bit wise not)
2271 $x->binc(); # increment x by 1
2272 $x->bdec(); # decrement x by 1
2274 $x->badd($y); # addition (add $y to $x)
2275 $x->bsub($y); # subtraction (subtract $y from $x)
2276 $x->bmul($y); # multiplication (multiply $x by $y)
2277 $x->bdiv($y); # divide, set $x to quotient
2278 # return (quo,rem) or quo if scalar
2280 $x->bmod($y); # modulus ($x % $y)
2281 $x->bpow($y); # power of arguments ($x ** $y)
2282 $x->blsft($y); # left shift
2283 $x->brsft($y); # right shift
2284 # return (quo,rem) or quo if scalar
2286 $x->blog(); # logarithm of $x to base e (Euler's number)
2287 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
2289 $x->band($y); # bit-wise and
2290 $x->bior($y); # bit-wise inclusive or
2291 $x->bxor($y); # bit-wise exclusive or
2292 $x->bnot(); # bit-wise not (two's complement)
2294 $x->bsqrt(); # calculate square-root
2295 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
2296 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2298 $x->bround($N); # accuracy: preserve $N digits
2299 $x->bfround($N); # precision: round to the $Nth digit
2301 $x->bfloor(); # return integer less or equal than $x
2302 $x->bceil(); # return integer greater or equal than $x
2304 # The following do not modify their arguments:
2306 bgcd(@values); # greatest common divisor
2307 blcm(@values); # lowest common multiplicator
2309 $x->bstr(); # return string
2310 $x->bsstr(); # return string in scientific notation
2312 $x->exponent(); # return exponent as BigInt
2313 $x->mantissa(); # return mantissa as BigInt
2314 $x->parts(); # return (mantissa,exponent) as BigInt
2316 $x->length(); # number of digits (w/o sign and '.')
2317 ($l,$f) = $x->length(); # number of digits, and length of fraction
2319 $x->precision(); # return P of $x (or global, if P of $x undef)
2320 $x->precision($n); # set P of $x to $n
2321 $x->accuracy(); # return A of $x (or global, if A of $x undef)
2322 $x->accuracy($n); # set A $x to $n
2324 # these get/set the appropriate global value for all BigFloat objects
2325 Math::BigFloat->precision(); # Precision
2326 Math::BigFloat->accuracy(); # Accuracy
2327 Math::BigFloat->round_mode(); # rounding mode
2331 All operators (inlcuding basic math operations) are overloaded if you
2332 declare your big floating point numbers as
2334 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2336 Operations with overloaded operators preserve the arguments, which is
2337 exactly what you expect.
2339 =head2 Canonical notation
2341 Input to these routines are either BigFloat objects, or strings of the
2342 following four forms:
2356 C</^[+-]\d+E[+-]?\d+$/>
2360 C</^[+-]\d*\.\d+E[+-]?\d+$/>
2364 all with optional leading and trailing zeros and/or spaces. Additonally,
2365 numbers are allowed to have an underscore between any two digits.
2367 Empty strings as well as other illegal numbers results in 'NaN'.
2369 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2370 are always stored in normalized form. On a string, it creates a BigFloat
2375 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2377 The string output will always have leading and trailing zeros stripped and drop
2378 a plus sign. C<bstr()> will give you always the form with a decimal point,
2379 while C<bsstr()> (s for scientific) gives you the scientific notation.
2381 Input bstr() bsstr()
2383 ' -123 123 123' '-123123123' '-123123123E0'
2384 '00.0123' '0.0123' '123E-4'
2385 '123.45E-2' '1.2345' '12345E-4'
2386 '10E+3' '10000' '1E4'
2388 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2389 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2390 return either undef, <0, 0 or >0 and are suited for sort.
2392 Actual math is done by using the class defined with C<with => Class;> (which
2393 defaults to BigInts) to represent the mantissa and exponent.
2395 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2396 represent the result when input arguments are not numbers, as well as
2397 the result of dividing by zero.
2399 =head2 C<mantissa()>, C<exponent()> and C<parts()>
2401 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2402 as BigInts such that:
2404 $m = $x->mantissa();
2405 $e = $x->exponent();
2406 $y = $m * ( 10 ** $e );
2407 print "ok\n" if $x == $y;
2409 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2411 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2413 Currently the mantissa is reduced as much as possible, favouring higher
2414 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2415 This might change in the future, so do not depend on it.
2417 =head2 Accuracy vs. Precision
2419 See also: L<Rounding|Rounding>.
2421 Math::BigFloat supports both precision and accuracy. For a full documentation,
2422 examples and tips on these topics please see the large section in
2425 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2426 a operation consumes all resources, each operation produces no more than
2427 the requested number of digits.
2429 Please refer to BigInt's documentation for the precedence rules of which
2430 accuracy/precision setting will be used.
2432 If there is no gloabl precision set, B<and> the operation inquestion was not
2433 called with a requested precision or accuracy, B<and> the input $x has no
2434 accuracy or precision set, then a fallback parameter will be used. For
2435 historical reasons, it is called C<div_scale> and can be accessed via:
2437 $d = Math::BigFloat->div_scale(); # query
2438 Math::BigFloat->div_scale($n); # set to $n digits
2440 The default value is 40 digits.
2442 In case the result of one operation has more precision than specified,
2443 it is rounded. The rounding mode taken is either the default mode, or the one
2444 supplied to the operation after the I<scale>:
2446 $x = Math::BigFloat->new(2);
2447 Math::BigFloat->precision(5); # 5 digits max
2448 $y = $x->copy()->bdiv(3); # will give 0.66666
2449 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2450 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
2451 Math::BigFloat->round_mode('zero');
2452 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2458 =item ffround ( +$scale )
2460 Rounds to the $scale'th place left from the '.', counting from the dot.
2461 The first digit is numbered 1.
2463 =item ffround ( -$scale )
2465 Rounds to the $scale'th place right from the '.', counting from the dot.
2469 Rounds to an integer.
2471 =item fround ( +$scale )
2473 Preserves accuracy to $scale digits from the left (aka significant digits)
2474 and pads the rest with zeros. If the number is between 1 and -1, the
2475 significant digits count from the first non-zero after the '.'
2477 =item fround ( -$scale ) and fround ( 0 )
2479 These are effectively no-ops.
2483 All rounding functions take as a second parameter a rounding mode from one of
2484 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2486 The default rounding mode is 'even'. By using
2487 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
2488 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2489 no longer supported.
2490 The second parameter to the round functions then overrides the default
2493 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
2494 'trunc' as rounding mode to make it equivalent to:
2499 You can override this by passing the desired rounding mode as parameter to
2502 $x = Math::BigFloat->new(2.5);
2503 $y = $x->as_number('odd'); # $y = 3
2509 =head1 Autocreating constants
2511 After C<use Math::BigFloat ':constant'> all the floating point constants
2512 in the given scope are converted to C<Math::BigFloat>. This conversion
2513 happens at compile time.
2517 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2519 prints the value of C<2E-100>. Note that without conversion of
2520 constants the expression 2E-100 will be calculated as normal floating point
2523 Please note that ':constant' does not affect integer constants, nor binary
2524 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2529 Math with the numbers is done (by default) by a module called
2530 Math::BigInt::Calc. This is equivalent to saying:
2532 use Math::BigFloat lib => 'Calc';
2534 You can change this by using:
2536 use Math::BigFloat lib => 'BitVect';
2538 The following would first try to find Math::BigInt::Foo, then
2539 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2541 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2543 Calc.pm uses as internal format an array of elements of some decimal base
2544 (usually 1e7, but this might be differen for some systems) with the least
2545 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2546 significant bit first. Other modules might use even different means of
2547 representing the numbers. See the respective module documentation for further
2550 Please note that Math::BigFloat does B<not> use the denoted library itself,
2551 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2554 use Math::BigInt lib => 'GMP';
2557 you can roll it all into one line:
2559 use Math::BigFloat lib => 'GMP';
2561 It is also possible to just require Math::BigFloat:
2563 require Math::BigFloat;
2565 This will load the neccessary things (like BigInt) when they are needed, and
2568 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2569 you ever wanted to know about loading a different library.
2571 =head2 Using Math::BigInt::Lite
2573 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2576 use Math::BigFloat with => 'Math::BigInt::Lite';
2578 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2579 can combine these if you want. For instance, you may want to use
2580 Math::BigInt objects in your main script, too.
2584 use Math::BigFloat with => 'Math::BigInt::Lite';
2586 Of course, you can combine this with the C<lib> parameter.
2589 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2591 There is no need for a "use Math::BigInt;" statement, even if you want to
2592 use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2593 always loads it. But if you add it, add it B<before>:
2597 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2599 Notice that the module with the last C<lib> will "win" and thus
2600 it's lib will be used if the lib is available:
2603 use Math::BigInt lib => 'Bar,Baz';
2604 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2606 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2607 words, Math::BigFloat will try to retain previously loaded libs when you
2608 don't specify it onem but if you specify one, it will try to load them.
2610 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2611 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2612 same as trying the latter load alone, except for the fact that one of Bar or
2613 Baz might be loaded needlessly in an intermidiate step (and thus hang around
2614 and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2615 will still be tried to be loaded, but this is not as time/memory consuming as
2616 actually loading one of them. Still, this type of usage is not recommended due
2619 The old way (loading the lib only in BigInt) still works though:
2622 use Math::BigInt lib => 'Bar,Baz';
2625 You can even load Math::BigInt afterwards:
2629 use Math::BigInt lib => 'Bar,Baz';
2631 But this has the same problems like #5, it will first load Calc
2632 (Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2633 Baz, depending on which of them works and is usable/loadable. Since this
2634 loads Calc unnecc., it is not recommended.
2636 Since it also possible to just require Math::BigFloat, this poses the question
2637 about what libary this will use:
2639 require Math::BigFloat;
2640 my $x = Math::BigFloat->new(123); $x += 123;
2642 It will use Calc. Please note that the call to import() is still done, but
2643 only when you use for the first time some Math::BigFloat math (it is triggered
2644 via any constructor, so the first time you create a Math::BigFloat, the load
2645 will happen in the background). This means:
2647 require Math::BigFloat;
2648 Math::BigFloat->import ( lib => 'Foo,Bar' );
2650 would be the same as:
2652 use Math::BigFloat lib => 'Foo, Bar';
2654 But don't try to be clever to insert some operations in between:
2656 require Math::BigFloat;
2657 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
2658 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
2659 $x = Math::BigFloat->bone()+4; # now use Pari
2661 While this works, it loads Calc needlessly. But maybe you just wanted that?
2663 B<Examples #3 is highly recommended> for daily usage.
2667 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
2673 =item stringify, bstr()
2675 Both stringify and bstr() now drop the leading '+'. The old code would return
2676 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2677 reasoning and details.
2681 The following will probably not do what you expect:
2683 print $c->bdiv(123.456),"\n";
2685 It prints both quotient and reminder since print works in list context. Also,
2686 bdiv() will modify $c, so be carefull. You probably want to use
2688 print $c / 123.456,"\n";
2689 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2693 =item Modifying and =
2697 $x = Math::BigFloat->new(5);
2700 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2701 a second reference to the B<same> object and stores it in $y. Thus anything
2702 that modifies $x will modify $y (except overloaded math operators), and vice
2703 versa. See L<Math::BigInt> for details and how to avoid that.
2707 C<bpow()> now modifies the first argument, unlike the old code which left
2708 it alone and only returned the result. This is to be consistent with
2709 C<badd()> etc. The first will modify $x, the second one won't:
2711 print bpow($x,$i),"\n"; # modify $x
2712 print $x->bpow($i),"\n"; # ditto
2713 print $x ** $i,"\n"; # leave $x alone
2719 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
2720 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
2722 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
2723 because they solve the autoupgrading/downgrading issue, at least partly.
2726 L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
2727 more documentation including a full version history, testcases, empty
2728 subclass files and benchmarks.
2732 This program is free software; you may redistribute it and/or modify it under
2733 the same terms as Perl itself.
2737 Mark Biggar, overloaded interface by Ilya Zakharevich.
2738 Completely rewritten by Tels http://bloodgate.com in 2001, 2002, and still