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: +,-,"NaN" if not a number
13 # _f: flags, used to signal MBI not to touch our private parts
20 @ISA = qw( Exporter Math::BigInt);
23 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode/;
24 use vars qw/$upgrade $downgrade/;
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 accessory
37 use constant MB_NEVER_ROUND => 0x0001;
41 # constant for easier life
44 # class constants, use Class->constant_name() to access
45 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
52 my $MBI = 'Math::BigInt'; # the package we are using for our private parts
53 # changable by use Math::BigFloat with => 'package'
55 ##############################################################################
56 # the old code had $rnd_mode, so we need to support it, too
58 sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
59 sub FETCH { return $round_mode; }
60 sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
65 tie $rnd_mode, 'Math::BigFloat';
68 ##############################################################################
70 # in case we call SUPER::->foo() and this wants to call modify()
71 # sub modify () { 0; }
74 # valid method aliases for AUTOLOAD
75 my %methods = map { $_ => 1 }
76 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
77 fint facmp fcmp fzero fnan finf finc fdec flog ffac
78 fceil ffloor frsft flsft fone flog
80 # valid method's that can be hand-ed up (for AUTOLOAD)
81 my %hand_ups = map { $_ => 1 }
82 qw / is_nan is_inf is_negative is_positive
83 accuracy precision div_scale round_mode fneg fabs babs fnot
84 objectify upgrade downgrade
88 sub method_alias { return exists $methods{$_[0]||''}; }
89 sub method_hand_up { return exists $hand_ups{$_[0]||''}; }
92 ##############################################################################
97 # create a new BigFloat object from a string or another bigfloat object.
100 # sign => sign (+/-), or "NaN"
102 my ($class,$wanted,@r) = @_;
104 # avoid numify-calls by not using || on $wanted!
105 return $class->bzero() if !defined $wanted; # default to 0
106 return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
108 my $self = {}; bless $self, $class;
109 # shortcut for bigints and its subclasses
110 if ((ref($wanted)) && (ref($wanted) ne $class))
112 $self->{_m} = $wanted->as_number(); # get us a bigint copy
113 $self->{_e} = $MBI->bzero();
115 $self->{sign} = $wanted->sign();
116 return $self->bnorm();
119 # handle '+inf', '-inf' first
120 if ($wanted =~ /^[+-]?inf$/)
122 return $downgrade->new($wanted) if $downgrade;
124 $self->{_e} = $MBI->bzero();
125 $self->{_m} = $MBI->bzero();
126 $self->{sign} = $wanted;
127 $self->{sign} = '+inf' if $self->{sign} eq 'inf';
128 return $self->bnorm();
130 #print "new string '$wanted'\n";
131 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
134 die "$wanted is not a number initialized to $class" if !$NaNOK;
136 return $downgrade->bnan() if $downgrade;
138 $self->{_e} = $MBI->bzero();
139 $self->{_m} = $MBI->bzero();
140 $self->{sign} = $nan;
144 # make integer from mantissa by adjusting exp, then convert to bigint
145 # undef,undef to signal MBI that we don't need no bloody rounding
146 $self->{_e} = $MBI->new("$$es$$ev",undef,undef); # exponent
147 $self->{_m} = $MBI->new("$$miv$$mfv",undef,undef); # create mant.
148 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
149 $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;
150 $self->{sign} = $$mis;
152 # if downgrade, inf, NaN or integers go down
154 if ($downgrade && $self->{_e}->{sign} eq '+')
156 # print "downgrading $$miv$$mfv"."E$$es$$ev";
157 if ($self->{_e}->is_zero())
159 $self->{_m}->{sign} = $$mis; # negative if wanted
160 return $downgrade->new($self->{_m});
162 return $downgrade->new("$$mis$$miv$$mfv"."E$$es$$ev");
164 # print "mbf new $self->{sign} $self->{_m} e $self->{_e} ",ref($self),"\n";
165 $self->bnorm()->round(@r); # first normalize, then round
170 # used by parent class bone() to initialize number to 1
172 $self->{_m} = $MBI->bzero();
173 $self->{_e} = $MBI->bzero();
178 # used by parent class bone() to initialize number to 1
180 $self->{_m} = $MBI->bzero();
181 $self->{_e} = $MBI->bzero();
186 # used by parent class bone() to initialize number to 1
188 $self->{_m} = $MBI->bone();
189 $self->{_e} = $MBI->bzero();
194 # used by parent class bone() to initialize number to 1
196 $self->{_m} = $MBI->bzero();
197 $self->{_e} = $MBI->bone();
202 my ($self,$class) = @_;
203 return if $class =~ /^Math::BigInt/; # we aren't one of these
204 UNIVERSAL::isa($self,$class);
209 # return (later set?) configuration data as hash ref
210 my $class = shift || 'Math::BigFloat';
212 my $cfg = $MBI->config();
215 $cfg->{class} = $class;
218 qw/upgrade downgrade precision accuracy round_mode VERSION div_scale/)
220 $cfg->{lc($_)} = ${"${class}::$_"};
225 ##############################################################################
226 # string conversation
230 # (ref to BFLOAT or num_str ) return num_str
231 # Convert number from internal format to (non-scientific) string format.
232 # internal format is always normalized (no leading zeros, "-0" => "+0")
233 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
234 #my $x = shift; my $class = ref($x) || $x;
235 #$x = $class->new(shift) unless ref($x);
237 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
238 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
239 if ($x->{sign} !~ /^[+-]$/)
241 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
245 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
247 my $not_zero = ! $x->is_zero();
250 $es = $x->{_m}->bstr();
251 $len = CORE::length($es);
252 if (!$x->{_e}->is_zero())
254 if ($x->{_e}->sign() eq '-')
257 if ($x->{_e} <= -$len)
259 # print "style: 0.xxxx\n";
260 my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
261 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
265 # print "insert '.' at $x->{_e} in '$es'\n";
266 substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
272 $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
276 $es = $x->{sign}.$es if $x->{sign} eq '-';
277 # if set accuracy or precision, pad with zeros
278 if ((defined $x->{_a}) && ($not_zero))
280 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
281 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
282 $zeros = $x->{_a} - $len if $cad != $len;
283 $es .= $dot.'0' x $zeros if $zeros > 0;
285 elsif ($x->{_p} || 0 < 0)
287 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
288 my $zeros = -$x->{_p} + $cad;
289 $es .= $dot.'0' x $zeros if $zeros > 0;
296 # (ref to BFLOAT or num_str ) return num_str
297 # Convert number from internal format to scientific string format.
298 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
299 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
300 #my $x = shift; my $class = ref($x) || $x;
301 #$x = $class->new(shift) unless ref($x);
303 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
304 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
305 if ($x->{sign} !~ /^[+-]$/)
307 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
310 my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
312 $x->{_m}->bstr().$sep.$x->{_e}->bstr();
317 # Make a number from a BigFloat object
318 # simple return string and let Perl's atoi()/atof() handle the rest
319 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
323 ##############################################################################
324 # public stuff (usually prefixed with "b")
327 # todo: this must be overwritten and return NaN for non-integer values
328 # band(), bior(), bxor(), too
331 # $class->SUPER::bnot($class,@_);
336 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
337 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
338 my ($self,$x,$y) = objectify(2,@_);
340 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
342 # handle +-inf and NaN
343 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
344 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
345 return +1 if $x->{sign} eq '+inf';
346 return -1 if $x->{sign} eq '-inf';
347 return -1 if $y->{sign} eq '+inf';
351 # check sign for speed first
352 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
353 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
356 my $xz = $x->is_zero();
357 my $yz = $y->is_zero();
358 return 0 if $xz && $yz; # 0 <=> 0
359 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
360 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
362 # adjust so that exponents are equal
363 my $lxm = $x->{_m}->length();
364 my $lym = $y->{_m}->length();
365 # the numify somewhat limits our length, but makes it much faster
366 my $lx = $lxm + $x->{_e}->numify();
367 my $ly = $lym + $y->{_e}->numify();
368 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
369 return $l <=> 0 if $l != 0;
371 # lengths (corrected by exponent) are equal
372 # so make mantissa equal length by padding with zero (shift left)
373 my $diff = $lxm - $lym;
374 my $xm = $x->{_m}; # not yet copy it
378 $ym = $y->{_m}->copy()->blsft($diff,10);
382 $xm = $x->{_m}->copy()->blsft(-$diff,10);
384 my $rc = $xm->bacmp($ym);
385 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
391 # Compares 2 values, ignoring their signs.
392 # Returns one of undef, <0, =0, >0. (suitable for sort)
393 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
394 my ($self,$x,$y) = objectify(2,@_);
396 # handle +-inf and NaN's
397 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
399 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
400 return 0 if ($x->is_inf() && $y->is_inf());
401 return 1 if ($x->is_inf() && !$y->is_inf());
406 my $xz = $x->is_zero();
407 my $yz = $y->is_zero();
408 return 0 if $xz && $yz; # 0 <=> 0
409 return -1 if $xz && !$yz; # 0 <=> +y
410 return 1 if $yz && !$xz; # +x <=> 0
412 # adjust so that exponents are equal
413 my $lxm = $x->{_m}->length();
414 my $lym = $y->{_m}->length();
415 # the numify somewhat limits our length, but makes it much faster
416 my $lx = $lxm + $x->{_e}->numify();
417 my $ly = $lym + $y->{_e}->numify();
419 return $l <=> 0 if $l != 0;
421 # lengths (corrected by exponent) are equal
422 # so make mantissa equal-length by padding with zero (shift left)
423 my $diff = $lxm - $lym;
424 my $xm = $x->{_m}; # not yet copy it
428 $ym = $y->{_m}->copy()->blsft($diff,10);
432 $xm = $x->{_m}->copy()->blsft(-$diff,10);
434 $xm->bacmp($ym) <=> 0;
439 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
440 # return result as BFLOAT
441 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
443 #print "mbf badd $x $y\n";
444 # inf and NaN handling
445 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
448 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
450 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
452 # +inf++inf or -inf+-inf => same, rest is NaN
453 return $x if $x->{sign} eq $y->{sign};
456 # +-inf + something => +inf; something +-inf => +-inf
457 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
461 return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
462 ((!$x->isa($self)) || (!$y->isa($self)));
464 # speed: no add for 0+y or x+0
465 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
466 if ($x->is_zero()) # 0+y
468 # make copy, clobbering up x (modify in place!)
469 $x->{_e} = $y->{_e}->copy();
470 $x->{_m} = $y->{_m}->copy();
471 $x->{sign} = $y->{sign} || $nan;
472 return $x->round($a,$p,$r,$y);
475 # take lower of the two e's and adapt m1 to it to match m2
477 $e = $MBI->bzero() if !defined $e; # if no BFLOAT ?
478 $e = $e->copy(); # make copy (didn't do it yet)
480 my $add = $y->{_m}->copy();
481 if ($e->{sign} eq '-') # < 0
483 my $e1 = $e->copy()->babs();
484 #$x->{_m} *= (10 ** $e1);
485 $x->{_m}->blsft($e1,10);
486 $x->{_e} += $e; # need the sign of e
488 elsif (!$e->is_zero()) # > 0
493 # else: both e are the same, so just leave them
494 $x->{_m}->{sign} = $x->{sign}; # fiddle with signs
495 $add->{sign} = $y->{sign};
496 $x->{_m} += $add; # finally do add/sub
497 $x->{sign} = $x->{_m}->{sign}; # re-adjust signs
498 $x->{_m}->{sign} = '+'; # mantissa always positiv
499 # delete trailing zeros, then round
500 return $x->bnorm()->round($a,$p,$r,$y);
505 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
506 # subtract second arg from first, modify first
507 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
509 if ($y->is_zero()) # still round for not adding zero
511 return $x->round($a,$p,$r);
514 $y->{sign} =~ tr/+\-/-+/; # does nothing for NaN
515 $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
516 $y->{sign} =~ tr/+\-/-+/; # refix $y (does nothing for NaN)
517 $x; # already rounded by badd()
522 # increment arg by one
523 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
525 if ($x->{_e}->sign() eq '-')
527 return $x->badd($self->bone(),$a,$p,$r); # digits after dot
530 if (!$x->{_e}->is_zero())
532 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
536 if ($x->{sign} eq '+')
539 return $x->bnorm()->bround($a,$p,$r);
541 elsif ($x->{sign} eq '-')
544 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
545 return $x->bnorm()->bround($a,$p,$r);
547 # inf, nan handling etc
548 $x->badd($self->__one(),$a,$p,$r); # does round
553 # decrement arg by one
554 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
556 if ($x->{_e}->sign() eq '-')
558 return $x->badd($self->bone('-'),$a,$p,$r); # digits after dot
561 if (!$x->{_e}->is_zero())
563 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
567 my $zero = $x->is_zero();
569 if (($x->{sign} eq '-') || $zero)
572 $x->{sign} = '-' if $zero; # 0 => 1 => -1
573 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
574 return $x->bnorm()->round($a,$p,$r);
577 elsif ($x->{sign} eq '+')
580 return $x->bnorm()->round($a,$p,$r);
582 # inf, nan handling etc
583 $x->badd($self->bone('-'),$a,$p,$r); # does round
588 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(2,@_);
590 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
594 # Taylor: | u 1 u^3 1 u^5 |
595 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
596 # |_ v 3 v^3 5 v^5 _|
598 # This takes much more steps to calculate the result:
601 # Taylor: | u 1 u^2 1 u^3 |
602 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
603 # |_ x 2 x^2 3 x^3 _|
605 # we need to limit the accuracy to protect against overflow
608 my @params = $x->_find_round_parameters($a,$p,$r);
610 # no rounding at all, so must use fallback
611 if (scalar @params == 1)
613 # simulate old behaviour
614 $params[1] = $self->div_scale(); # and round to it as accuracy
615 $scale = $params[1]+4; # at least four more for proper round
616 $params[3] = $r; # round mode by caller or undef
617 $fallback = 1; # to clear a/p afterwards
621 # the 4 below is empirical, and there might be cases where it is not
623 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
626 return $x->bzero(@params) if $x->is_one();
627 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
628 #return $x->bone('+',@params) if $x->bcmp($base) == 0;
630 # when user set globals, they would interfere with our calculation, so
631 # disable then and later re-enable them
633 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
634 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
635 # we also need to disable any set A or P on $x (_find_round_parameters took
636 # them already into account), since these would interfere, too
637 delete $x->{_a}; delete $x->{_p};
638 # need to disable $upgrade in BigInt, to avoid deep recursion
639 local $Math::BigInt::upgrade = undef;
641 my ($case,$limit,$v,$u,$below,$factor,$two,$next,$over,$f);
644 #if ($x <= Math::BigFloat->new("0.5"))
647 # print "case $case $x < 0.5\n";
648 $v = $x->copy(); $v->binc(); # v = x+1
649 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
650 $x->bdiv($v,$scale); # first term: u/v
653 $u *= $u; $v *= $v; # u^2, v^2
654 $below->bmul($v); # u^3, v^3
656 $factor = $self->new(3); $f = $self->new(2);
661 # print "case 1 $x > 0.5\n";
662 # $v = $x->copy(); # v = x
663 # $u = $x->copy(); $u->bdec(); # u = x-1;
664 # $x->bdec(); $x->bdiv($v,$scale); # first term: x-1/x
665 # $below = $v->copy();
666 # $over = $u->copy();
667 # $below->bmul($v); # u^2, v^2
669 # $factor = $self->new(2); $f = $self->bone();
671 $limit = $self->new("1E-". ($scale-1));
675 # we calculate the next term, and add it to the last
676 # when the next term is below our limit, it won't affect the outcome
677 # anymore, so we stop
678 $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
679 last if $next->bcmp($limit) <= 0;
681 # print "step $steps $x\n";
682 # calculate things for the next term
683 $over *= $u; $below *= $v; $factor->badd($f);
686 $x->bmul(2) if $case == 0;
687 #print "took $steps steps\n";
689 # shortcut to not run trough _find_round_parameters again
690 if (defined $params[1])
692 $x->bround($params[1],$params[3]); # then round accordingly
696 $x->bfround($params[2],$params[3]); # then round accordingly
700 # clear a/p after round, since user did not request it
701 $x->{_a} = undef; $x->{_p} = undef;
704 $$abr = $ab; $$pbr = $pb;
711 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
712 # does not modify arguments, but returns new object
713 # Lowest Common Multiplicator
715 my ($self,@arg) = objectify(0,@_);
716 my $x = $self->new(shift @arg);
717 while (@arg) { $x = _lcm($x,shift @arg); }
723 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
724 # does not modify arguments, but returns new object
725 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
727 my ($self,@arg) = objectify(0,@_);
728 my $x = $self->new(shift @arg);
729 while (@arg) { $x = _gcd($x,shift @arg); }
733 ###############################################################################
734 # is_foo methods (is_negative, is_positive are inherited from BigInt)
738 # return true if arg (BFLOAT or num_str) is an integer
739 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
741 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
742 $x->{_e}->{sign} eq '+'; # 1e-1 => no integer
748 # return true if arg (BFLOAT or num_str) is zero
749 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
751 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
757 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
758 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
760 my $sign = shift || ''; $sign = '+' if $sign ne '-';
762 if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
768 # return true if arg (BFLOAT or num_str) is odd or false if even
769 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
771 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
772 ($x->{_e}->is_zero() && $x->{_m}->is_odd());
778 # return true if arg (BINT or num_str) is even or false if odd
779 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
781 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
782 return 1 if ($x->{_e}->{sign} eq '+' # 123.45 is never
783 && $x->{_m}->is_even()); # but 1200 is
789 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
790 # (BINT or num_str, BINT or num_str) return BINT
791 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
793 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
796 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
798 return $x->bnan() if $x->is_zero() || $y->is_zero();
799 # result will always be +-inf:
800 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
801 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
802 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
803 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
804 return $x->binf('-');
807 return $x->bzero() if $x->is_zero() || $y->is_zero();
809 return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
810 ((!$x->isa($self)) || (!$y->isa($self)));
812 # aEb * cEd = (a*c)E(b+d)
813 $x->{_m}->bmul($y->{_m});
814 $x->{_e}->badd($y->{_e});
816 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
817 return $x->bnorm()->round($a,$p,$r,$y);
822 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
823 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
824 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
826 return $self->_div_inf($x,$y)
827 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
829 # x== 0 # also: or y == 1 or y == -1
830 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
833 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
835 # we need to limit the accuracy to protect against overflow
838 my @params = $x->_find_round_parameters($a,$p,$r,$y);
840 # no rounding at all, so must use fallback
841 if (scalar @params == 1)
843 # simulate old behaviour
844 $params[1] = $self->div_scale(); # and round to it as accuracy
845 $scale = $params[1]+4; # at least four more for proper round
846 $params[3] = $r; # round mode by caller or undef
847 $fallback = 1; # to clear a/p afterwards
851 # the 4 below is empirical, and there might be cases where it is not
853 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
855 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
856 $scale = $lx if $lx > $scale;
857 $scale = $ly if $ly > $scale;
858 my $diff = $ly - $lx;
859 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
861 # make copy of $x in case of list context for later reminder calculation
863 if (wantarray && !$y->is_one())
868 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
870 # check for / +-1 ( +/- 1E0)
873 # promote BigInts and it's subclasses (except when already a BigFloat)
874 $y = $self->new($y) unless $y->isa('Math::BigFloat');
876 #print "bdiv $y ",ref($y),"\n";
877 # need to disable $upgrade in BigInt, to avoid deep recursion
878 local $Math::BigInt::upgrade = undef; # should be parent class vs MBI
880 # calculate the result to $scale digits and then round it
881 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
882 $x->{_m}->blsft($scale,10);
883 $x->{_m}->bdiv( $y->{_m} ); # a/c
884 $x->{_e}->bsub( $y->{_e} ); # b-d
885 $x->{_e}->bsub($scale); # correct for 10**scale
886 $x->bnorm(); # remove trailing 0's
889 # shortcut to not run trough _find_round_parameters again
890 if (defined $params[1])
892 $x->bround($params[1],$params[3]); # then round accordingly
896 $x->bfround($params[2],$params[3]); # then round accordingly
900 # clear a/p after round, since user did not request it
901 $x->{_a} = undef; $x->{_p} = undef;
908 $rem->bmod($y,$params[1],$params[2],$params[3]); # copy already done
912 $rem = $self->bzero();
916 # clear a/p after round, since user did not request it
917 $rem->{_a} = undef; $rem->{_p} = undef;
926 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
927 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
929 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
931 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
932 return $re->round($a,$p,$r,$y);
934 return $x->bnan() if $x->is_zero() && $y->is_zero();
935 return $x if $y->is_zero();
936 return $x->bnan() if $x->is_nan() || $y->is_nan();
937 return $x->bzero() if $y->is_one() || $x->is_zero();
939 # inf handling is missing here
941 my $cmp = $x->bacmp($y); # equal or $x < $y?
942 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
944 # only $y of the operands negative?
945 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
947 $x->{sign} = $y->{sign}; # calc sign first
948 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
950 my $ym = $y->{_m}->copy();
953 $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
955 # if $y has digits after dot
956 my $shifty = 0; # correct _e of $x by this
957 if ($y->{_e}->{sign} eq '-') # has digits after dot
959 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
960 $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
961 $x->blsft($shifty,10); # 123 => 1230, $y->{_m} is already 25
963 # $ym is now mantissa of $y based on exponent 0
965 my $shiftx = 0; # correct _e of $x by this
966 if ($x->{_e}->{sign} eq '-') # has digits after dot
968 # 123.4 % 20 => 1234 % 200
969 $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
970 $ym->blsft($shiftx,10);
972 # 123e1 % 20 => 1230 % 20
973 if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
975 $x->{_m}->blsft($x->{_e},10);
977 $x->{_e} = $MBI->bzero() unless $x->{_e}->is_zero();
979 $x->{_e}->bsub($shiftx) if $shiftx != 0;
980 $x->{_e}->bsub($shifty) if $shifty != 0;
982 # now mantissas are equalized, exponent of $x is adjusted, so calc result
986 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
989 if ($neg != 0) # one of them negative => correct in place
994 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
998 $x->round($a,$p,$r,$y); # round and return
1003 # calculate square root; this should probably
1004 # use a different test to see whether the accuracy we want is...
1005 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1007 return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
1008 return $x if $x->{sign} eq '+inf'; # +inf
1009 return $x if $x->is_zero() || $x->is_one();
1011 # we need to limit the accuracy to protect against overflow
1014 my @params = $x->_find_round_parameters($a,$p,$r);
1016 # no rounding at all, so must use fallback
1017 if (scalar @params == 1)
1019 # simulate old behaviour
1020 $params[1] = $self->div_scale(); # and round to it as accuracy
1021 $scale = $params[1]+4; # at least four more for proper round
1022 $params[3] = $r; # round mode by caller or undef
1023 $fallback = 1; # to clear a/p afterwards
1027 # the 4 below is empirical, and there might be cases where it is not
1029 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1032 # when user set globals, they would interfere with our calculation, so
1033 # disable them and later re-enable them
1035 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1036 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1037 # we also need to disable any set A or P on $x (_find_round_parameters took
1038 # them already into account), since these would interfere, too
1039 delete $x->{_a}; delete $x->{_p};
1040 # need to disable $upgrade in BigInt, to avoid deep recursion
1041 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1043 my $xas = $x->as_number();
1044 my $gs = $xas->copy()->bsqrt(); # some guess
1046 # print "guess $gs\n";
1047 if (($x->{_e}->{sign} ne '-') # guess can't be accurate if there are
1048 # digits after the dot
1049 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
1052 $x->{_m} = $gs; $x->{_e} = $MBI->bzero(); $x->bnorm();
1053 # shortcut to not run trough _find_round_parameters again
1054 if (defined $params[1])
1056 $x->bround($params[1],$params[3]); # then round accordingly
1060 $x->bfround($params[2],$params[3]); # then round accordingly
1064 # clear a/p after round, since user did not request it
1065 $x->{_a} = undef; $x->{_p} = undef;
1067 # re-enable A and P, upgrade is taken care of by "local"
1068 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1071 $gs = $self->new( $gs ); # BigInt to BigFloat
1073 my $lx = $x->{_m}->length();
1074 $scale = $lx if $scale < $lx;
1075 my $e = $self->new("1E-$scale"); # make test variable
1078 my $two = $self->new(2);
1080 # promote BigInts and it's subclasses (except when already a BigFloat)
1081 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1084 while ($diff->bacmp($e) >= 0)
1086 $rem = $y->copy()->bdiv($gs,$scale);
1087 $rem = $y->copy()->bdiv($gs,$scale)->badd($gs)->bdiv($two,$scale);
1088 $diff = $rem->copy()->bsub($gs);
1091 # copy over to modify $x
1092 $x->{_m} = $rem->{_m}; $x->{_e} = $rem->{_e};
1094 # shortcut to not run trough _find_round_parameters again
1095 if (defined $params[1])
1097 $x->bround($params[1],$params[3]); # then round accordingly
1101 $x->bfround($params[2],$params[3]); # then round accordingly
1105 # clear a/p after round, since user did not request it
1106 $x->{_a} = undef; $x->{_p} = undef;
1109 $$abr = $ab; $$pbr = $pb;
1115 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1116 # compute factorial numbers
1117 # modifies first argument
1118 my ($self,$x,@r) = objectify(1,@_);
1121 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1122 ($x->{_e}->{sign} ne '+')); # digits after dot?
1124 return $x->bone(@r) if $x->is_zero() || $x->is_one(); # 0 or 1 => 1
1126 # use BigInt's bfac() for faster calc
1127 $x->{_m}->blsft($x->{_e},10); # un-norm m
1128 $x->{_e}->bzero(); # norm $x again
1129 $x->{_m}->bfac(); # factorial
1130 $x->bnorm()->round(@r);
1135 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1136 my ($x,$y,$a,$p,$r) = @_;
1139 # we need to limit the accuracy to protect against overflow
1142 my @params = $x->_find_round_parameters($a,$p,$r);
1144 # no rounding at all, so must use fallback
1145 if (scalar @params == 1)
1147 # simulate old behaviour
1148 $params[1] = $self->div_scale(); # and round to it as accuracy
1149 $scale = $params[1]+4; # at least four more for proper round
1150 $params[3] = $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[1] || $params[2]) + 4; # take whatever is defined
1160 # when user set globals, they would interfere with our calculation, so
1161 # disable then and later re-enable them
1163 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1164 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1165 # we also need to disable any set A or P on $x (_find_round_parameters took
1166 # them already into account), since these would interfere, too
1167 delete $x->{_a}; delete $x->{_p};
1168 # need to disable $upgrade in BigInt, to avoid deep recursion
1169 local $Math::BigInt::upgrade = undef;
1171 # split the second argument into its integer and fraction part
1172 # we calculate the result then from these two parts, like in
1173 # 2 ** 2.4 == (2 ** 2) * (2 ** 0.4)
1174 my $c = $self->new($y->as_number()); # integer part
1175 my $d = $y-$c; # fractional part
1176 my $xc = $x->copy(); # a temp. copy
1178 # now calculate binary fraction from the decimal fraction on the fly
1180 # 0.654 * 2 = 1.308 > 1 => 0.1 ( 1.308 - 1 = 0.308)
1181 # 0.308 * 2 = 0.616 < 1 => 0.10
1182 # 0.616 * 2 = 1.232 > 1 => 0.101 ( 1.232 - 1 = 0.232)
1184 # The process stops when the result is exactly one, or when we have
1187 # From the binary fraction we calculate the result as follows:
1188 # we assume the fraction ends in 1, and we remove this one first.
1189 # For each digit after the dot, assume 1 eq R and 0 eq XR, where R means
1190 # take square root and X multiply with the original X.
1196 last if $d->is_one(); # == 1
1200 $x->bsqrt(); $x->bmul($xc); $d->bdec(); # 1
1204 # assume fraction ends in 1
1208 $x->bmul( $xc->bpow($c) );
1210 elsif (!$c->is_zero())
1216 # shortcut to not run trough _find_round_parameters again
1217 if (defined $params[1])
1219 $x->bround($params[1],$params[3]); # then round accordingly
1223 $x->bfround($params[2],$params[3]); # then round accordingly
1227 # clear a/p after round, since user did not request it
1228 $x->{_a} = undef; $x->{_p} = undef;
1231 $$abr = $ab; $$pbr = $pb;
1237 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1238 my ($x,$y,$a,$p,$r) = @_;
1241 # if $y == 0.5, it is sqrt($x)
1242 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1246 # Taylor: | u u^2 u^3 |
1247 # x ** y = 1 + | --- + --- + * ----- + ... |
1250 # we need to limit the accuracy to protect against overflow
1253 my @params = $x->_find_round_parameters($a,$p,$r);
1255 # no rounding at all, so must use fallback
1256 if (scalar @params == 1)
1258 # simulate old behaviour
1259 $params[1] = $self->div_scale(); # and round to it as accuracy
1260 $scale = $params[1]+4; # at least four more for proper round
1261 $params[3] = $r; # round mode by caller or undef
1262 $fallback = 1; # to clear a/p afterwards
1266 # the 4 below is empirical, and there might be cases where it is not
1268 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1271 # when user set globals, they would interfere with our calculation, so
1272 # disable then and later re-enable them
1274 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1275 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1276 # we also need to disable any set A or P on $x (_find_round_parameters took
1277 # them already into account), since these would interfere, too
1278 delete $x->{_a}; delete $x->{_p};
1279 # need to disable $upgrade in BigInt, to avoid deep recursion
1280 local $Math::BigInt::upgrade = undef;
1282 my ($limit,$v,$u,$below,$factor,$next,$over);
1284 $u = $x->copy()->blog($scale)->bmul($y);
1285 $v = $self->bone(); # 1
1286 $factor = $self->new(2); # 2
1287 $x->bone(); # first term: 1
1289 $below = $v->copy();
1292 $limit = $self->new("1E-". ($scale-1));
1296 # we calculate the next term, and add it to the last
1297 # when the next term is below our limit, it won't affect the outcome
1298 # anymore, so we stop
1299 $next = $over->copy()->bdiv($below,$scale);
1300 last if $next->bcmp($limit) <= 0;
1303 # calculate things for the next term
1304 $over *= $u; $below *= $factor; $factor->binc();
1308 # shortcut to not run trough _find_round_parameters again
1309 if (defined $params[1])
1311 $x->bround($params[1],$params[3]); # then round accordingly
1315 $x->bfround($params[2],$params[3]); # then round accordingly
1319 # clear a/p after round, since user did not request it
1320 $x->{_a} = undef; $x->{_p} = undef;
1323 $$abr = $ab; $$pbr = $pb;
1329 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1330 # compute power of two numbers, second arg is used as integer
1331 # modifies first argument
1333 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1335 return $x if $x->{sign} =~ /^[+-]inf$/;
1336 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1337 return $x->bone() if $y->is_zero();
1338 return $x if $x->is_one() || $y->is_one();
1340 return $x->_pow2($y,$a,$p,$r) if !$y->is_int(); # non-integer power
1342 my $y1 = $y->as_number(); # make bigint
1344 if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1346 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
1347 return $y1->is_odd() ? $x : $x->babs(1);
1351 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1352 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1356 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1358 $x->{_m}->bpow($y1);
1359 $x->{_e}->bmul($y1);
1360 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1362 if ($y->{sign} eq '-')
1364 # modify $x in place!
1365 my $z = $x->copy(); $x->bzero()->binc();
1366 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
1368 $x->round($a,$p,$r,$y);
1371 ###############################################################################
1372 # rounding functions
1376 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1377 # $n == 0 means round to integer
1378 # expects and returns normalized numbers!
1379 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1381 return $x if $x->modify('bfround');
1383 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1384 return $x if !defined $scale; # no-op
1386 # never round a 0, +-inf, NaN
1389 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1392 return $x if $x->{sign} !~ /^[+-]$/;
1393 # print "MBF bfround $x to scale $scale mode $mode\n";
1395 # don't round if x already has lower precision
1396 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1398 $x->{_p} = $scale; # remember round in any case
1399 $x->{_a} = undef; # and clear A
1402 # print "bfround scale $scale e $x->{_e}\n";
1403 # round right from the '.'
1404 return $x if $x->{_e} >= 0; # nothing to round
1405 $scale = -$scale; # positive for simplicity
1406 my $len = $x->{_m}->length(); # length of mantissa
1407 my $dad = -$x->{_e}; # digits after dot
1408 my $zad = 0; # zeros after dot
1409 $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
1410 #print "scale $scale dad $dad zad $zad len $len\n";
1412 # number bsstr len zad dad
1413 # 0.123 123e-3 3 0 3
1414 # 0.0123 123e-4 3 1 4
1417 # 1.2345 12345e-4 5 0 4
1419 # do not round after/right of the $dad
1420 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1422 # round to zero if rounding inside the $zad, but not for last zero like:
1423 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1424 return $x->bzero() if $scale < $zad;
1425 if ($scale == $zad) # for 0.006, scale -3 and trunc
1431 # adjust round-point to be inside mantissa
1434 $scale = $scale-$zad;
1438 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
1439 $scale = $dbd+$scale;
1442 # print "round to $x->{_m} to $scale\n";
1446 # 123 => 100 means length(123) = 3 - $scale (2) => 1
1448 my $dbt = $x->{_m}->length();
1450 my $dbd = $dbt + $x->{_e};
1451 # should be the same, so treat it as this
1452 $scale = 1 if $scale == 0;
1453 # shortcut if already integer
1454 return $x if $scale == 1 && $dbt <= $dbd;
1455 # maximum digits before dot
1460 # not enough digits before dot, so round to zero
1463 elsif ( $scale == $dbd )
1470 $scale = $dbd - $scale;
1474 # print "using $scale for $x->{_m} with '$mode'\n";
1475 # pass sign to bround for rounding modes '+inf' and '-inf'
1476 $x->{_m}->{sign} = $x->{sign};
1477 $x->{_m}->bround($scale,$mode);
1478 $x->{_m}->{sign} = '+'; # fix sign back
1484 # accuracy: preserve $N digits, and overwrite the rest with 0's
1485 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1487 die ('bround() needs positive accuracy') if ($_[0] || 0) < 0;
1489 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1490 return $x if !defined $scale; # no-op
1492 return $x if $x->modify('bround');
1494 # scale is now either $x->{_a}, $accuracy, or the user parameter
1495 # test whether $x already has lower accuracy, do nothing in this case
1496 # but do round if the accuracy is the same, since a math operation might
1497 # want to round a number with A=5 to 5 digits afterwards again
1498 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1500 # scale < 0 makes no sense
1501 # never round a +-inf, NaN
1502 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
1504 # 1: $scale == 0 => keep all digits
1505 # 2: never round a 0
1506 # 3: if we should keep more digits than the mantissa has, do nothing
1507 if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1509 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1513 # pass sign to bround for '+inf' and '-inf' rounding modes
1514 $x->{_m}->{sign} = $x->{sign};
1515 $x->{_m}->bround($scale,$mode); # round mantissa
1516 $x->{_m}->{sign} = '+'; # fix sign back
1517 # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1518 $x->{_a} = $scale; # remember rounding
1519 $x->{_p} = undef; # and clear P
1520 $x->bnorm(); # del trailing zeros gen. by bround()
1525 # return integer less or equal then $x
1526 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1528 return $x if $x->modify('bfloor');
1530 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1532 # if $x has digits after dot
1533 if ($x->{_e}->{sign} eq '-')
1535 #$x->{_m}->brsft(-$x->{_e},10);
1537 #$x-- if $x->{sign} eq '-';
1539 $x->{_e}->{sign} = '+'; # negate e
1540 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1541 $x->{_e}->bzero(); # trunc/norm
1542 $x->{_m}->binc() if $x->{sign} eq '-'; # decrement if negative
1544 $x->round($a,$p,$r);
1549 # return integer greater or equal then $x
1550 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1552 return $x if $x->modify('bceil');
1553 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1555 # if $x has digits after dot
1556 if ($x->{_e}->{sign} eq '-')
1558 #$x->{_m}->brsft(-$x->{_e},10);
1560 #$x++ if $x->{sign} eq '+';
1562 $x->{_e}->{sign} = '+'; # negate e
1563 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1564 $x->{_e}->bzero(); # trunc/norm
1565 $x->{_m}->binc() if $x->{sign} eq '+'; # decrement if negative
1567 $x->round($a,$p,$r);
1572 # shift right by $y (divide by power of 2)
1573 my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1575 return $x if $x->modify('brsft');
1576 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1578 $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1579 $x->bdiv($n ** $y,$a,$p,$r,$y);
1584 # shift right by $y (divide by power of 2)
1585 my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1587 return $x if $x->modify('brsft');
1588 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1590 $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1591 $x->bmul($n ** $y,$a,$p,$r,$y);
1594 ###############################################################################
1598 # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1603 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
1604 # or falling back to MBI::bxxx()
1605 my $name = $AUTOLOAD;
1607 $name =~ s/.*:://; # split package
1609 if (!method_alias($name))
1613 # delayed load of Carp and avoid recursion
1615 Carp::croak ("Can't call a method without name");
1617 if (!method_hand_up($name))
1619 # delayed load of Carp and avoid recursion
1621 Carp::croak ("Can't call $class\-\>$name, not a valid method");
1623 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1625 return &{"$MBI"."::$name"}(@_);
1627 my $bname = $name; $bname =~ s/^f/b/;
1628 *{$class."::$name"} = \&$bname;
1634 # return a copy of the exponent
1635 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1637 if ($x->{sign} !~ /^[+-]$/)
1639 my $s = $x->{sign}; $s =~ s/^[+-]//;
1640 return $self->new($s); # -inf, +inf => +inf
1642 return $x->{_e}->copy();
1647 # return a copy of the mantissa
1648 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1650 if ($x->{sign} !~ /^[+-]$/)
1652 my $s = $x->{sign}; $s =~ s/^[+]//;
1653 return $self->new($s); # -inf, +inf => +inf
1655 my $m = $x->{_m}->copy(); # faster than going via bstr()
1656 $m->bneg() if $x->{sign} eq '-';
1663 # return a copy of both the exponent and the mantissa
1664 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1666 if ($x->{sign} !~ /^[+-]$/)
1668 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
1669 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
1671 my $m = $x->{_m}->copy(); # faster than going via bstr()
1672 $m->bneg() if $x->{sign} eq '-';
1673 return ($m,$x->{_e}->copy());
1676 ##############################################################################
1677 # private stuff (internal use only)
1683 my $lib = ''; my @a;
1684 for ( my $i = 0; $i < $l ; $i++)
1686 # print "at $_[$i] (",$_[$i+1]||'undef',")\n";
1687 if ( $_[$i] eq ':constant' )
1689 # this rest causes overlord er load to step in
1690 # print "overload @_\n";
1691 overload::constant float => sub { $self->new(shift); };
1693 elsif ($_[$i] eq 'upgrade')
1695 # this causes upgrading
1696 $upgrade = $_[$i+1]; # or undef to disable
1699 elsif ($_[$i] eq 'downgrade')
1701 # this causes downgrading
1702 $downgrade = $_[$i+1]; # or undef to disable
1705 elsif ($_[$i] eq 'lib')
1707 $lib = $_[$i+1] || ''; # default Calc
1710 elsif ($_[$i] eq 'with')
1712 $MBI = $_[$i+1] || 'Math::BigInt'; # default Math::BigInt
1722 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
1723 my $mbilib = eval { Math::BigInt->config()->{lib} };
1724 if ((defined $mbilib) && ($MBI eq 'Math::BigInt'))
1726 # MBI already loaded
1727 $MBI->import('lib',"$lib,$mbilib", 'objectify');
1731 # MBI not loaded, or with ne "Math::BigInt"
1732 $lib .= ",$mbilib" if defined $mbilib;
1734 # my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
1735 # my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
1736 # $file = File::Spec->catfile (@parts, $file);
1740 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
1741 # used in the same script, or eval inside import().
1742 my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
1743 my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
1744 $file = File::Spec->catfile (@parts, $file);
1745 eval { require $file; $MBI->import( lib => '$lib', 'objectify' ); }
1749 my $rc = "use $MBI lib => '$lib', 'objectify';";
1753 die ("Couldn't load $MBI: $! $@") if $@;
1755 # any non :constant stuff is handled by our parent, Exporter
1756 # even if @_ is empty, to give it a chance
1757 $self->SUPER::import(@a); # for subclasses
1758 $self->export_to_level(1,$self,@a); # need this, too
1763 # adjust m and e so that m is smallest possible
1764 # round number according to accuracy and precision settings
1765 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1767 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
1769 # if (!$x->{_m}->is_odd())
1771 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
1774 $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros);
1776 # for something like 0Ey, set y to 1, and -0 => +0
1777 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
1779 # this is to prevent automatically rounding when MBI's globals are set
1780 $x->{_m}->{_f} = MB_NEVER_ROUND;
1781 $x->{_e}->{_f} = MB_NEVER_ROUND;
1782 # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
1783 $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
1784 $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
1785 $x; # MBI bnorm is no-op, so dont call it
1788 ##############################################################################
1789 # internal calculation routines
1793 # return copy as a bigint representation of this BigFloat number
1794 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1796 my $z = $x->{_m}->copy();
1797 if ($x->{_e}->{sign} eq '-') # < 0
1799 $x->{_e}->{sign} = '+'; # flip
1800 $z->brsft($x->{_e},10);
1801 $x->{_e}->{sign} = '-'; # flip back
1803 elsif (!$x->{_e}->is_zero()) # > 0
1805 $z->blsft($x->{_e},10);
1807 $z->{sign} = $x->{sign};
1814 my $class = ref($x) || $x;
1815 $x = $class->new(shift) unless ref($x);
1817 return 1 if $x->{_m}->is_zero();
1818 my $len = $x->{_m}->length();
1819 $len += $x->{_e} if $x->{_e}->sign() eq '+';
1822 my $t = $MBI->bzero();
1823 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1834 Math::BigFloat - Arbitrary size floating point math package
1841 $x = Math::BigFloat->new($str); # defaults to 0
1842 $nan = Math::BigFloat->bnan(); # create a NotANumber
1843 $zero = Math::BigFloat->bzero(); # create a +0
1844 $inf = Math::BigFloat->binf(); # create a +inf
1845 $inf = Math::BigFloat->binf('-'); # create a -inf
1846 $one = Math::BigFloat->bone(); # create a +1
1847 $one = Math::BigFloat->bone('-'); # create a -1
1850 $x->is_zero(); # true if arg is +0
1851 $x->is_nan(); # true if arg is NaN
1852 $x->is_one(); # true if arg is +1
1853 $x->is_one('-'); # true if arg is -1
1854 $x->is_odd(); # true if odd, false for even
1855 $x->is_even(); # true if even, false for odd
1856 $x->is_positive(); # true if >= 0
1857 $x->is_negative(); # true if < 0
1858 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
1860 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
1861 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
1862 $x->sign(); # return the sign, either +,- or NaN
1863 $x->digit($n); # return the nth digit, counting from right
1864 $x->digit(-$n); # return the nth digit, counting from left
1866 # The following all modify their first argument:
1869 $x->bzero(); # set $i to 0
1870 $x->bnan(); # set $i to NaN
1871 $x->bone(); # set $x to +1
1872 $x->bone('-'); # set $x to -1
1873 $x->binf(); # set $x to inf
1874 $x->binf('-'); # set $x to -inf
1876 $x->bneg(); # negation
1877 $x->babs(); # absolute value
1878 $x->bnorm(); # normalize (no-op)
1879 $x->bnot(); # two's complement (bit wise not)
1880 $x->binc(); # increment x by 1
1881 $x->bdec(); # decrement x by 1
1883 $x->badd($y); # addition (add $y to $x)
1884 $x->bsub($y); # subtraction (subtract $y from $x)
1885 $x->bmul($y); # multiplication (multiply $x by $y)
1886 $x->bdiv($y); # divide, set $i to quotient
1887 # return (quo,rem) or quo if scalar
1889 $x->bmod($y); # modulus
1890 $x->bpow($y); # power of arguments (a**b)
1891 $x->blsft($y); # left shift
1892 $x->brsft($y); # right shift
1893 # return (quo,rem) or quo if scalar
1895 $x->blog($base); # logarithm of $x, base defaults to e
1896 # (other bases than e not supported yet)
1898 $x->band($y); # bit-wise and
1899 $x->bior($y); # bit-wise inclusive or
1900 $x->bxor($y); # bit-wise exclusive or
1901 $x->bnot(); # bit-wise not (two's complement)
1903 $x->bsqrt(); # calculate square-root
1904 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
1906 $x->bround($N); # accuracy: preserver $N digits
1907 $x->bfround($N); # precision: round to the $Nth digit
1909 # The following do not modify their arguments:
1910 bgcd(@values); # greatest common divisor
1911 blcm(@values); # lowest common multiplicator
1913 $x->bstr(); # return string
1914 $x->bsstr(); # return string in scientific notation
1916 $x->bfloor(); # return integer less or equal than $x
1917 $x->bceil(); # return integer greater or equal than $x
1919 $x->exponent(); # return exponent as BigInt
1920 $x->mantissa(); # return mantissa as BigInt
1921 $x->parts(); # return (mantissa,exponent) as BigInt
1923 $x->length(); # number of digits (w/o sign and '.')
1924 ($l,$f) = $x->length(); # number of digits, and length of fraction
1928 All operators (inlcuding basic math operations) are overloaded if you
1929 declare your big floating point numbers as
1931 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1933 Operations with overloaded operators preserve the arguments, which is
1934 exactly what you expect.
1936 =head2 Canonical notation
1938 Input to these routines are either BigFloat objects, or strings of the
1939 following four forms:
1953 C</^[+-]\d+E[+-]?\d+$/>
1957 C</^[+-]\d*\.\d+E[+-]?\d+$/>
1961 all with optional leading and trailing zeros and/or spaces. Additonally,
1962 numbers are allowed to have an underscore between any two digits.
1964 Empty strings as well as other illegal numbers results in 'NaN'.
1966 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1967 are always stored in normalized form. On a string, it creates a BigFloat
1972 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1974 The string output will always have leading and trailing zeros stripped and drop
1975 a plus sign. C<bstr()> will give you always the form with a decimal point,
1976 while C<bsstr()> (for scientific) gives you the scientific notation.
1978 Input bstr() bsstr()
1980 ' -123 123 123' '-123123123' '-123123123E0'
1981 '00.0123' '0.0123' '123E-4'
1982 '123.45E-2' '1.2345' '12345E-4'
1983 '10E+3' '10000' '1E4'
1985 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1986 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1987 return either undef, <0, 0 or >0 and are suited for sort.
1989 Actual math is done by using BigInts to represent the mantissa and exponent.
1990 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1991 represent the result when input arguments are not numbers, as well as
1992 the result of dividing by zero.
1994 =head2 C<mantissa()>, C<exponent()> and C<parts()>
1996 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1997 as BigInts such that:
1999 $m = $x->mantissa();
2000 $e = $x->exponent();
2001 $y = $m * ( 10 ** $e );
2002 print "ok\n" if $x == $y;
2004 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2006 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2008 Currently the mantissa is reduced as much as possible, favouring higher
2009 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2010 This might change in the future, so do not depend on it.
2012 =head2 Accuracy vs. Precision
2014 See also: L<Rounding|Rounding>.
2016 Math::BigFloat supports both precision and accuracy. For a full documentation,
2017 examples and tips on these topics please see the large section in
2020 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2021 a operation consumes all resources, each operation produces no more than
2022 C<Math::BigFloat::precision()> digits.
2024 In case the result of one operation has more precision than specified,
2025 it is rounded. The rounding mode taken is either the default mode, or the one
2026 supplied to the operation after the I<scale>:
2028 $x = Math::BigFloat->new(2);
2029 Math::BigFloat::precision(5); # 5 digits max
2030 $y = $x->copy()->bdiv(3); # will give 0.66666
2031 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2032 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
2033 Math::BigFloat::round_mode('zero');
2034 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2040 =item ffround ( +$scale )
2042 Rounds to the $scale'th place left from the '.', counting from the dot.
2043 The first digit is numbered 1.
2045 =item ffround ( -$scale )
2047 Rounds to the $scale'th place right from the '.', counting from the dot.
2051 Rounds to an integer.
2053 =item fround ( +$scale )
2055 Preserves accuracy to $scale digits from the left (aka significant digits)
2056 and pads the rest with zeros. If the number is between 1 and -1, the
2057 significant digits count from the first non-zero after the '.'
2059 =item fround ( -$scale ) and fround ( 0 )
2061 These are effetively no-ops.
2065 All rounding functions take as a second parameter a rounding mode from one of
2066 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2068 The default rounding mode is 'even'. By using
2069 C<< Math::BigFloat::round_mode($round_mode); >> you can get and set the default
2070 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2071 no longer supported.
2072 The second parameter to the round functions then overrides the default
2075 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
2076 'trunc' as rounding mode to make it equivalent to:
2081 You can override this by passing the desired rounding mode as parameter to
2084 $x = Math::BigFloat->new(2.5);
2085 $y = $x->as_number('odd'); # $y = 3
2091 =head1 Autocreating constants
2093 After C<use Math::BigFloat ':constant'> all the floating point constants
2094 in the given scope are converted to C<Math::BigFloat>. This conversion
2095 happens at compile time.
2099 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2101 prints the value of C<2E-100>. Note that without conversion of
2102 constants the expression 2E-100 will be calculated as normal floating point
2105 Please note that ':constant' does not affect integer constants, nor binary
2106 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2111 Math with the numbers is done (by default) by a module called
2112 Math::BigInt::Calc. This is equivalent to saying:
2114 use Math::BigFloat lib => 'Calc';
2116 You can change this by using:
2118 use Math::BigFloat lib => 'BitVect';
2120 The following would first try to find Math::BigInt::Foo, then
2121 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2123 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2125 Calc.pm uses as internal format an array of elements of some decimal base
2126 (usually 1e7, but this might be differen for some systems) with the least
2127 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2128 significant bit first. Other modules might use even different means of
2129 representing the numbers. See the respective module documentation for further
2132 Please note that Math::BigFloat does B<not> use the denoted library itself,
2133 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2136 use Math::BigInt lib => 'GMP';
2139 you can roll it all into one line:
2141 use Math::BigFloat lib => 'GMP';
2143 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details.
2145 =head2 Using Math::BigInt::Lite
2147 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2150 use Math::BigFloat with => 'Math::BigInt::Lite';
2152 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2153 can combine these if you want. For instance, you may want to use
2154 Math::BigInt objects in your main script, too.
2158 use Math::BigFloat with => 'Math::BigInt::Lite';
2160 Of course, you can combine this with the C<lib> parameter.
2163 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2165 If you want to use Math::BigInt's, too, simple add a Math::BigInt B<before>:
2169 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2171 Notice that the module with the last C<lib> will "win" and thus
2172 it's lib will be used if the lib is available:
2175 use Math::BigInt lib => 'Bar,Baz';
2176 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2178 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2179 words, Math::BigFloat will try to retain previously loaded libs when you
2180 don't specify it one.
2182 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2183 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2184 same as trying the latter load alone, except for the fact that Bar or Baz
2185 might be loaded needlessly in an intermidiate step
2187 The old way still works though:
2190 use Math::BigInt lib => 'Bar,Baz';
2193 But B<examples #3 and #4 are recommended> for usage.
2201 The following does not work yet:
2203 $m = $x->mantissa();
2204 $e = $x->exponent();
2205 $y = $m * ( 10 ** $e );
2206 print "ok\n" if $x == $y;
2210 There is no fmod() function yet.
2218 =item stringify, bstr()
2220 Both stringify and bstr() now drop the leading '+'. The old code would return
2221 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2222 reasoning and details.
2226 The following will probably not do what you expect:
2228 print $c->bdiv(123.456),"\n";
2230 It prints both quotient and reminder since print works in list context. Also,
2231 bdiv() will modify $c, so be carefull. You probably want to use
2233 print $c / 123.456,"\n";
2234 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2238 =item Modifying and =
2242 $x = Math::BigFloat->new(5);
2245 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2246 a second reference to the B<same> object and stores it in $y. Thus anything
2247 that modifies $x will modify $y, and vice versa.
2250 print "$x, $y\n"; # prints '10, 10'
2252 If you want a true copy of $x, use:
2256 See also the documentation in L<overload> regarding C<=>.
2260 C<bpow()> now modifies the first argument, unlike the old code which left
2261 it alone and only returned the result. This is to be consistent with
2262 C<badd()> etc. The first will modify $x, the second one won't:
2264 print bpow($x,$i),"\n"; # modify $x
2265 print $x->bpow($i),"\n"; # ditto
2266 print $x ** $i,"\n"; # leave $x alone
2272 This program is free software; you may redistribute it and/or modify it under
2273 the same terms as Perl itself.
2277 Mark Biggar, overloaded interface by Ilya Zakharevich.
2278 Completely rewritten by Tels http://bloodgate.com in 2001.