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 $MBI/;
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 $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);
207 ##############################################################################
208 # string conversation
212 # (ref to BFLOAT or num_str ) return num_str
213 # Convert number from internal format to (non-scientific) string format.
214 # internal format is always normalized (no leading zeros, "-0" => "+0")
215 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
216 #my $x = shift; my $class = ref($x) || $x;
217 #$x = $class->new(shift) unless ref($x);
219 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
220 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
221 if ($x->{sign} !~ /^[+-]$/)
223 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
227 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
229 my $not_zero = ! $x->is_zero();
232 $es = $x->{_m}->bstr();
233 $len = CORE::length($es);
234 if (!$x->{_e}->is_zero())
236 if ($x->{_e}->sign() eq '-')
239 if ($x->{_e} <= -$len)
241 # print "style: 0.xxxx\n";
242 my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
243 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
247 # print "insert '.' at $x->{_e} in '$es'\n";
248 substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
254 $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
258 $es = $x->{sign}.$es if $x->{sign} eq '-';
259 # if set accuracy or precision, pad with zeros
260 if ((defined $x->{_a}) && ($not_zero))
262 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
263 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
264 $zeros = $x->{_a} - $len if $cad != $len;
265 $es .= $dot.'0' x $zeros if $zeros > 0;
267 elsif ($x->{_p} || 0 < 0)
269 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
270 my $zeros = -$x->{_p} + $cad;
271 $es .= $dot.'0' x $zeros if $zeros > 0;
278 # (ref to BFLOAT or num_str ) return num_str
279 # Convert number from internal format to scientific string format.
280 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
281 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
282 #my $x = shift; my $class = ref($x) || $x;
283 #$x = $class->new(shift) unless ref($x);
285 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
286 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
287 if ($x->{sign} !~ /^[+-]$/)
289 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
292 my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
294 $x->{_m}->bstr().$sep.$x->{_e}->bstr();
299 # Make a number from a BigFloat object
300 # simple return string and let Perl's atoi()/atof() handle the rest
301 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
305 ##############################################################################
306 # public stuff (usually prefixed with "b")
309 # todo: this must be overwritten and return NaN for non-integer values
310 # band(), bior(), bxor(), too
313 # $class->SUPER::bnot($class,@_);
318 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
319 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
320 my ($self,$x,$y) = objectify(2,@_);
322 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
324 # handle +-inf and NaN
325 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
326 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
327 return +1 if $x->{sign} eq '+inf';
328 return -1 if $x->{sign} eq '-inf';
329 return -1 if $y->{sign} eq '+inf';
333 # check sign for speed first
334 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
335 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
338 my $xz = $x->is_zero();
339 my $yz = $y->is_zero();
340 return 0 if $xz && $yz; # 0 <=> 0
341 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
342 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
344 # adjust so that exponents are equal
345 my $lxm = $x->{_m}->length();
346 my $lym = $y->{_m}->length();
347 # the numify somewhat limits our length, but makes it much faster
348 my $lx = $lxm + $x->{_e}->numify();
349 my $ly = $lym + $y->{_e}->numify();
350 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
351 return $l <=> 0 if $l != 0;
353 # lengths (corrected by exponent) are equal
354 # so make mantissa equal length by padding with zero (shift left)
355 my $diff = $lxm - $lym;
356 my $xm = $x->{_m}; # not yet copy it
360 $ym = $y->{_m}->copy()->blsft($diff,10);
364 $xm = $x->{_m}->copy()->blsft(-$diff,10);
366 my $rc = $xm->bacmp($ym);
367 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
373 # Compares 2 values, ignoring their signs.
374 # Returns one of undef, <0, =0, >0. (suitable for sort)
375 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
376 my ($self,$x,$y) = objectify(2,@_);
378 # handle +-inf and NaN's
379 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
381 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
382 return 0 if ($x->is_inf() && $y->is_inf());
383 return 1 if ($x->is_inf() && !$y->is_inf());
388 my $xz = $x->is_zero();
389 my $yz = $y->is_zero();
390 return 0 if $xz && $yz; # 0 <=> 0
391 return -1 if $xz && !$yz; # 0 <=> +y
392 return 1 if $yz && !$xz; # +x <=> 0
394 # adjust so that exponents are equal
395 my $lxm = $x->{_m}->length();
396 my $lym = $y->{_m}->length();
397 # the numify somewhat limits our length, but makes it much faster
398 my $lx = $lxm + $x->{_e}->numify();
399 my $ly = $lym + $y->{_e}->numify();
401 return $l <=> 0 if $l != 0;
403 # lengths (corrected by exponent) are equal
404 # so make mantissa equal-length by padding with zero (shift left)
405 my $diff = $lxm - $lym;
406 my $xm = $x->{_m}; # not yet copy it
410 $ym = $y->{_m}->copy()->blsft($diff,10);
414 $xm = $x->{_m}->copy()->blsft(-$diff,10);
416 $xm->bacmp($ym) <=> 0;
421 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
422 # return result as BFLOAT
423 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
425 #print "mbf badd $x $y\n";
426 # inf and NaN handling
427 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
430 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
432 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
434 # +inf++inf or -inf+-inf => same, rest is NaN
435 return $x if $x->{sign} eq $y->{sign};
438 # +-inf + something => +inf; something +-inf => +-inf
439 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
443 # speed: no add for 0+y or x+0
444 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
445 if ($x->is_zero()) # 0+y
447 # make copy, clobbering up x (modify in place!)
448 $x->{_e} = $y->{_e}->copy();
449 $x->{_m} = $y->{_m}->copy();
450 $x->{sign} = $y->{sign} || $nan;
451 return $x->round($a,$p,$r,$y);
454 # take lower of the two e's and adapt m1 to it to match m2
456 $e = $MBI->bzero() if !defined $e; # if no BFLOAT ?
457 $e = $e->copy(); # make copy (didn't do it yet)
459 my $add = $y->{_m}->copy();
460 if ($e->{sign} eq '-') # < 0
462 my $e1 = $e->copy()->babs();
463 #$x->{_m} *= (10 ** $e1);
464 $x->{_m}->blsft($e1,10);
465 $x->{_e} += $e; # need the sign of e
467 elsif (!$e->is_zero()) # > 0
472 # else: both e are the same, so just leave them
473 $x->{_m}->{sign} = $x->{sign}; # fiddle with signs
474 $add->{sign} = $y->{sign};
475 $x->{_m} += $add; # finally do add/sub
476 $x->{sign} = $x->{_m}->{sign}; # re-adjust signs
477 $x->{_m}->{sign} = '+'; # mantissa always positiv
478 # delete trailing zeros, then round
479 return $x->bnorm()->round($a,$p,$r,$y);
484 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
485 # subtract second arg from first, modify first
486 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
488 if ($y->is_zero()) # still round for not adding zero
490 return $x->round($a,$p,$r);
493 $y->{sign} =~ tr/+\-/-+/; # does nothing for NaN
494 $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
495 $y->{sign} =~ tr/+\-/-+/; # refix $y (does nothing for NaN)
496 $x; # already rounded by badd()
501 # increment arg by one
502 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
504 if ($x->{_e}->sign() eq '-')
506 return $x->badd($self->bone(),$a,$p,$r); # digits after dot
509 if (!$x->{_e}->is_zero())
511 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
515 if ($x->{sign} eq '+')
518 return $x->bnorm()->bround($a,$p,$r);
520 elsif ($x->{sign} eq '-')
523 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
524 return $x->bnorm()->bround($a,$p,$r);
526 # inf, nan handling etc
527 $x->badd($self->__one(),$a,$p,$r); # does round
532 # decrement arg by one
533 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
535 if ($x->{_e}->sign() eq '-')
537 return $x->badd($self->bone('-'),$a,$p,$r); # digits after dot
540 if (!$x->{_e}->is_zero())
542 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
546 my $zero = $x->is_zero();
548 if (($x->{sign} eq '-') || $zero)
551 $x->{sign} = '-' if $zero; # 0 => 1 => -1
552 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
553 return $x->bnorm()->round($a,$p,$r);
556 elsif ($x->{sign} eq '+')
559 return $x->bnorm()->round($a,$p,$r);
561 # inf, nan handling etc
562 $x->badd($self->bone('-'),$a,$p,$r); # does round
567 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(2,@_);
569 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
573 # Taylor: | u 1 u^3 1 u^5 |
574 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
575 # |_ v 3 v^3 5 v^5 _|
577 # This takes much more steps to calculate the result:
580 # Taylor: | u 1 u^2 1 u^3 |
581 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
582 # |_ x 2 x^2 3 x^3 _|
584 # we need to limit the accuracy to protect against overflow
587 my @params = $x->_find_round_parameters($a,$p,$r);
589 # no rounding at all, so must use fallback
590 if (scalar @params == 1)
592 # simulate old behaviour
593 $params[1] = $self->div_scale(); # and round to it as accuracy
594 $scale = $params[1]+4; # at least four more for proper round
595 $params[3] = $r; # round mode by caller or undef
596 $fallback = 1; # to clear a/p afterwards
600 # the 4 below is empirical, and there might be cases where it is not
602 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
605 return $x->bzero(@params) if $x->is_one();
606 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
607 #return $x->bone('+',@params) if $x->bcmp($base) == 0;
609 # when user set globals, they would interfere with our calculation, so
610 # disable then and later re-enable them
612 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
613 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
614 # we also need to disable any set A or P on $x (_find_round_parameters took
615 # them already into account), since these would interfere, too
616 delete $x->{_a}; delete $x->{_p};
617 # need to disable $upgrade in BigInt, to avoid deep recursion
618 local $Math::BigInt::upgrade = undef;
620 my ($case,$limit,$v,$u,$below,$factor,$two,$next,$over,$f);
623 #if ($x <= Math::BigFloat->new("0.5"))
626 # print "case $case $x < 0.5\n";
627 $v = $x->copy(); $v->binc(); # v = x+1
628 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
629 $x->bdiv($v,$scale); # first term: u/v
632 $u *= $u; $v *= $v; # u^2, v^2
633 $below->bmul($v); # u^3, v^3
635 $factor = $self->new(3); $f = $self->new(2);
640 # print "case 1 $x > 0.5\n";
641 # $v = $x->copy(); # v = x
642 # $u = $x->copy(); $u->bdec(); # u = x-1;
643 # $x->bdec(); $x->bdiv($v,$scale); # first term: x-1/x
644 # $below = $v->copy();
645 # $over = $u->copy();
646 # $below->bmul($v); # u^2, v^2
648 # $factor = $self->new(2); $f = $self->bone();
650 $limit = $self->new("1E-". ($scale-1));
654 # we calculate the next term, and add it to the last
655 # when the next term is below our limit, it won't affect the outcome
656 # anymore, so we stop
657 $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
658 last if $next->bcmp($limit) <= 0;
660 # print "step $steps $x\n";
661 # calculate things for the next term
662 $over *= $u; $below *= $v; $factor->badd($f);
665 $x->bmul(2) if $case == 0;
666 #print "took $steps steps\n";
668 # shortcut to not run trough _find_round_parameters again
669 if (defined $params[1])
671 $x->bround($params[1],$params[3]); # then round accordingly
675 $x->bfround($params[2],$params[3]); # then round accordingly
679 # clear a/p after round, since user did not request it
680 $x->{_a} = undef; $x->{_p} = undef;
683 $$abr = $ab; $$pbr = $pb;
690 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
691 # does not modify arguments, but returns new object
692 # Lowest Common Multiplicator
694 my ($self,@arg) = objectify(0,@_);
695 my $x = $self->new(shift @arg);
696 while (@arg) { $x = _lcm($x,shift @arg); }
702 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
703 # does not modify arguments, but returns new object
704 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
706 my ($self,@arg) = objectify(0,@_);
707 my $x = $self->new(shift @arg);
708 while (@arg) { $x = _gcd($x,shift @arg); }
712 ###############################################################################
713 # is_foo methods (is_negative, is_positive are inherited from BigInt)
717 # return true if arg (BFLOAT or num_str) is an integer
718 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
720 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
721 $x->{_e}->{sign} eq '+'; # 1e-1 => no integer
727 # return true if arg (BFLOAT or num_str) is zero
728 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
730 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
736 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
737 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
739 my $sign = shift || ''; $sign = '+' if $sign ne '-';
741 if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
747 # return true if arg (BFLOAT or num_str) is odd or false if even
748 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
750 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
751 ($x->{_e}->is_zero() && $x->{_m}->is_odd());
757 # return true if arg (BINT or num_str) is even or false if odd
758 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
760 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
761 return 1 if ($x->{_e}->{sign} eq '+' # 123.45 is never
762 && $x->{_m}->is_even()); # but 1200 is
768 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
769 # (BINT or num_str, BINT or num_str) return BINT
770 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
772 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
775 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
777 return $x->bnan() if $x->is_zero() || $y->is_zero();
778 # result will always be +-inf:
779 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
780 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
781 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
782 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
783 return $x->binf('-');
786 return $x->bzero() if $x->is_zero() || $y->is_zero();
788 # aEb * cEd = (a*c)E(b+d)
789 $x->{_m}->bmul($y->{_m});
790 $x->{_e}->badd($y->{_e});
792 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
793 return $x->bnorm()->round($a,$p,$r,$y);
798 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
799 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
800 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
802 return $self->_div_inf($x,$y)
803 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
805 # x== 0 # also: or y == 1 or y == -1
806 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
809 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
811 # we need to limit the accuracy to protect against overflow
814 my @params = $x->_find_round_parameters($a,$p,$r,$y);
816 # no rounding at all, so must use fallback
817 if (scalar @params == 1)
819 # simulate old behaviour
820 $params[1] = $self->div_scale(); # and round to it as accuracy
821 $scale = $params[1]+4; # at least four more for proper round
822 $params[3] = $r; # round mode by caller or undef
823 $fallback = 1; # to clear a/p afterwards
827 # the 4 below is empirical, and there might be cases where it is not
829 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
831 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
832 $scale = $lx if $lx > $scale;
833 $scale = $ly if $ly > $scale;
834 my $diff = $ly - $lx;
835 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
837 # make copy of $x in case of list context for later reminder calculation
839 if (wantarray && !$y->is_one())
844 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
846 # check for / +-1 ( +/- 1E0)
849 # promote BigInts and it's subclasses (except when already a BigFloat)
850 $y = $self->new($y) unless $y->isa('Math::BigFloat');
852 #print "bdiv $y ",ref($y),"\n";
853 # need to disable $upgrade in BigInt, to avoid deep recursion
854 local $Math::BigInt::upgrade = undef; # should be parent class vs MBI
856 # calculate the result to $scale digits and then round it
857 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
858 $x->{_m}->blsft($scale,10);
859 $x->{_m}->bdiv( $y->{_m} ); # a/c
860 $x->{_e}->bsub( $y->{_e} ); # b-d
861 $x->{_e}->bsub($scale); # correct for 10**scale
862 $x->bnorm(); # remove trailing 0's
865 # shortcut to not run trough _find_round_parameters again
866 if (defined $params[1])
868 $x->bround($params[1],$params[3]); # then round accordingly
872 $x->bfround($params[2],$params[3]); # then round accordingly
876 # clear a/p after round, since user did not request it
877 $x->{_a} = undef; $x->{_p} = undef;
884 $rem->bmod($y,$params[1],$params[2],$params[3]); # copy already done
888 $rem = $self->bzero();
892 # clear a/p after round, since user did not request it
893 $rem->{_a} = undef; $rem->{_p} = undef;
902 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
903 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
905 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
907 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
908 return $re->round($a,$p,$r,$y);
910 return $x->bnan() if $x->is_zero() && $y->is_zero();
911 return $x if $y->is_zero();
912 return $x->bnan() if $x->is_nan() || $y->is_nan();
913 return $x->bzero() if $y->is_one() || $x->is_zero();
915 # inf handling is missing here
917 my $cmp = $x->bacmp($y); # equal or $x < $y?
918 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
920 # only $y of the operands negative?
921 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
923 $x->{sign} = $y->{sign}; # calc sign first
924 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
926 my $ym = $y->{_m}->copy();
929 $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
931 # if $y has digits after dot
932 my $shifty = 0; # correct _e of $x by this
933 if ($y->{_e}->{sign} eq '-') # has digits after dot
935 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
936 $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
937 $x->blsft($shifty,10); # 123 => 1230, $y->{_m} is already 25
939 # $ym is now mantissa of $y based on exponent 0
941 my $shiftx = 0; # correct _e of $x by this
942 if ($x->{_e}->{sign} eq '-') # has digits after dot
944 # 123.4 % 20 => 1234 % 200
945 $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
946 $ym->blsft($shiftx,10);
948 # 123e1 % 20 => 1230 % 20
949 if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
951 $x->{_m}->blsft($x->{_e},10);
953 $x->{_e} = $MBI->bzero() unless $x->{_e}->is_zero();
955 $x->{_e}->bsub($shiftx) if $shiftx != 0;
956 $x->{_e}->bsub($shifty) if $shifty != 0;
958 # now mantissas are equalized, exponent of $x is adjusted, so calc result
962 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
965 if ($neg != 0) # one of them negative => correct in place
970 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
974 $x->round($a,$p,$r,$y); # round and return
979 # calculate square root; this should probably
980 # use a different test to see whether the accuracy we want is...
981 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
983 return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
984 return $x if $x->{sign} eq '+inf'; # +inf
985 return $x if $x->is_zero() || $x->is_one();
987 # we need to limit the accuracy to protect against overflow
990 my @params = $x->_find_round_parameters($a,$p,$r);
992 # no rounding at all, so must use fallback
993 if (scalar @params == 1)
995 # simulate old behaviour
996 $params[1] = $self->div_scale(); # and round to it as accuracy
997 $scale = $params[1]+4; # at least four more for proper round
998 $params[3] = $r; # round mode by caller or undef
999 $fallback = 1; # to clear a/p afterwards
1003 # the 4 below is empirical, and there might be cases where it is not
1005 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1008 # when user set globals, they would interfere with our calculation, so
1009 # disable them and later re-enable them
1011 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1012 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1013 # we also need to disable any set A or P on $x (_find_round_parameters took
1014 # them already into account), since these would interfere, too
1015 delete $x->{_a}; delete $x->{_p};
1016 # need to disable $upgrade in BigInt, to avoid deep recursion
1017 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1019 my $xas = $x->as_number();
1020 my $gs = $xas->copy()->bsqrt(); # some guess
1022 # print "guess $gs\n";
1023 if (($x->{_e}->{sign} ne '-') # guess can't be accurate if there are
1024 # digits after the dot
1025 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
1028 $x->{_m} = $gs; $x->{_e} = $MBI->bzero(); $x->bnorm();
1029 # shortcut to not run trough _find_round_parameters again
1030 if (defined $params[1])
1032 $x->bround($params[1],$params[3]); # then round accordingly
1036 $x->bfround($params[2],$params[3]); # then round accordingly
1040 # clear a/p after round, since user did not request it
1041 $x->{_a} = undef; $x->{_p} = undef;
1043 # re-enable A and P, upgrade is taken care of by "local"
1044 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1047 $gs = $self->new( $gs ); # BigInt to BigFloat
1049 my $lx = $x->{_m}->length();
1050 $scale = $lx if $scale < $lx;
1051 my $e = $self->new("1E-$scale"); # make test variable
1054 my $two = $self->new(2);
1056 # promote BigInts and it's subclasses (except when already a BigFloat)
1057 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1060 while ($diff->bacmp($e) >= 0)
1062 $rem = $y->copy()->bdiv($gs,$scale);
1063 $rem = $y->copy()->bdiv($gs,$scale)->badd($gs)->bdiv($two,$scale);
1064 $diff = $rem->copy()->bsub($gs);
1067 # copy over to modify $x
1068 $x->{_m} = $rem->{_m}; $x->{_e} = $rem->{_e};
1070 # shortcut to not run trough _find_round_parameters again
1071 if (defined $params[1])
1073 $x->bround($params[1],$params[3]); # then round accordingly
1077 $x->bfround($params[2],$params[3]); # then round accordingly
1081 # clear a/p after round, since user did not request it
1082 $x->{_a} = undef; $x->{_p} = undef;
1085 $$abr = $ab; $$pbr = $pb;
1091 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1092 # compute factorial numbers
1093 # modifies first argument
1094 my ($self,$x,@r) = objectify(1,@_);
1097 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1098 ($x->{_e}->{sign} ne '+')); # digits after dot?
1100 return $x->bone(@r) if $x->is_zero() || $x->is_one(); # 0 or 1 => 1
1102 # use BigInt's bfac() for faster calc
1103 $x->{_m}->blsft($x->{_e},10); # un-norm m
1104 $x->{_e}->bzero(); # norm $x again
1105 $x->{_m}->bfac(); # factorial
1106 $x->bnorm()->round(@r);
1111 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1112 my ($x,$y,$a,$p,$r) = @_;
1115 # we need to limit the accuracy to protect against overflow
1118 my @params = $x->_find_round_parameters($a,$p,$r);
1120 # no rounding at all, so must use fallback
1121 if (scalar @params == 1)
1123 # simulate old behaviour
1124 $params[1] = $self->div_scale(); # and round to it as accuracy
1125 $scale = $params[1]+4; # at least four more for proper round
1126 $params[3] = $r; # round mode by caller or undef
1127 $fallback = 1; # to clear a/p afterwards
1131 # the 4 below is empirical, and there might be cases where it is not
1133 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1136 # when user set globals, they would interfere with our calculation, so
1137 # disable then and later re-enable them
1139 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1140 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1141 # we also need to disable any set A or P on $x (_find_round_parameters took
1142 # them already into account), since these would interfere, too
1143 delete $x->{_a}; delete $x->{_p};
1144 # need to disable $upgrade in BigInt, to avoid deep recursion
1145 local $Math::BigInt::upgrade = undef;
1147 # split the second argument into its integer and fraction part
1148 # we calculate the result then from these two parts, like in
1149 # 2 ** 2.4 == (2 ** 2) * (2 ** 0.4)
1150 my $c = $self->new($y->as_number()); # integer part
1151 my $d = $y-$c; # fractional part
1152 my $xc = $x->copy(); # a temp. copy
1154 # now calculate binary fraction from the decimal fraction on the fly
1156 # 0.654 * 2 = 1.308 > 1 => 0.1 ( 1.308 - 1 = 0.308)
1157 # 0.308 * 2 = 0.616 < 1 => 0.10
1158 # 0.616 * 2 = 1.232 > 1 => 0.101 ( 1.232 - 1 = 0.232)
1160 # The process stops when the result is exactly one, or when we have
1163 # From the binary fraction we calculate the result as follows:
1164 # we assume the fraction ends in 1, and we remove this one first.
1165 # For each digit after the dot, assume 1 eq R and 0 eq XR, where R means
1166 # take square root and X multiply with the original X.
1172 last if $d->is_one(); # == 1
1176 $x->bsqrt(); $x->bmul($xc); $d->bdec(); # 1
1180 # assume fraction ends in 1
1184 $x->bmul( $xc->bpow($c) );
1186 elsif (!$c->is_zero())
1192 # shortcut to not run trough _find_round_parameters again
1193 if (defined $params[1])
1195 $x->bround($params[1],$params[3]); # then round accordingly
1199 $x->bfround($params[2],$params[3]); # then round accordingly
1203 # clear a/p after round, since user did not request it
1204 $x->{_a} = undef; $x->{_p} = undef;
1207 $$abr = $ab; $$pbr = $pb;
1213 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1214 my ($x,$y,$a,$p,$r) = @_;
1217 # if $y == 0.5, it is sqrt($x)
1218 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1222 # Taylor: | u u^2 u^3 |
1223 # x ** y = 1 + | --- + --- + * ----- + ... |
1226 # we need to limit the accuracy to protect against overflow
1229 my @params = $x->_find_round_parameters($a,$p,$r);
1231 # no rounding at all, so must use fallback
1232 if (scalar @params == 1)
1234 # simulate old behaviour
1235 $params[1] = $self->div_scale(); # and round to it as accuracy
1236 $scale = $params[1]+4; # at least four more for proper round
1237 $params[3] = $r; # round mode by caller or undef
1238 $fallback = 1; # to clear a/p afterwards
1242 # the 4 below is empirical, and there might be cases where it is not
1244 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1247 # when user set globals, they would interfere with our calculation, so
1248 # disable then and later re-enable them
1250 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1251 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1252 # we also need to disable any set A or P on $x (_find_round_parameters took
1253 # them already into account), since these would interfere, too
1254 delete $x->{_a}; delete $x->{_p};
1255 # need to disable $upgrade in BigInt, to avoid deep recursion
1256 local $Math::BigInt::upgrade = undef;
1258 my ($limit,$v,$u,$below,$factor,$next,$over);
1260 $u = $x->copy()->blog($scale)->bmul($y);
1261 $v = $self->bone(); # 1
1262 $factor = $self->new(2); # 2
1263 $x->bone(); # first term: 1
1265 $below = $v->copy();
1268 $limit = $self->new("1E-". ($scale-1));
1272 # we calculate the next term, and add it to the last
1273 # when the next term is below our limit, it won't affect the outcome
1274 # anymore, so we stop
1275 $next = $over->copy()->bdiv($below,$scale);
1276 last if $next->bcmp($limit) <= 0;
1279 # calculate things for the next term
1280 $over *= $u; $below *= $factor; $factor->binc();
1284 # shortcut to not run trough _find_round_parameters again
1285 if (defined $params[1])
1287 $x->bround($params[1],$params[3]); # then round accordingly
1291 $x->bfround($params[2],$params[3]); # then round accordingly
1295 # clear a/p after round, since user did not request it
1296 $x->{_a} = undef; $x->{_p} = undef;
1299 $$abr = $ab; $$pbr = $pb;
1305 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1306 # compute power of two numbers, second arg is used as integer
1307 # modifies first argument
1309 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1311 return $x if $x->{sign} =~ /^[+-]inf$/;
1312 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1313 return $x->bone() if $y->is_zero();
1314 return $x if $x->is_one() || $y->is_one();
1316 return $x->_pow2($y,$a,$p,$r) if !$y->is_int(); # non-integer power
1318 my $y1 = $y->as_number(); # make bigint
1320 if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1322 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
1323 return $y1->is_odd() ? $x : $x->babs(1);
1327 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1328 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1332 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1334 $x->{_m}->bpow($y1);
1335 $x->{_e}->bmul($y1);
1336 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1338 if ($y->{sign} eq '-')
1340 # modify $x in place!
1341 my $z = $x->copy(); $x->bzero()->binc();
1342 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
1344 $x->round($a,$p,$r,$y);
1347 ###############################################################################
1348 # rounding functions
1352 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1353 # $n == 0 means round to integer
1354 # expects and returns normalized numbers!
1355 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1357 return $x if $x->modify('bfround');
1359 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1360 return $x if !defined $scale; # no-op
1362 # never round a 0, +-inf, NaN
1365 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1368 return $x if $x->{sign} !~ /^[+-]$/;
1369 # print "MBF bfround $x to scale $scale mode $mode\n";
1371 # don't round if x already has lower precision
1372 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1374 $x->{_p} = $scale; # remember round in any case
1375 $x->{_a} = undef; # and clear A
1378 # print "bfround scale $scale e $x->{_e}\n";
1379 # round right from the '.'
1380 return $x if $x->{_e} >= 0; # nothing to round
1381 $scale = -$scale; # positive for simplicity
1382 my $len = $x->{_m}->length(); # length of mantissa
1383 my $dad = -$x->{_e}; # digits after dot
1384 my $zad = 0; # zeros after dot
1385 $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
1386 #print "scale $scale dad $dad zad $zad len $len\n";
1388 # number bsstr len zad dad
1389 # 0.123 123e-3 3 0 3
1390 # 0.0123 123e-4 3 1 4
1393 # 1.2345 12345e-4 5 0 4
1395 # do not round after/right of the $dad
1396 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1398 # round to zero if rounding inside the $zad, but not for last zero like:
1399 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1400 return $x->bzero() if $scale < $zad;
1401 if ($scale == $zad) # for 0.006, scale -3 and trunc
1407 # adjust round-point to be inside mantissa
1410 $scale = $scale-$zad;
1414 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
1415 $scale = $dbd+$scale;
1418 # print "round to $x->{_m} to $scale\n";
1422 # 123 => 100 means length(123) = 3 - $scale (2) => 1
1424 my $dbt = $x->{_m}->length();
1426 my $dbd = $dbt + $x->{_e};
1427 # should be the same, so treat it as this
1428 $scale = 1 if $scale == 0;
1429 # shortcut if already integer
1430 return $x if $scale == 1 && $dbt <= $dbd;
1431 # maximum digits before dot
1436 # not enough digits before dot, so round to zero
1439 elsif ( $scale == $dbd )
1446 $scale = $dbd - $scale;
1450 # print "using $scale for $x->{_m} with '$mode'\n";
1451 # pass sign to bround for rounding modes '+inf' and '-inf'
1452 $x->{_m}->{sign} = $x->{sign};
1453 $x->{_m}->bround($scale,$mode);
1454 $x->{_m}->{sign} = '+'; # fix sign back
1460 # accuracy: preserve $N digits, and overwrite the rest with 0's
1461 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1463 die ('bround() needs positive accuracy') if ($_[0] || 0) < 0;
1465 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1466 return $x if !defined $scale; # no-op
1468 return $x if $x->modify('bround');
1470 # scale is now either $x->{_a}, $accuracy, or the user parameter
1471 # test whether $x already has lower accuracy, do nothing in this case
1472 # but do round if the accuracy is the same, since a math operation might
1473 # want to round a number with A=5 to 5 digits afterwards again
1474 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1476 # scale < 0 makes no sense
1477 # never round a +-inf, NaN
1478 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
1480 # 1: $scale == 0 => keep all digits
1481 # 2: never round a 0
1482 # 3: if we should keep more digits than the mantissa has, do nothing
1483 if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1485 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1489 # pass sign to bround for '+inf' and '-inf' rounding modes
1490 $x->{_m}->{sign} = $x->{sign};
1491 $x->{_m}->bround($scale,$mode); # round mantissa
1492 $x->{_m}->{sign} = '+'; # fix sign back
1493 # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1494 $x->{_a} = $scale; # remember rounding
1495 $x->{_p} = undef; # and clear P
1496 $x->bnorm(); # del trailing zeros gen. by bround()
1501 # return integer less or equal then $x
1502 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1504 return $x if $x->modify('bfloor');
1506 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1508 # if $x has digits after dot
1509 if ($x->{_e}->{sign} eq '-')
1511 #$x->{_m}->brsft(-$x->{_e},10);
1513 #$x-- if $x->{sign} eq '-';
1515 $x->{_e}->{sign} = '+'; # negate e
1516 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1517 $x->{_e}->bzero(); # trunc/norm
1518 $x->{_m}->binc() if $x->{sign} eq '-'; # decrement if negative
1520 $x->round($a,$p,$r);
1525 # return integer greater or equal then $x
1526 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1528 return $x if $x->modify('bceil');
1529 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1531 # if $x has digits after dot
1532 if ($x->{_e}->{sign} eq '-')
1534 #$x->{_m}->brsft(-$x->{_e},10);
1536 #$x++ if $x->{sign} eq '+';
1538 $x->{_e}->{sign} = '+'; # negate e
1539 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1540 $x->{_e}->bzero(); # trunc/norm
1541 $x->{_m}->binc() if $x->{sign} eq '+'; # decrement if negative
1543 $x->round($a,$p,$r);
1548 # shift right by $y (divide by power of 2)
1549 my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1551 return $x if $x->modify('brsft');
1552 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1554 $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1555 $x->bdiv($n ** $y,$a,$p,$r,$y);
1560 # shift right by $y (divide by power of 2)
1561 my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1563 return $x if $x->modify('brsft');
1564 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1566 $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1567 $x->bmul($n ** $y,$a,$p,$r,$y);
1570 ###############################################################################
1574 # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1579 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
1580 # or falling back to MBI::bxxx()
1581 my $name = $AUTOLOAD;
1583 $name =~ s/.*:://; # split package
1585 if (!method_alias($name))
1589 # delayed load of Carp and avoid recursion
1591 Carp::croak ("Can't call a method without name");
1593 if (!method_hand_up($name))
1595 # delayed load of Carp and avoid recursion
1597 Carp::croak ("Can't call $class\-\>$name, not a valid method");
1599 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1601 return &{"$MBI"."::$name"}(@_);
1603 my $bname = $name; $bname =~ s/^f/b/;
1604 *{$class."::$name"} = \&$bname;
1610 # return a copy of the exponent
1611 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1613 if ($x->{sign} !~ /^[+-]$/)
1615 my $s = $x->{sign}; $s =~ s/^[+-]//;
1616 return $self->new($s); # -inf, +inf => +inf
1618 return $x->{_e}->copy();
1623 # return a copy of the mantissa
1624 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1626 if ($x->{sign} !~ /^[+-]$/)
1628 my $s = $x->{sign}; $s =~ s/^[+]//;
1629 return $self->new($s); # -inf, +inf => +inf
1631 my $m = $x->{_m}->copy(); # faster than going via bstr()
1632 $m->bneg() if $x->{sign} eq '-';
1639 # return a copy of both the exponent and the mantissa
1640 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1642 if ($x->{sign} !~ /^[+-]$/)
1644 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
1645 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
1647 my $m = $x->{_m}->copy(); # faster than going via bstr()
1648 $m->bneg() if $x->{sign} eq '-';
1649 return ($m,$x->{_e}->copy());
1652 ##############################################################################
1653 # private stuff (internal use only)
1658 my $l = scalar @_; my $j = 0; my @a = @_;
1660 for ( my $i = 0; $i < $l ; $i++, $j++)
1662 if ( $_[$i] eq ':constant' )
1664 # this rest causes overlord er load to step in
1665 # print "overload @_\n";
1666 overload::constant float => sub { $self->new(shift); };
1667 splice @a, $j, 1; $j--;
1669 elsif ($_[$i] eq 'upgrade')
1671 # this causes upgrading
1672 $upgrade = $_[$i+1]; # or undef to disable
1673 my $s = 2; $s = 1 if @a-$j < 2; # avoid "can not modify non-existant..."
1674 splice @a, $j, $s; $j -= $s;
1676 elsif ($_[$i] eq 'downgrade')
1678 # this causes downgrading
1679 $downgrade = $_[$i+1]; # or undef to disable
1680 my $s = 2; $s = 1 if @a-$j < 2; # avoid "can not modify non-existant..."
1681 splice @a, $j, $s; $j -= $s;
1683 elsif ($_[$i] eq 'lib')
1685 $lib = $_[$i+1] || ''; # default Calc
1686 my $s = 2; $s = 1 if @a-$j < 2; # avoid "can not modify non-existant..."
1687 splice @a, $j, $s; $j -= $s;
1689 elsif ($_[$i] eq 'with')
1691 $MBI = $_[$i+1] || 'Math::BigInt'; # default Math::BigInt
1692 my $s = 2; $s = 1 if @a-$j < 2; # avoid "can not modify non-existant..."
1693 splice @a, $j, $s; $j -= $s;
1696 my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
1697 my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
1698 $file = File::Spec->catdir (@parts, $file);
1699 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
1700 my $mbilib = eval { Math::BigInt->config()->{lib} };
1701 $lib .= ",$mbilib" if defined $mbilib;
1703 $MBI->import ( lib => $lib, 'objectify' );
1705 # any non :constant stuff is handled by our parent, Exporter
1706 # even if @_ is empty, to give it a chance
1707 $self->SUPER::import(@a); # for subclasses
1708 $self->export_to_level(1,$self,@a); # need this, too
1713 # adjust m and e so that m is smallest possible
1714 # round number according to accuracy and precision settings
1715 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1717 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
1719 # if (!$x->{_m}->is_odd())
1721 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
1724 $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros);
1726 # for something like 0Ey, set y to 1, and -0 => +0
1727 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
1729 # this is to prevent automatically rounding when MBI's globals are set
1730 $x->{_m}->{_f} = MB_NEVER_ROUND;
1731 $x->{_e}->{_f} = MB_NEVER_ROUND;
1732 # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
1733 $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
1734 $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
1735 $x; # MBI bnorm is no-op, so dont call it
1738 ##############################################################################
1739 # internal calculation routines
1743 # return copy as a bigint representation of this BigFloat number
1744 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1746 my $z = $x->{_m}->copy();
1747 if ($x->{_e}->{sign} eq '-') # < 0
1749 $x->{_e}->{sign} = '+'; # flip
1750 $z->brsft($x->{_e},10);
1751 $x->{_e}->{sign} = '-'; # flip back
1753 elsif (!$x->{_e}->is_zero()) # > 0
1755 $z->blsft($x->{_e},10);
1757 $z->{sign} = $x->{sign};
1764 my $class = ref($x) || $x;
1765 $x = $class->new(shift) unless ref($x);
1767 return 1 if $x->{_m}->is_zero();
1768 my $len = $x->{_m}->length();
1769 $len += $x->{_e} if $x->{_e}->sign() eq '+';
1772 my $t = $MBI->bzero();
1773 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1784 Math::BigFloat - Arbitrary size floating point math package
1791 $x = Math::BigFloat->new($str); # defaults to 0
1792 $nan = Math::BigFloat->bnan(); # create a NotANumber
1793 $zero = Math::BigFloat->bzero(); # create a +0
1794 $inf = Math::BigFloat->binf(); # create a +inf
1795 $inf = Math::BigFloat->binf('-'); # create a -inf
1796 $one = Math::BigFloat->bone(); # create a +1
1797 $one = Math::BigFloat->bone('-'); # create a -1
1800 $x->is_zero(); # true if arg is +0
1801 $x->is_nan(); # true if arg is NaN
1802 $x->is_one(); # true if arg is +1
1803 $x->is_one('-'); # true if arg is -1
1804 $x->is_odd(); # true if odd, false for even
1805 $x->is_even(); # true if even, false for odd
1806 $x->is_positive(); # true if >= 0
1807 $x->is_negative(); # true if < 0
1808 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
1810 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
1811 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
1812 $x->sign(); # return the sign, either +,- or NaN
1813 $x->digit($n); # return the nth digit, counting from right
1814 $x->digit(-$n); # return the nth digit, counting from left
1816 # The following all modify their first argument:
1819 $x->bzero(); # set $i to 0
1820 $x->bnan(); # set $i to NaN
1821 $x->bone(); # set $x to +1
1822 $x->bone('-'); # set $x to -1
1823 $x->binf(); # set $x to inf
1824 $x->binf('-'); # set $x to -inf
1826 $x->bneg(); # negation
1827 $x->babs(); # absolute value
1828 $x->bnorm(); # normalize (no-op)
1829 $x->bnot(); # two's complement (bit wise not)
1830 $x->binc(); # increment x by 1
1831 $x->bdec(); # decrement x by 1
1833 $x->badd($y); # addition (add $y to $x)
1834 $x->bsub($y); # subtraction (subtract $y from $x)
1835 $x->bmul($y); # multiplication (multiply $x by $y)
1836 $x->bdiv($y); # divide, set $i to quotient
1837 # return (quo,rem) or quo if scalar
1839 $x->bmod($y); # modulus
1840 $x->bpow($y); # power of arguments (a**b)
1841 $x->blsft($y); # left shift
1842 $x->brsft($y); # right shift
1843 # return (quo,rem) or quo if scalar
1845 $x->blog($base); # logarithm of $x, base defaults to e
1846 # (other bases than e not supported yet)
1848 $x->band($y); # bit-wise and
1849 $x->bior($y); # bit-wise inclusive or
1850 $x->bxor($y); # bit-wise exclusive or
1851 $x->bnot(); # bit-wise not (two's complement)
1853 $x->bsqrt(); # calculate square-root
1854 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
1856 $x->bround($N); # accuracy: preserver $N digits
1857 $x->bfround($N); # precision: round to the $Nth digit
1859 # The following do not modify their arguments:
1860 bgcd(@values); # greatest common divisor
1861 blcm(@values); # lowest common multiplicator
1863 $x->bstr(); # return string
1864 $x->bsstr(); # return string in scientific notation
1866 $x->bfloor(); # return integer less or equal than $x
1867 $x->bceil(); # return integer greater or equal than $x
1869 $x->exponent(); # return exponent as BigInt
1870 $x->mantissa(); # return mantissa as BigInt
1871 $x->parts(); # return (mantissa,exponent) as BigInt
1873 $x->length(); # number of digits (w/o sign and '.')
1874 ($l,$f) = $x->length(); # number of digits, and length of fraction
1878 All operators (inlcuding basic math operations) are overloaded if you
1879 declare your big floating point numbers as
1881 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1883 Operations with overloaded operators preserve the arguments, which is
1884 exactly what you expect.
1886 =head2 Canonical notation
1888 Input to these routines are either BigFloat objects, or strings of the
1889 following four forms:
1903 C</^[+-]\d+E[+-]?\d+$/>
1907 C</^[+-]\d*\.\d+E[+-]?\d+$/>
1911 all with optional leading and trailing zeros and/or spaces. Additonally,
1912 numbers are allowed to have an underscore between any two digits.
1914 Empty strings as well as other illegal numbers results in 'NaN'.
1916 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1917 are always stored in normalized form. On a string, it creates a BigFloat
1922 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1924 The string output will always have leading and trailing zeros stripped and drop
1925 a plus sign. C<bstr()> will give you always the form with a decimal point,
1926 while C<bsstr()> (for scientific) gives you the scientific notation.
1928 Input bstr() bsstr()
1930 ' -123 123 123' '-123123123' '-123123123E0'
1931 '00.0123' '0.0123' '123E-4'
1932 '123.45E-2' '1.2345' '12345E-4'
1933 '10E+3' '10000' '1E4'
1935 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1936 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1937 return either undef, <0, 0 or >0 and are suited for sort.
1939 Actual math is done by using BigInts to represent the mantissa and exponent.
1940 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1941 represent the result when input arguments are not numbers, as well as
1942 the result of dividing by zero.
1944 =head2 C<mantissa()>, C<exponent()> and C<parts()>
1946 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1947 as BigInts such that:
1949 $m = $x->mantissa();
1950 $e = $x->exponent();
1951 $y = $m * ( 10 ** $e );
1952 print "ok\n" if $x == $y;
1954 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
1956 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
1958 Currently the mantissa is reduced as much as possible, favouring higher
1959 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
1960 This might change in the future, so do not depend on it.
1962 =head2 Accuracy vs. Precision
1964 See also: L<Rounding|Rounding>.
1966 Math::BigFloat supports both precision and accuracy. For a full documentation,
1967 examples and tips on these topics please see the large section in
1970 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
1971 a operation consumes all resources, each operation produces no more than
1972 C<Math::BigFloat::precision()> digits.
1974 In case the result of one operation has more precision than specified,
1975 it is rounded. The rounding mode taken is either the default mode, or the one
1976 supplied to the operation after the I<scale>:
1978 $x = Math::BigFloat->new(2);
1979 Math::BigFloat::precision(5); # 5 digits max
1980 $y = $x->copy()->bdiv(3); # will give 0.66666
1981 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1982 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
1983 Math::BigFloat::round_mode('zero');
1984 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1990 =item ffround ( +$scale )
1992 Rounds to the $scale'th place left from the '.', counting from the dot.
1993 The first digit is numbered 1.
1995 =item ffround ( -$scale )
1997 Rounds to the $scale'th place right from the '.', counting from the dot.
2001 Rounds to an integer.
2003 =item fround ( +$scale )
2005 Preserves accuracy to $scale digits from the left (aka significant digits)
2006 and pads the rest with zeros. If the number is between 1 and -1, the
2007 significant digits count from the first non-zero after the '.'
2009 =item fround ( -$scale ) and fround ( 0 )
2011 These are effetively no-ops.
2015 All rounding functions take as a second parameter a rounding mode from one of
2016 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2018 The default rounding mode is 'even'. By using
2019 C<< Math::BigFloat::round_mode($round_mode); >> you can get and set the default
2020 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2021 no longer supported.
2022 The second parameter to the round functions then overrides the default
2025 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
2026 'trunc' as rounding mode to make it equivalent to:
2031 You can override this by passing the desired rounding mode as parameter to
2034 $x = Math::BigFloat->new(2.5);
2035 $y = $x->as_number('odd'); # $y = 3
2041 =head1 Autocreating constants
2043 After C<use Math::BigFloat ':constant'> all the floating point constants
2044 in the given scope are converted to C<Math::BigFloat>. This conversion
2045 happens at compile time.
2049 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2051 prints the value of C<2E-100>. Note that without conversion of
2052 constants the expression 2E-100 will be calculated as normal floating point
2055 Please note that ':constant' does not affect integer constants, nor binary
2056 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2061 Math with the numbers is done (by default) by a module called
2062 Math::BigInt::Calc. This is equivalent to saying:
2064 use Math::BigFloat lib => 'Calc';
2066 You can change this by using:
2068 use Math::BigFloat lib => 'BitVect';
2070 The following would first try to find Math::BigInt::Foo, then
2071 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2073 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2075 Calc.pm uses as internal format an array of elements of some decimal base
2076 (usually 1e7, but this might be differen for some systems) with the least
2077 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2078 significant bit first. Other modules might use even different means of
2079 representing the numbers. See the respective module documentation for further
2082 Please note that Math::BigFloat does B<not> use the denoted library itself,
2083 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2086 use Math::BigInt lib => 'GMP';
2089 you can roll it all into one line:
2091 use Math::BigFloat lib => 'GMP';
2093 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details.
2095 =head2 Using Math::BigInt::Lite
2097 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2100 use Math::BigFloat with => 'Math::BigInt::Lite';
2102 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2103 can combine these if you want. For instance, you may want to use
2104 Math::BigInt objects in your main script, too.
2108 use Math::BigFloat with => 'Math::BigInt::Lite';
2110 Of course, you can combine this with the C<lib> parameter.
2113 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2115 If you want to use Math::BigInt's, too, simple add a Math::BigInt B<before>:
2119 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2121 Notice that the module with the last C<lib> will "win" and thus
2122 it's lib will be used if the lib is available:
2125 use Math::BigInt lib => 'Bar,Baz';
2126 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2128 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2129 words, Math::BigFloat will try to retain previously loaded libs when you
2130 don't specify it one.
2132 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2133 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2134 same as trying the latter load alone, except for the fact that Bar or Baz
2135 might be loaded needlessly in an intermidiate step
2137 The old way still works though:
2140 use Math::BigInt lib => 'Bar,Baz';
2143 But B<examples #3 and #4 are recommended> for usage.
2151 The following does not work yet:
2153 $m = $x->mantissa();
2154 $e = $x->exponent();
2155 $y = $m * ( 10 ** $e );
2156 print "ok\n" if $x == $y;
2160 There is no fmod() function yet.
2168 =item stringify, bstr()
2170 Both stringify and bstr() now drop the leading '+'. The old code would return
2171 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2172 reasoning and details.
2176 The following will probably not do what you expect:
2178 print $c->bdiv(123.456),"\n";
2180 It prints both quotient and reminder since print works in list context. Also,
2181 bdiv() will modify $c, so be carefull. You probably want to use
2183 print $c / 123.456,"\n";
2184 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2188 =item Modifying and =
2192 $x = Math::BigFloat->new(5);
2195 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2196 a second reference to the B<same> object and stores it in $y. Thus anything
2197 that modifies $x will modify $y, and vice versa.
2200 print "$x, $y\n"; # prints '10, 10'
2202 If you want a true copy of $x, use:
2206 See also the documentation in L<overload> regarding C<=>.
2210 C<bpow()> now modifies the first argument, unlike the old code which left
2211 it alone and only returned the result. This is to be consistent with
2212 C<badd()> etc. The first will modify $x, the second one won't:
2214 print bpow($x,$i),"\n"; # modify $x
2215 print $x->bpow($i),"\n"; # ditto
2216 print $x ** $i,"\n"; # leave $x alone
2222 This program is free software; you may redistribute it and/or modify it under
2223 the same terms as Perl itself.
2227 Mark Biggar, overloaded interface by Ilya Zakharevich.
2228 Completely rewritten by Tels http://bloodgate.com in 2001.