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
18 use Math::BigInt qw/objectify/;
19 @ISA = qw( Exporter Math::BigInt);
22 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode/;
23 use vars qw/$upgrade $downgrade/;
24 my $class = "Math::BigFloat";
27 '<=>' => sub { $_[2] ?
28 ref($_[0])->bcmp($_[1],$_[0]) :
29 ref($_[0])->bcmp($_[0],$_[1])},
30 'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
33 ##############################################################################
34 # global constants, flags and accessory
36 use constant MB_NEVER_ROUND => 0x0001;
40 # constant for easier life
43 # class constants, use Class->constant_name() to access
44 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
52 ##############################################################################
53 # the old code had $rnd_mode, so we need to support it, too
56 sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
57 sub FETCH { return $round_mode; }
58 sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
60 BEGIN { tie $rnd_mode, 'Math::BigFloat'; }
62 ##############################################################################
64 # in case we call SUPER::->foo() and this wants to call modify()
65 # sub modify () { 0; }
68 # valid method aliases for AUTOLOAD
69 my %methods = map { $_ => 1 }
70 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
71 fint facmp fcmp fzero fnan finf finc fdec flog ffac
72 fceil ffloor frsft flsft fone flog
74 # valid method's that can be hand-ed up (for AUTOLOAD)
75 my %hand_ups = map { $_ => 1 }
76 qw / is_nan is_inf is_negative is_positive
77 accuracy precision div_scale round_mode fneg fabs babs fnot
78 objectify upgrade downgrade
82 sub method_alias { return exists $methods{$_[0]||''}; }
83 sub method_hand_up { return exists $hand_ups{$_[0]||''}; }
86 ##############################################################################
91 # create a new BigFloat object from a string or another bigfloat object.
94 # sign => sign (+/-), or "NaN"
96 my ($class,$wanted,@r) = @_;
98 # avoid numify-calls by not using || on $wanted!
99 return $class->bzero() if !defined $wanted; # default to 0
100 return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
102 my $self = {}; bless $self, $class;
103 # shortcut for bigints and its subclasses
104 if ((ref($wanted)) && (ref($wanted) ne $class))
106 $self->{_m} = $wanted->as_number(); # get us a bigint copy
107 $self->{_e} = Math::BigInt->bzero();
109 $self->{sign} = $wanted->sign();
110 return $self->bnorm();
113 # handle '+inf', '-inf' first
114 if ($wanted =~ /^[+-]?inf$/)
116 return $downgrade->new($wanted) if $downgrade;
118 $self->{_e} = Math::BigInt->bzero();
119 $self->{_m} = Math::BigInt->bzero();
120 $self->{sign} = $wanted;
121 $self->{sign} = '+inf' if $self->{sign} eq 'inf';
122 return $self->bnorm();
124 #print "new string '$wanted'\n";
125 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
128 die "$wanted is not a number initialized to $class" if !$NaNOK;
130 return $downgrade->bnan() if $downgrade;
132 $self->{_e} = Math::BigInt->bzero();
133 $self->{_m} = Math::BigInt->bzero();
134 $self->{sign} = $nan;
138 # make integer from mantissa by adjusting exp, then convert to bigint
139 # undef,undef to signal MBI that we don't need no bloody rounding
140 $self->{_e} = Math::BigInt->new("$$es$$ev",undef,undef); # exponent
141 $self->{_m} = Math::BigInt->new("$$miv$$mfv",undef,undef); # create mant.
142 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
143 $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;
144 $self->{sign} = $$mis;
146 # if downgrade, inf, NaN or integers go down
148 if ($downgrade && $self->{_e}->{sign} eq '+')
150 # print "downgrading $$miv$$mfv"."E$$es$$ev";
151 if ($self->{_e}->is_zero())
153 $self->{_m}->{sign} = $$mis; # negative if wanted
154 return $downgrade->new($self->{_m});
156 return $downgrade->new("$$mis$$miv$$mfv"."E$$es$$ev");
158 # print "mbf new $self->{sign} $self->{_m} e $self->{_e} ",ref($self),"\n";
159 $self->bnorm()->round(@r); # first normalize, then round
164 # used by parent class bone() to initialize number to 1
166 $self->{_m} = Math::BigInt->bzero();
167 $self->{_e} = Math::BigInt->bzero();
172 # used by parent class bone() to initialize number to 1
174 $self->{_m} = Math::BigInt->bzero();
175 $self->{_e} = Math::BigInt->bzero();
180 # used by parent class bone() to initialize number to 1
182 $self->{_m} = Math::BigInt->bone();
183 $self->{_e} = Math::BigInt->bzero();
188 # used by parent class bone() to initialize number to 1
190 $self->{_m} = Math::BigInt->bzero();
191 $self->{_e} = Math::BigInt->bone();
196 my ($self,$class) = @_;
197 return if $class eq 'Math::BigInt'; # we aren't
198 return UNIVERSAL::isa($self,$class);
201 ##############################################################################
202 # string conversation
206 # (ref to BFLOAT or num_str ) return num_str
207 # Convert number from internal format to (non-scientific) string format.
208 # internal format is always normalized (no leading zeros, "-0" => "+0")
209 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
210 #my $x = shift; my $class = ref($x) || $x;
211 #$x = $class->new(shift) unless ref($x);
213 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
214 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
215 if ($x->{sign} !~ /^[+-]$/)
217 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
221 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
223 my $not_zero = ! $x->is_zero();
226 $es = $x->{_m}->bstr();
227 $len = CORE::length($es);
228 if (!$x->{_e}->is_zero())
230 if ($x->{_e}->sign() eq '-')
233 if ($x->{_e} <= -$len)
235 # print "style: 0.xxxx\n";
236 my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
237 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
241 # print "insert '.' at $x->{_e} in '$es'\n";
242 substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
248 $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
252 $es = $x->{sign}.$es if $x->{sign} eq '-';
253 # if set accuracy or precision, pad with zeros
254 if ((defined $x->{_a}) && ($not_zero))
256 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
257 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
258 $zeros = $x->{_a} - $len if $cad != $len;
259 $es .= $dot.'0' x $zeros if $zeros > 0;
261 elsif ($x->{_p} || 0 < 0)
263 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
264 my $zeros = -$x->{_p} + $cad;
265 $es .= $dot.'0' x $zeros if $zeros > 0;
272 # (ref to BFLOAT or num_str ) return num_str
273 # Convert number from internal format to scientific string format.
274 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
275 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
276 #my $x = shift; my $class = ref($x) || $x;
277 #$x = $class->new(shift) unless ref($x);
279 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
280 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
281 if ($x->{sign} !~ /^[+-]$/)
283 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
286 my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
288 return $x->{_m}->bstr().$sep.$x->{_e}->bstr();
293 # Make a number from a BigFloat object
294 # simple return string and let Perl's atoi()/atof() handle the rest
295 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
299 ##############################################################################
300 # public stuff (usually prefixed with "b")
303 # todo: this must be overwritten and return NaN for non-integer values
304 # band(), bior(), bxor(), too
307 # $class->SUPER::bnot($class,@_);
312 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
313 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
314 my ($self,$x,$y) = objectify(2,@_);
316 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
318 # handle +-inf and NaN
319 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
320 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
321 return +1 if $x->{sign} eq '+inf';
322 return -1 if $x->{sign} eq '-inf';
323 return -1 if $y->{sign} eq '+inf';
327 # check sign for speed first
328 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
329 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
332 my $xz = $x->is_zero();
333 my $yz = $y->is_zero();
334 return 0 if $xz && $yz; # 0 <=> 0
335 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
336 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
338 # adjust so that exponents are equal
339 my $lxm = $x->{_m}->length();
340 my $lym = $y->{_m}->length();
341 # the numify somewhat limits our length, but makes it much faster
342 my $lx = $lxm + $x->{_e}->numify();
343 my $ly = $lym + $y->{_e}->numify();
344 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
345 return $l <=> 0 if $l != 0;
347 # lengths (corrected by exponent) are equal
348 # so make mantissa equal length by padding with zero (shift left)
349 my $diff = $lxm - $lym;
350 my $xm = $x->{_m}; # not yet copy it
354 $ym = $y->{_m}->copy()->blsft($diff,10);
358 $xm = $x->{_m}->copy()->blsft(-$diff,10);
360 my $rc = $xm->bacmp($ym);
361 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
367 # Compares 2 values, ignoring their signs.
368 # Returns one of undef, <0, =0, >0. (suitable for sort)
369 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
370 my ($self,$x,$y) = objectify(2,@_);
372 # handle +-inf and NaN's
373 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
375 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
376 return 0 if ($x->is_inf() && $y->is_inf());
377 return 1 if ($x->is_inf() && !$y->is_inf());
382 my $xz = $x->is_zero();
383 my $yz = $y->is_zero();
384 return 0 if $xz && $yz; # 0 <=> 0
385 return -1 if $xz && !$yz; # 0 <=> +y
386 return 1 if $yz && !$xz; # +x <=> 0
388 # adjust so that exponents are equal
389 my $lxm = $x->{_m}->length();
390 my $lym = $y->{_m}->length();
391 # the numify somewhat limits our length, but makes it much faster
392 my $lx = $lxm + $x->{_e}->numify();
393 my $ly = $lym + $y->{_e}->numify();
395 return $l <=> 0 if $l != 0;
397 # lengths (corrected by exponent) are equal
398 # so make mantissa equal-length by padding with zero (shift left)
399 my $diff = $lxm - $lym;
400 my $xm = $x->{_m}; # not yet copy it
404 $ym = $y->{_m}->copy()->blsft($diff,10);
408 $xm = $x->{_m}->copy()->blsft(-$diff,10);
410 $xm->bacmp($ym) <=> 0;
415 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
416 # return result as BFLOAT
417 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
419 #print "mbf badd $x $y\n";
420 # inf and NaN handling
421 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
424 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
426 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
428 # +inf++inf or -inf+-inf => same, rest is NaN
429 return $x if $x->{sign} eq $y->{sign};
432 # +-inf + something => +inf
433 # something +-inf => +-inf
434 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
438 # speed: no add for 0+y or x+0
439 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
440 if ($x->is_zero()) # 0+y
442 # make copy, clobbering up x (modify in place!)
443 $x->{_e} = $y->{_e}->copy();
444 $x->{_m} = $y->{_m}->copy();
445 $x->{sign} = $y->{sign} || $nan;
446 return $x->round($a,$p,$r,$y);
449 # take lower of the two e's and adapt m1 to it to match m2
451 $e = Math::BigInt::bzero() if !defined $e; # if no BFLOAT ?
452 $e = $e->copy(); # make copy (didn't do it yet)
454 my $add = $y->{_m}->copy();
456 if ($e->{sign} eq '-') # < 0
458 my $e1 = $e->copy()->babs();
459 #$x->{_m} *= (10 ** $e1);
460 $x->{_m}->blsft($e1,10);
461 $x->{_e} += $e; # need the sign of e
464 elsif (!$e->is_zero()) # > 0
469 # else: both e are the same, so just leave them
470 $x->{_m}->{sign} = $x->{sign}; # fiddle with signs
471 $add->{sign} = $y->{sign};
472 $x->{_m} += $add; # finally do add/sub
473 $x->{sign} = $x->{_m}->{sign}; # re-adjust signs
474 $x->{_m}->{sign} = '+'; # mantissa always positiv
475 # delete trailing zeros, then round
476 return $x->bnorm()->round($a,$p,$r,$y);
481 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
482 # subtract second arg from first, modify first
483 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
485 if ($y->is_zero()) # still round for not adding zero
487 return $x->round($a,$p,$r);
490 $y->{sign} =~ tr/+\-/-+/; # does nothing for NaN
491 $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
492 $y->{sign} =~ tr/+\-/-+/; # refix $y (does nothing for NaN)
493 $x; # already rounded by badd()
498 # increment arg by one
499 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
501 if ($x->{_e}->sign() eq '-')
503 return $x->badd($self->bone(),$a,$p,$r); # digits after dot
506 if (!$x->{_e}->is_zero())
508 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
512 if ($x->{sign} eq '+')
515 return $x->bnorm()->bround($a,$p,$r);
517 elsif ($x->{sign} eq '-')
520 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
521 return $x->bnorm()->bround($a,$p,$r);
523 # inf, nan handling etc
524 $x->badd($self->__one(),$a,$p,$r); # does round
529 # decrement arg by one
530 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
532 if ($x->{_e}->sign() eq '-')
534 return $x->badd($self->bone('-'),$a,$p,$r); # digits after dot
537 if (!$x->{_e}->is_zero())
539 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
543 my $zero = $x->is_zero();
545 if (($x->{sign} eq '-') || $zero)
548 $x->{sign} = '-' if $zero; # 0 => 1 => -1
549 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
550 return $x->bnorm()->round($a,$p,$r);
553 elsif ($x->{sign} eq '+')
556 return $x->bnorm()->round($a,$p,$r);
558 # inf, nan handling etc
559 $x->badd($self->bone('-'),$a,$p,$r); # does round
564 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(2,@_);
566 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
570 # Taylor: | u 1 u^3 1 u^5 |
571 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
572 # |_ v 3 v^3 5 v^5 _|
574 # This takes much more steps to calculate the result:
577 # Taylor: | u 1 u^2 1 u^3 |
578 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
579 # |_ x 2 x^2 3 x^3 _|
581 # we need to limit the accuracy to protect against overflow
584 my @params = $x->_find_round_parameters($a,$p,$r);
586 # no rounding at all, so must use fallback
587 if (scalar @params == 1)
589 # simulate old behaviour
590 $params[1] = $self->div_scale(); # and round to it as accuracy
591 $scale = $params[1]+4; # at least four more for proper round
592 $params[3] = $r; # round mode by caller or undef
593 $fallback = 1; # to clear a/p afterwards
597 # the 4 below is empirical, and there might be cases where it is not
599 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
602 return $x->bzero(@params) if $x->is_one();
603 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
604 #return $x->bone('+',@params) if $x->bcmp($base) == 0;
606 # when user set globals, they would interfere with our calculation, so
607 # disable then and later re-enable them
609 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
610 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
611 # we also need to disable any set A or P on $x (_find_round_parameters took
612 # them already into account), since these would interfere, too
613 delete $x->{_a}; delete $x->{_p};
614 # need to disable $upgrade in BigInt, to avoid deep recursion
615 local $Math::BigInt::upgrade = undef;
617 my ($case,$limit,$v,$u,$below,$factor,$two,$next,$over,$f);
620 #if ($x <= Math::BigFloat->new("0.5"))
623 # print "case $case $x < 0.5\n";
624 $v = $x->copy(); $v->binc(); # v = x+1
625 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
626 $x->bdiv($v,$scale); # first term: u/v
629 $u *= $u; $v *= $v; # u^2, v^2
630 $below->bmul($v); # u^3, v^3
632 $factor = $self->new(3); $f = $self->new(2);
637 # print "case 1 $x > 0.5\n";
638 # $v = $x->copy(); # v = x
639 # $u = $x->copy(); $u->bdec(); # u = x-1;
640 # $x->bdec(); $x->bdiv($v,$scale); # first term: x-1/x
641 # $below = $v->copy();
642 # $over = $u->copy();
643 # $below->bmul($v); # u^2, v^2
645 # $factor = $self->new(2); $f = $self->bone();
647 $limit = $self->new("1E-". ($scale-1));
651 # we calculate the next term, and add it to the last
652 # when the next term is below our limit, it won't affect the outcome
653 # anymore, so we stop
654 $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
655 last if $next->bcmp($limit) <= 0;
657 # print "step $steps $x\n";
658 # calculate things for the next term
659 $over *= $u; $below *= $v; $factor->badd($f);
662 $x->bmul(2) if $case == 0;
663 #print "took $steps steps\n";
665 # shortcut to not run trough _find_round_parameters again
666 if (defined $params[1])
668 $x->bround($params[1],$params[3]); # then round accordingly
672 $x->bfround($params[2],$params[3]); # then round accordingly
676 # clear a/p after round, since user did not request it
677 $x->{_a} = undef; $x->{_p} = undef;
680 $$abr = $ab; $$pbr = $pb;
687 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
688 # does not modify arguments, but returns new object
689 # Lowest Common Multiplicator
691 my ($self,@arg) = objectify(0,@_);
692 my $x = $self->new(shift @arg);
693 while (@arg) { $x = _lcm($x,shift @arg); }
699 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
700 # does not modify arguments, but returns new object
701 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
703 my ($self,@arg) = objectify(0,@_);
704 my $x = $self->new(shift @arg);
705 while (@arg) { $x = _gcd($x,shift @arg); }
709 ###############################################################################
710 # is_foo methods (is_negative, is_positive are inherited from BigInt)
714 # return true if arg (BFLOAT or num_str) is an integer
715 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
717 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
718 $x->{_e}->{sign} eq '+'; # 1e-1 => no integer
724 # return true if arg (BFLOAT or num_str) is zero
725 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
727 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
733 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
734 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
736 my $sign = shift || ''; $sign = '+' if $sign ne '-';
738 if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
744 # return true if arg (BFLOAT or num_str) is odd or false if even
745 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
747 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
748 ($x->{_e}->is_zero() && $x->{_m}->is_odd());
754 # return true if arg (BINT or num_str) is even or false if odd
755 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
757 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
758 return 1 if ($x->{_e}->{sign} eq '+' # 123.45 is never
759 && $x->{_m}->is_even()); # but 1200 is
765 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
766 # (BINT or num_str, BINT or num_str) return BINT
767 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
769 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
772 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
774 return $x->bnan() if $x->is_zero() || $y->is_zero();
775 # result will always be +-inf:
776 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
777 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
778 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
779 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
780 return $x->binf('-');
783 return $x->bzero() if $x->is_zero() || $y->is_zero();
785 # aEb * cEd = (a*c)E(b+d)
786 $x->{_m}->bmul($y->{_m});
787 $x->{_e}->badd($y->{_e});
789 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
790 return $x->bnorm()->round($a,$p,$r,$y);
795 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
796 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
797 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
799 return $self->_div_inf($x,$y)
800 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
802 # x== 0 # also: or y == 1 or y == -1
803 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
806 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
808 # we need to limit the accuracy to protect against overflow
811 my @params = $x->_find_round_parameters($a,$p,$r,$y);
813 # no rounding at all, so must use fallback
814 if (scalar @params == 1)
816 # simulate old behaviour
817 $params[1] = $self->div_scale(); # and round to it as accuracy
818 $scale = $params[1]+4; # at least four more for proper round
819 $params[3] = $r; # round mode by caller or undef
820 $fallback = 1; # to clear a/p afterwards
824 # the 4 below is empirical, and there might be cases where it is not
826 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
828 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
829 $scale = $lx if $lx > $scale;
830 $scale = $ly if $ly > $scale;
831 my $diff = $ly - $lx;
832 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
834 # make copy of $x in case of list context for later reminder calculation
836 if (wantarray && !$y->is_one())
841 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
843 # check for / +-1 ( +/- 1E0)
846 # promote BigInts and it's subclasses (except when already a BigFloat)
847 $y = $self->new($y) unless $y->isa('Math::BigFloat');
849 #print "bdiv $y ",ref($y),"\n";
850 # need to disable $upgrade in BigInt, to avoid deep recursion
851 local $Math::BigInt::upgrade = undef; # should be parent class vs MBI
853 # calculate the result to $scale digits and then round it
854 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
855 $x->{_m}->blsft($scale,10);
856 $x->{_m}->bdiv( $y->{_m} ); # a/c
857 $x->{_e}->bsub( $y->{_e} ); # b-d
858 $x->{_e}->bsub($scale); # correct for 10**scale
859 $x->bnorm(); # remove trailing 0's
862 # shortcut to not run trough _find_round_parameters again
863 if (defined $params[1])
865 $x->bround($params[1],$params[3]); # then round accordingly
869 $x->bfround($params[2],$params[3]); # then round accordingly
873 # clear a/p after round, since user did not request it
874 $x->{_a} = undef; $x->{_p} = undef;
881 $rem->bmod($y,$params[1],$params[2],$params[3]); # copy already done
885 $rem = $self->bzero();
889 # clear a/p after round, since user did not request it
890 $rem->{_a} = undef; $rem->{_p} = undef;
899 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
900 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
902 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
904 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
905 return $re->round($a,$p,$r,$y);
907 return $x->bnan() if $x->is_zero() && $y->is_zero();
908 return $x if $y->is_zero();
909 return $x->bnan() if $x->is_nan() || $y->is_nan();
910 return $x->bzero() if $y->is_one() || $x->is_zero();
912 # inf handling is missing here
914 my $cmp = $x->bacmp($y); # equal or $x < $y?
915 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
917 # only $y of the operands negative?
918 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
920 $x->{sign} = $y->{sign}; # calc sign first
921 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
923 my $ym = $y->{_m}->copy();
926 $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
928 # if $y has digits after dot
929 my $shifty = 0; # correct _e of $x by this
930 if ($y->{_e}->{sign} eq '-') # has digits after dot
932 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
933 $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
934 $x->blsft($shifty,10); # 123 => 1230, $y->{_m} is already 25
936 # $ym is now mantissa of $y based on exponent 0
938 my $shiftx = 0; # correct _e of $x by this
939 if ($x->{_e}->{sign} eq '-') # has digits after dot
941 # 123.4 % 20 => 1234 % 200
942 $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
943 $ym->blsft($shiftx,10);
945 # 123e1 % 20 => 1230 % 20
946 if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
948 $x->{_m}->blsft($x->{_e},10);
950 $x->{_e} = Math::BigInt->bzero() unless $x->{_e}->is_zero();
952 $x->{_e}->bsub($shiftx) if $shiftx != 0;
953 $x->{_e}->bsub($shifty) if $shifty != 0;
955 # now mantissas are equalized, exponent of $x is adjusted, so calc result
956 # $ym->{sign} = '-' if $neg; # bmod() will make the correction for us
960 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
963 if ($neg != 0) # one of them negative => correct in place
968 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
972 $x->round($a,$p,$r,$y); # round and return
977 # calculate square root; this should probably
978 # use a different test to see whether the accuracy we want is...
979 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
981 return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
982 return $x if $x->{sign} eq '+inf'; # +inf
983 return $x if $x->is_zero() || $x->is_one();
985 # we need to limit the accuracy to protect against overflow
988 my @params = $x->_find_round_parameters($a,$p,$r);
990 # no rounding at all, so must use fallback
991 if (scalar @params == 1)
993 # simulate old behaviour
994 $params[1] = $self->div_scale(); # and round to it as accuracy
995 $scale = $params[1]+4; # at least four more for proper round
996 $params[3] = $r; # round mode by caller or undef
997 $fallback = 1; # to clear a/p afterwards
1001 # the 4 below is empirical, and there might be cases where it is not
1003 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1006 # when user set globals, they would interfere with our calculation, so
1007 # disable them and later re-enable them
1009 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1010 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1011 # we also need to disable any set A or P on $x (_find_round_parameters took
1012 # them already into account), since these would interfere, too
1013 delete $x->{_a}; delete $x->{_p};
1014 # need to disable $upgrade in BigInt, to avoid deep recursion
1015 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1017 my $xas = $x->as_number();
1018 my $gs = $xas->copy()->bsqrt(); # some guess
1020 # print "guess $gs\n";
1021 if (($x->{_e}->{sign} ne '-') # guess can't be accurate if there are
1022 # digits after the dot
1023 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
1026 $x->{_m} = $gs; $x->{_e} = Math::BigInt->bzero(); $x->bnorm();
1027 # shortcut to not run trough _find_round_parameters again
1028 if (defined $params[1])
1030 $x->bround($params[1],$params[3]); # then round accordingly
1034 $x->bfround($params[2],$params[3]); # then round accordingly
1038 # clear a/p after round, since user did not request it
1039 $x->{_a} = undef; $x->{_p} = undef;
1041 # re-enable A and P, upgrade is taken care of by "local"
1042 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1045 $gs = $self->new( $gs ); # BigInt to BigFloat
1047 my $lx = $x->{_m}->length();
1048 $scale = $lx if $scale < $lx;
1049 my $e = $self->new("1E-$scale"); # make test variable
1052 my $two = $self->new(2);
1054 # promote BigInts and it's subclasses (except when already a BigFloat)
1055 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1058 while ($diff->bacmp($e) >= 0)
1060 $rem = $y->copy()->bdiv($gs,$scale);
1061 $rem = $y->copy()->bdiv($gs,$scale)->badd($gs)->bdiv($two,$scale);
1062 $diff = $rem->copy()->bsub($gs);
1065 # copy over to modify $x
1066 $x->{_m} = $rem->{_m}; $x->{_e} = $rem->{_e};
1068 # shortcut to not run trough _find_round_parameters again
1069 if (defined $params[1])
1071 $x->bround($params[1],$params[3]); # then round accordingly
1075 $x->bfround($params[2],$params[3]); # then round accordingly
1079 # clear a/p after round, since user did not request it
1080 $x->{_a} = undef; $x->{_p} = undef;
1083 $$abr = $ab; $$pbr = $pb;
1089 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1090 # compute factorial numbers
1091 # modifies first argument
1092 my ($self,$x,@r) = objectify(1,@_);
1095 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1096 ($x->{_e}->{sign} ne '+')); # digits after dot?
1098 return $x->bone(@r) if $x->is_zero() || $x->is_one(); # 0 or 1 => 1
1100 # use BigInt's bfac() for faster calc
1101 $x->{_m}->blsft($x->{_e},10); # un-norm m
1102 $x->{_e}->bzero(); # norm $x again
1103 $x->{_m}->bfac(); # factorial
1104 $x->bnorm()->round(@r);
1109 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1110 my ($x,$y,$a,$p,$r) = @_;
1113 # if $y == 0.5, it is sqrt($x)
1114 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1118 # Taylor: | u u^2 u^3 |
1119 # x ** y = 1 + | --- + --- + * ----- + ... |
1122 # we need to limit the accuracy to protect against overflow
1125 my @params = $x->_find_round_parameters($a,$p,$r);
1127 # no rounding at all, so must use fallback
1128 if (scalar @params == 1)
1130 # simulate old behaviour
1131 $params[1] = $self->div_scale(); # and round to it as accuracy
1132 $scale = $params[1]+4; # at least four more for proper round
1133 $params[3] = $r; # round mode by caller or undef
1134 $fallback = 1; # to clear a/p afterwards
1138 # the 4 below is empirical, and there might be cases where it is not
1140 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1143 # when user set globals, they would interfere with our calculation, so
1144 # disable then and later re-enable them
1146 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1147 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1148 # we also need to disable any set A or P on $x (_find_round_parameters took
1149 # them already into account), since these would interfere, too
1150 delete $x->{_a}; delete $x->{_p};
1151 # need to disable $upgrade in BigInt, to avoid deep recursion
1152 local $Math::BigInt::upgrade = undef;
1154 my ($limit,$v,$u,$below,$factor,$next,$over);
1156 $u = $x->copy()->blog($scale)->bmul($y);
1157 $v = $self->bone(); # 1
1158 $factor = $self->new(2); # 2
1159 $x->bone(); # first term: 1
1161 $below = $v->copy();
1164 $limit = $self->new("1E-". ($scale-1));
1168 # we calculate the next term, and add it to the last
1169 # when the next term is below our limit, it won't affect the outcome
1170 # anymore, so we stop
1171 $next = $over->copy()->bdiv($below,$scale);
1172 last if $next->bcmp($limit) <= 0;
1175 # calculate things for the next term
1176 $over *= $u; $below *= $factor; $factor->binc();
1180 # shortcut to not run trough _find_round_parameters again
1181 if (defined $params[1])
1183 $x->bround($params[1],$params[3]); # then round accordingly
1187 $x->bfround($params[2],$params[3]); # then round accordingly
1191 # clear a/p after round, since user did not request it
1192 $x->{_a} = undef; $x->{_p} = undef;
1195 $$abr = $ab; $$pbr = $pb;
1201 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1202 # compute power of two numbers, second arg is used as integer
1203 # modifies first argument
1205 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1207 return $x if $x->{sign} =~ /^[+-]inf$/;
1208 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1209 return $x->bone() if $y->is_zero();
1210 return $x if $x->is_one() || $y->is_one();
1212 return $x->_pow($y,$a,$p,$r) if !$y->is_int(); # non-integer power
1214 my $y1 = $y->as_number(); # make bigint
1216 if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1218 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
1219 return $y1->is_odd() ? $x : $x->babs(1);
1223 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1224 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1228 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1230 $x->{_m}->bpow($y1);
1231 $x->{_e}->bmul($y1);
1232 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1234 if ($y->{sign} eq '-')
1236 # modify $x in place!
1237 my $z = $x->copy(); $x->bzero()->binc();
1238 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
1240 $x->round($a,$p,$r,$y);
1243 ###############################################################################
1244 # rounding functions
1248 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1249 # $n == 0 means round to integer
1250 # expects and returns normalized numbers!
1251 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1253 return $x if $x->modify('bfround');
1255 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1256 return $x if !defined $scale; # no-op
1258 # never round a 0, +-inf, NaN
1261 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1264 return $x if $x->{sign} !~ /^[+-]$/;
1265 # print "MBF bfround $x to scale $scale mode $mode\n";
1267 # don't round if x already has lower precision
1268 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1270 $x->{_p} = $scale; # remember round in any case
1271 $x->{_a} = undef; # and clear A
1274 # print "bfround scale $scale e $x->{_e}\n";
1275 # round right from the '.'
1276 return $x if $x->{_e} >= 0; # nothing to round
1277 $scale = -$scale; # positive for simplicity
1278 my $len = $x->{_m}->length(); # length of mantissa
1279 my $dad = -$x->{_e}; # digits after dot
1280 my $zad = 0; # zeros after dot
1281 $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
1282 #print "scale $scale dad $dad zad $zad len $len\n";
1284 # number bsstr len zad dad
1285 # 0.123 123e-3 3 0 3
1286 # 0.0123 123e-4 3 1 4
1289 # 1.2345 12345e-4 5 0 4
1291 # do not round after/right of the $dad
1292 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1294 # round to zero if rounding inside the $zad, but not for last zero like:
1295 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1296 return $x->bzero() if $scale < $zad;
1297 if ($scale == $zad) # for 0.006, scale -3 and trunc
1303 # adjust round-point to be inside mantissa
1306 $scale = $scale-$zad;
1310 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
1311 $scale = $dbd+$scale;
1314 # print "round to $x->{_m} to $scale\n";
1318 # 123 => 100 means length(123) = 3 - $scale (2) => 1
1320 my $dbt = $x->{_m}->length();
1322 my $dbd = $dbt + $x->{_e};
1323 # should be the same, so treat it as this
1324 $scale = 1 if $scale == 0;
1325 # shortcut if already integer
1326 return $x if $scale == 1 && $dbt <= $dbd;
1327 # maximum digits before dot
1332 # not enough digits before dot, so round to zero
1335 elsif ( $scale == $dbd )
1342 $scale = $dbd - $scale;
1346 # print "using $scale for $x->{_m} with '$mode'\n";
1347 # pass sign to bround for rounding modes '+inf' and '-inf'
1348 $x->{_m}->{sign} = $x->{sign};
1349 $x->{_m}->bround($scale,$mode);
1350 $x->{_m}->{sign} = '+'; # fix sign back
1356 # accuracy: preserve $N digits, and overwrite the rest with 0's
1357 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1359 die ('bround() needs positive accuracy') if ($_[0] || 0) < 0;
1361 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1362 return $x if !defined $scale; # no-op
1364 return $x if $x->modify('bround');
1366 # scale is now either $x->{_a}, $accuracy, or the user parameter
1367 # test whether $x already has lower accuracy, do nothing in this case
1368 # but do round if the accuracy is the same, since a math operation might
1369 # want to round a number with A=5 to 5 digits afterwards again
1370 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1372 # scale < 0 makes no sense
1373 # never round a +-inf, NaN
1374 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
1376 # 1: $scale == 0 => keep all digits
1377 # 2: never round a 0
1378 # 3: if we should keep more digits than the mantissa has, do nothing
1379 if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1381 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1385 # pass sign to bround for '+inf' and '-inf' rounding modes
1386 $x->{_m}->{sign} = $x->{sign};
1387 $x->{_m}->bround($scale,$mode); # round mantissa
1388 $x->{_m}->{sign} = '+'; # fix sign back
1389 # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1390 $x->{_a} = $scale; # remember rounding
1391 $x->{_p} = undef; # and clear P
1392 $x->bnorm(); # del trailing zeros gen. by bround()
1397 # return integer less or equal then $x
1398 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1400 return $x if $x->modify('bfloor');
1402 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1404 # if $x has digits after dot
1405 if ($x->{_e}->{sign} eq '-')
1407 #$x->{_m}->brsft(-$x->{_e},10);
1409 #$x-- if $x->{sign} eq '-';
1411 $x->{_e}->{sign} = '+'; # negate e
1412 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1413 $x->{_e}->bzero(); # trunc/norm
1414 $x->{_m}->binc() if $x->{sign} eq '-'; # decrement if negative
1416 $x->round($a,$p,$r);
1421 # return integer greater or equal then $x
1422 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1424 return $x if $x->modify('bceil');
1425 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1427 # if $x has digits after dot
1428 if ($x->{_e}->{sign} eq '-')
1430 #$x->{_m}->brsft(-$x->{_e},10);
1432 #$x++ if $x->{sign} eq '+';
1434 $x->{_e}->{sign} = '+'; # negate e
1435 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1436 $x->{_e}->bzero(); # trunc/norm
1437 $x->{_m}->binc() if $x->{sign} eq '+'; # decrement if negative
1439 $x->round($a,$p,$r);
1444 # shift right by $y (divide by power of 2)
1445 my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1447 return $x if $x->modify('brsft');
1448 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1450 $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1451 $x->bdiv($n ** $y,$a,$p,$r,$y);
1456 # shift right by $y (divide by power of 2)
1457 my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1459 return $x if $x->modify('brsft');
1460 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1462 $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1463 $x->bmul($n ** $y,$a,$p,$r,$y);
1466 ###############################################################################
1470 # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1475 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
1476 # or falling back to MBI::bxxx()
1477 my $name = $AUTOLOAD;
1479 $name =~ s/.*:://; # split package
1481 if (!method_alias($name))
1485 # delayed load of Carp and avoid recursion
1487 Carp::croak ("Can't call a method without name");
1489 if (!method_hand_up($name))
1491 # delayed load of Carp and avoid recursion
1493 Carp::croak ("Can't call $class\-\>$name, not a valid method");
1495 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1497 return &{'Math::BigInt'."::$name"}(@_);
1499 my $bname = $name; $bname =~ s/^f/b/;
1500 *{$class."::$name"} = \&$bname;
1506 # return a copy of the exponent
1507 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1509 if ($x->{sign} !~ /^[+-]$/)
1511 my $s = $x->{sign}; $s =~ s/^[+-]//;
1512 return $self->new($s); # -inf, +inf => +inf
1514 return $x->{_e}->copy();
1519 # return a copy of the mantissa
1520 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1522 if ($x->{sign} !~ /^[+-]$/)
1524 my $s = $x->{sign}; $s =~ s/^[+]//;
1525 return $self->new($s); # -inf, +inf => +inf
1527 my $m = $x->{_m}->copy(); # faster than going via bstr()
1528 $m->bneg() if $x->{sign} eq '-';
1535 # return a copy of both the exponent and the mantissa
1536 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1538 if ($x->{sign} !~ /^[+-]$/)
1540 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
1541 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
1543 my $m = $x->{_m}->copy(); # faster than going via bstr()
1544 $m->bneg() if $x->{sign} eq '-';
1545 return ($m,$x->{_e}->copy());
1548 ##############################################################################
1549 # private stuff (internal use only)
1554 my $l = scalar @_; my $j = 0; my @a = @_;
1555 for ( my $i = 0; $i < $l ; $i++, $j++)
1557 if ( $_[$i] eq ':constant' )
1559 # this rest causes overlord er load to step in
1560 # print "overload @_\n";
1561 overload::constant float => sub { $self->new(shift); };
1562 splice @a, $j, 1; $j--;
1564 elsif ($_[$i] eq 'upgrade')
1566 # this causes upgrading
1567 $upgrade = $_[$i+1]; # or undef to disable
1568 my $s = 2; $s = 1 if @a-$j < 2; # avoid "can not modify non-existant..."
1569 splice @a, $j, $s; $j -= $s;
1571 elsif ($_[$i] eq 'downgrade')
1573 # this causes downgrading
1574 $downgrade = $_[$i+1]; # or undef to disable
1575 my $s = 2; $s = 1 if @a-$j < 2; # avoid "can not modify non-existant..."
1576 splice @a, $j, $s; $j -= $s;
1579 # any non :constant stuff is handled by our parent, Exporter
1580 # even if @_ is empty, to give it a chance
1581 $self->SUPER::import(@a); # for subclasses
1582 $self->export_to_level(1,$self,@a); # need this, too
1587 # adjust m and e so that m is smallest possible
1588 # round number according to accuracy and precision settings
1589 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1591 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
1593 # if (!$x->{_m}->is_odd())
1595 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
1598 $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros);
1600 # for something like 0Ey, set y to 1, and -0 => +0
1601 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
1603 # this is to prevent automatically rounding when MBI's globals are set
1604 $x->{_m}->{_f} = MB_NEVER_ROUND;
1605 $x->{_e}->{_f} = MB_NEVER_ROUND;
1606 # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
1607 $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
1608 $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
1609 $x; # MBI bnorm is no-op, so dont call it
1612 ##############################################################################
1613 # internal calculation routines
1617 # return copy as a bigint representation of this BigFloat number
1618 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1620 my $z = $x->{_m}->copy();
1621 if ($x->{_e}->{sign} eq '-') # < 0
1623 $x->{_e}->{sign} = '+'; # flip
1624 $z->brsft($x->{_e},10);
1625 $x->{_e}->{sign} = '-'; # flip back
1627 elsif (!$x->{_e}->is_zero()) # > 0
1629 $z->blsft($x->{_e},10);
1631 $z->{sign} = $x->{sign};
1638 my $class = ref($x) || $x;
1639 $x = $class->new(shift) unless ref($x);
1641 return 1 if $x->{_m}->is_zero();
1642 my $len = $x->{_m}->length();
1643 $len += $x->{_e} if $x->{_e}->sign() eq '+';
1646 my $t = Math::BigInt::bzero();
1647 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1658 Math::BigFloat - Arbitrary size floating point math package
1665 $x = Math::BigFloat->new($str); # defaults to 0
1666 $nan = Math::BigFloat->bnan(); # create a NotANumber
1667 $zero = Math::BigFloat->bzero(); # create a +0
1668 $inf = Math::BigFloat->binf(); # create a +inf
1669 $inf = Math::BigFloat->binf('-'); # create a -inf
1670 $one = Math::BigFloat->bone(); # create a +1
1671 $one = Math::BigFloat->bone('-'); # create a -1
1674 $x->is_zero(); # true if arg is +0
1675 $x->is_nan(); # true if arg is NaN
1676 $x->is_one(); # true if arg is +1
1677 $x->is_one('-'); # true if arg is -1
1678 $x->is_odd(); # true if odd, false for even
1679 $x->is_even(); # true if even, false for odd
1680 $x->is_positive(); # true if >= 0
1681 $x->is_negative(); # true if < 0
1682 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
1684 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
1685 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
1686 $x->sign(); # return the sign, either +,- or NaN
1687 $x->digit($n); # return the nth digit, counting from right
1688 $x->digit(-$n); # return the nth digit, counting from left
1690 # The following all modify their first argument:
1693 $x->bzero(); # set $i to 0
1694 $x->bnan(); # set $i to NaN
1695 $x->bone(); # set $x to +1
1696 $x->bone('-'); # set $x to -1
1697 $x->binf(); # set $x to inf
1698 $x->binf('-'); # set $x to -inf
1700 $x->bneg(); # negation
1701 $x->babs(); # absolute value
1702 $x->bnorm(); # normalize (no-op)
1703 $x->bnot(); # two's complement (bit wise not)
1704 $x->binc(); # increment x by 1
1705 $x->bdec(); # decrement x by 1
1707 $x->badd($y); # addition (add $y to $x)
1708 $x->bsub($y); # subtraction (subtract $y from $x)
1709 $x->bmul($y); # multiplication (multiply $x by $y)
1710 $x->bdiv($y); # divide, set $i to quotient
1711 # return (quo,rem) or quo if scalar
1713 $x->bmod($y); # modulus
1714 $x->bpow($y); # power of arguments (a**b)
1715 $x->blsft($y); # left shift
1716 $x->brsft($y); # right shift
1717 # return (quo,rem) or quo if scalar
1719 $x->blog($base); # logarithm of $x, base defaults to e
1720 # (other bases than e not supported yet)
1722 $x->band($y); # bit-wise and
1723 $x->bior($y); # bit-wise inclusive or
1724 $x->bxor($y); # bit-wise exclusive or
1725 $x->bnot(); # bit-wise not (two's complement)
1727 $x->bsqrt(); # calculate square-root
1728 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
1730 $x->bround($N); # accuracy: preserver $N digits
1731 $x->bfround($N); # precision: round to the $Nth digit
1733 # The following do not modify their arguments:
1734 bgcd(@values); # greatest common divisor
1735 blcm(@values); # lowest common multiplicator
1737 $x->bstr(); # return string
1738 $x->bsstr(); # return string in scientific notation
1740 $x->bfloor(); # return integer less or equal than $x
1741 $x->bceil(); # return integer greater or equal than $x
1743 $x->exponent(); # return exponent as BigInt
1744 $x->mantissa(); # return mantissa as BigInt
1745 $x->parts(); # return (mantissa,exponent) as BigInt
1747 $x->length(); # number of digits (w/o sign and '.')
1748 ($l,$f) = $x->length(); # number of digits, and length of fraction
1752 All operators (inlcuding basic math operations) are overloaded if you
1753 declare your big floating point numbers as
1755 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1757 Operations with overloaded operators preserve the arguments, which is
1758 exactly what you expect.
1760 =head2 Canonical notation
1762 Input to these routines are either BigFloat objects, or strings of the
1763 following four forms:
1777 C</^[+-]\d+E[+-]?\d+$/>
1781 C</^[+-]\d*\.\d+E[+-]?\d+$/>
1785 all with optional leading and trailing zeros and/or spaces. Additonally,
1786 numbers are allowed to have an underscore between any two digits.
1788 Empty strings as well as other illegal numbers results in 'NaN'.
1790 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1791 are always stored in normalized form. On a string, it creates a BigFloat
1796 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1798 The string output will always have leading and trailing zeros stripped and drop
1799 a plus sign. C<bstr()> will give you always the form with a decimal point,
1800 while C<bsstr()> (for scientific) gives you the scientific notation.
1802 Input bstr() bsstr()
1804 ' -123 123 123' '-123123123' '-123123123E0'
1805 '00.0123' '0.0123' '123E-4'
1806 '123.45E-2' '1.2345' '12345E-4'
1807 '10E+3' '10000' '1E4'
1809 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1810 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1811 return either undef, <0, 0 or >0 and are suited for sort.
1813 Actual math is done by using BigInts to represent the mantissa and exponent.
1814 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1815 represent the result when input arguments are not numbers, as well as
1816 the result of dividing by zero.
1818 =head2 C<mantissa()>, C<exponent()> and C<parts()>
1820 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1821 as BigInts such that:
1823 $m = $x->mantissa();
1824 $e = $x->exponent();
1825 $y = $m * ( 10 ** $e );
1826 print "ok\n" if $x == $y;
1828 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
1830 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
1832 Currently the mantissa is reduced as much as possible, favouring higher
1833 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
1834 This might change in the future, so do not depend on it.
1836 =head2 Accuracy vs. Precision
1838 See also: L<Rounding|Rounding>.
1840 Math::BigFloat supports both precision and accuracy. For a full documentation,
1841 examples and tips on these topics please see the large section in
1844 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
1845 a operation consumes all resources, each operation produces no more than
1846 C<Math::BigFloat::precision()> digits.
1848 In case the result of one operation has more precision than specified,
1849 it is rounded. The rounding mode taken is either the default mode, or the one
1850 supplied to the operation after the I<scale>:
1852 $x = Math::BigFloat->new(2);
1853 Math::BigFloat::precision(5); # 5 digits max
1854 $y = $x->copy()->bdiv(3); # will give 0.66666
1855 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1856 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
1857 Math::BigFloat::round_mode('zero');
1858 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1864 =item ffround ( +$scale )
1866 Rounds to the $scale'th place left from the '.', counting from the dot.
1867 The first digit is numbered 1.
1869 =item ffround ( -$scale )
1871 Rounds to the $scale'th place right from the '.', counting from the dot.
1875 Rounds to an integer.
1877 =item fround ( +$scale )
1879 Preserves accuracy to $scale digits from the left (aka significant digits)
1880 and pads the rest with zeros. If the number is between 1 and -1, the
1881 significant digits count from the first non-zero after the '.'
1883 =item fround ( -$scale ) and fround ( 0 )
1885 These are effetively no-ops.
1889 All rounding functions take as a second parameter a rounding mode from one of
1890 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
1892 The default rounding mode is 'even'. By using
1893 C<< Math::BigFloat::round_mode($round_mode); >> you can get and set the default
1894 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
1895 no longer supported.
1896 The second parameter to the round functions then overrides the default
1899 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
1900 'trunc' as rounding mode to make it equivalent to:
1905 You can override this by passing the desired rounding mode as parameter to
1908 $x = Math::BigFloat->new(2.5);
1909 $y = $x->as_number('odd'); # $y = 3
1915 =head1 Autocreating constants
1917 After C<use Math::BigFloat ':constant'> all the floating point constants
1918 in the given scope are converted to C<Math::BigFloat>. This conversion
1919 happens at compile time.
1923 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
1925 prints the value of C<2E-100>. Note that without conversion of
1926 constants the expression 2E-100 will be calculated as normal floating point
1935 The following does not work yet:
1937 $m = $x->mantissa();
1938 $e = $x->exponent();
1939 $y = $m * ( 10 ** $e );
1940 print "ok\n" if $x == $y;
1944 There is no fmod() function yet.
1952 =item stringify, bstr()
1954 Both stringify and bstr() now drop the leading '+'. The old code would return
1955 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
1956 reasoning and details.
1960 The following will probably not do what you expect:
1962 print $c->bdiv(123.456),"\n";
1964 It prints both quotient and reminder since print works in list context. Also,
1965 bdiv() will modify $c, so be carefull. You probably want to use
1967 print $c / 123.456,"\n";
1968 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
1972 =item Modifying and =
1976 $x = Math::BigFloat->new(5);
1979 It will not do what you think, e.g. making a copy of $x. Instead it just makes
1980 a second reference to the B<same> object and stores it in $y. Thus anything
1981 that modifies $x will modify $y, and vice versa.
1984 print "$x, $y\n"; # prints '10, 10'
1986 If you want a true copy of $x, use:
1990 See also the documentation in L<overload> regarding C<=>.
1994 C<bpow()> now modifies the first argument, unlike the old code which left
1995 it alone and only returned the result. This is to be consistent with
1996 C<badd()> etc. The first will modify $x, the second one won't:
1998 print bpow($x,$i),"\n"; # modify $x
1999 print $x->bpow($i),"\n"; # ditto
2000 print $x ** $i,"\n"; # leave $x alone
2006 This program is free software; you may redistribute it and/or modify it under
2007 the same terms as Perl itself.
2011 Mark Biggar, overloaded interface by Ilya Zakharevich.
2012 Completely rewritten by Tels http://bloodgate.com in 2001.