1 package Math::BigFloat;
4 # Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
7 # The following hash values are internally used:
8 # _e: exponent (BigInt)
9 # _m: mantissa (absolute BigInt)
10 # sign: +,-,"NaN" if not a number
13 # _f: flags, used to signal MBI not to touch our private parts
20 @ISA = qw( Exporter Math::BigInt);
23 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode/;
24 use vars qw/$upgrade $downgrade/;
25 my $class = "Math::BigFloat";
28 '<=>' => sub { $_[2] ?
29 ref($_[0])->bcmp($_[1],$_[0]) :
30 ref($_[0])->bcmp($_[0],$_[1])},
31 'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
34 ##############################################################################
35 # global constants, flags and accessory
37 use constant MB_NEVER_ROUND => 0x0001;
41 # constant for easier life
44 # class constants, use Class->constant_name() to access
45 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
52 my $MBI = 'Math::BigInt'; # the package we are using for our private parts
53 # changable by use Math::BigFloat with => 'package'
55 ##############################################################################
56 # the old code had $rnd_mode, so we need to support it, too
58 sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
59 sub FETCH { return $round_mode; }
60 sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
65 tie $rnd_mode, 'Math::BigFloat';
68 ##############################################################################
70 # in case we call SUPER::->foo() and this wants to call modify()
71 # sub modify () { 0; }
74 # valid method aliases for AUTOLOAD
75 my %methods = map { $_ => 1 }
76 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
77 fint facmp fcmp fzero fnan finf finc fdec flog ffac
78 fceil ffloor frsft flsft fone flog
80 # valid method's that can be hand-ed up (for AUTOLOAD)
81 my %hand_ups = map { $_ => 1 }
82 qw / is_nan is_inf is_negative is_positive
83 accuracy precision div_scale round_mode fneg fabs babs fnot
84 objectify upgrade downgrade
88 sub method_alias { return exists $methods{$_[0]||''}; }
89 sub method_hand_up { return exists $hand_ups{$_[0]||''}; }
92 ##############################################################################
97 # create a new BigFloat object from a string or another bigfloat object.
100 # sign => sign (+/-), or "NaN"
102 my ($class,$wanted,@r) = @_;
104 # avoid numify-calls by not using || on $wanted!
105 return $class->bzero() if !defined $wanted; # default to 0
106 return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
108 my $self = {}; bless $self, $class;
109 # shortcut for bigints and its subclasses
110 if ((ref($wanted)) && (ref($wanted) ne $class))
112 $self->{_m} = $wanted->as_number(); # get us a bigint copy
113 $self->{_e} = $MBI->bzero();
115 $self->{sign} = $wanted->sign();
116 return $self->bnorm();
119 # handle '+inf', '-inf' first
120 if ($wanted =~ /^[+-]?inf$/)
122 return $downgrade->new($wanted) if $downgrade;
124 $self->{_e} = $MBI->bzero();
125 $self->{_m} = $MBI->bzero();
126 $self->{sign} = $wanted;
127 $self->{sign} = '+inf' if $self->{sign} eq 'inf';
128 return $self->bnorm();
130 #print "new string '$wanted'\n";
131 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
134 die "$wanted is not a number initialized to $class" if !$NaNOK;
136 return $downgrade->bnan() if $downgrade;
138 $self->{_e} = $MBI->bzero();
139 $self->{_m} = $MBI->bzero();
140 $self->{sign} = $nan;
144 # make integer from mantissa by adjusting exp, then convert to bigint
145 # undef,undef to signal MBI that we don't need no bloody rounding
146 $self->{_e} = $MBI->new("$$es$$ev",undef,undef); # exponent
147 $self->{_m} = $MBI->new("$$miv$$mfv",undef,undef); # create mant.
148 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
149 $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;
150 $self->{sign} = $$mis;
152 # if downgrade, inf, NaN or integers go down
154 if ($downgrade && $self->{_e}->{sign} eq '+')
156 # print "downgrading $$miv$$mfv"."E$$es$$ev";
157 if ($self->{_e}->is_zero())
159 $self->{_m}->{sign} = $$mis; # negative if wanted
160 return $downgrade->new($self->{_m});
162 return $downgrade->new("$$mis$$miv$$mfv"."E$$es$$ev");
164 # print "mbf new $self->{sign} $self->{_m} e $self->{_e} ",ref($self),"\n";
165 $self->bnorm()->round(@r); # first normalize, then round
170 # used by parent class bone() to initialize number to 1
172 $self->{_m} = $MBI->bzero();
173 $self->{_e} = $MBI->bzero();
178 # used by parent class bone() to initialize number to 1
180 $self->{_m} = $MBI->bzero();
181 $self->{_e} = $MBI->bzero();
186 # used by parent class bone() to initialize number to 1
188 $self->{_m} = $MBI->bone();
189 $self->{_e} = $MBI->bzero();
194 # used by parent class bone() to initialize number to 1
196 $self->{_m} = $MBI->bzero();
197 $self->{_e} = $MBI->bone();
202 my ($self,$class) = @_;
203 return if $class =~ /^Math::BigInt/; # we aren't one of these
204 UNIVERSAL::isa($self,$class);
209 # return (later set?) configuration data as hash ref
210 my $class = shift || 'Math::BigFloat';
212 my $cfg = $MBI->config();
215 $cfg->{class} = $class;
218 qw/upgrade downgrade precision accuracy round_mode VERSION div_scale/)
220 $cfg->{lc($_)} = ${"${class}::$_"};
225 ##############################################################################
226 # string conversation
230 # (ref to BFLOAT or num_str ) return num_str
231 # Convert number from internal format to (non-scientific) string format.
232 # internal format is always normalized (no leading zeros, "-0" => "+0")
233 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
234 #my $x = shift; my $class = ref($x) || $x;
235 #$x = $class->new(shift) unless ref($x);
237 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
238 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
239 if ($x->{sign} !~ /^[+-]$/)
241 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
245 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
247 my $not_zero = ! $x->is_zero();
250 $es = $x->{_m}->bstr();
251 $len = CORE::length($es);
252 if (!$x->{_e}->is_zero())
254 if ($x->{_e}->sign() eq '-')
257 if ($x->{_e} <= -$len)
259 # print "style: 0.xxxx\n";
260 my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
261 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
265 # print "insert '.' at $x->{_e} in '$es'\n";
266 substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
272 $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
276 $es = $x->{sign}.$es if $x->{sign} eq '-';
277 # if set accuracy or precision, pad with zeros
278 if ((defined $x->{_a}) && ($not_zero))
280 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
281 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
282 $zeros = $x->{_a} - $len if $cad != $len;
283 $es .= $dot.'0' x $zeros if $zeros > 0;
285 elsif ($x->{_p} || 0 < 0)
287 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
288 my $zeros = -$x->{_p} + $cad;
289 $es .= $dot.'0' x $zeros if $zeros > 0;
296 # (ref to BFLOAT or num_str ) return num_str
297 # Convert number from internal format to scientific string format.
298 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
299 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
300 #my $x = shift; my $class = ref($x) || $x;
301 #$x = $class->new(shift) unless ref($x);
303 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
304 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
305 if ($x->{sign} !~ /^[+-]$/)
307 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
310 my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
312 $x->{_m}->bstr().$sep.$x->{_e}->bstr();
317 # Make a number from a BigFloat object
318 # simple return string and let Perl's atoi()/atof() handle the rest
319 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
323 ##############################################################################
324 # public stuff (usually prefixed with "b")
327 # todo: this must be overwritten and return NaN for non-integer values
328 # band(), bior(), bxor(), too
331 # $class->SUPER::bnot($class,@_);
336 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
337 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
338 my ($self,$x,$y) = objectify(2,@_);
340 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
342 # handle +-inf and NaN
343 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
344 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
345 return +1 if $x->{sign} eq '+inf';
346 return -1 if $x->{sign} eq '-inf';
347 return -1 if $y->{sign} eq '+inf';
351 # check sign for speed first
352 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
353 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
356 my $xz = $x->is_zero();
357 my $yz = $y->is_zero();
358 return 0 if $xz && $yz; # 0 <=> 0
359 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
360 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
362 # adjust so that exponents are equal
363 my $lxm = $x->{_m}->length();
364 my $lym = $y->{_m}->length();
365 # the numify somewhat limits our length, but makes it much faster
366 my $lx = $lxm + $x->{_e}->numify();
367 my $ly = $lym + $y->{_e}->numify();
368 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
369 return $l <=> 0 if $l != 0;
371 # lengths (corrected by exponent) are equal
372 # so make mantissa equal length by padding with zero (shift left)
373 my $diff = $lxm - $lym;
374 my $xm = $x->{_m}; # not yet copy it
378 $ym = $y->{_m}->copy()->blsft($diff,10);
382 $xm = $x->{_m}->copy()->blsft(-$diff,10);
384 my $rc = $xm->bacmp($ym);
385 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
391 # Compares 2 values, ignoring their signs.
392 # Returns one of undef, <0, =0, >0. (suitable for sort)
393 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
394 my ($self,$x,$y) = objectify(2,@_);
396 # handle +-inf and NaN's
397 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
399 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
400 return 0 if ($x->is_inf() && $y->is_inf());
401 return 1 if ($x->is_inf() && !$y->is_inf());
406 my $xz = $x->is_zero();
407 my $yz = $y->is_zero();
408 return 0 if $xz && $yz; # 0 <=> 0
409 return -1 if $xz && !$yz; # 0 <=> +y
410 return 1 if $yz && !$xz; # +x <=> 0
412 # adjust so that exponents are equal
413 my $lxm = $x->{_m}->length();
414 my $lym = $y->{_m}->length();
415 # the numify somewhat limits our length, but makes it much faster
416 my $lx = $lxm + $x->{_e}->numify();
417 my $ly = $lym + $y->{_e}->numify();
419 return $l <=> 0 if $l != 0;
421 # lengths (corrected by exponent) are equal
422 # so make mantissa equal-length by padding with zero (shift left)
423 my $diff = $lxm - $lym;
424 my $xm = $x->{_m}; # not yet copy it
428 $ym = $y->{_m}->copy()->blsft($diff,10);
432 $xm = $x->{_m}->copy()->blsft(-$diff,10);
434 $xm->bacmp($ym) <=> 0;
439 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
440 # return result as BFLOAT
441 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
443 # inf and NaN handling
444 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
447 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
449 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
451 # +inf++inf or -inf+-inf => same, rest is NaN
452 return $x if $x->{sign} eq $y->{sign};
455 # +-inf + something => +inf; something +-inf => +-inf
456 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
460 return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
461 ((!$x->isa($self)) || (!$y->isa($self)));
463 # speed: no add for 0+y or x+0
464 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
465 if ($x->is_zero()) # 0+y
467 # make copy, clobbering up x (modify in place!)
468 $x->{_e} = $y->{_e}->copy();
469 $x->{_m} = $y->{_m}->copy();
470 $x->{sign} = $y->{sign} || $nan;
471 return $x->round($a,$p,$r,$y);
474 # take lower of the two e's and adapt m1 to it to match m2
476 $e = $MBI->bzero() if !defined $e; # if no BFLOAT ?
477 $e = $e->copy(); # make copy (didn't do it yet)
479 my $add = $y->{_m}->copy();
480 if ($e->{sign} eq '-') # < 0
482 my $e1 = $e->copy()->babs();
483 #$x->{_m} *= (10 ** $e1);
484 $x->{_m}->blsft($e1,10);
485 $x->{_e} += $e; # need the sign of e
487 elsif (!$e->is_zero()) # > 0
492 # else: both e are the same, so just leave them
493 $x->{_m}->{sign} = $x->{sign}; # fiddle with signs
494 $add->{sign} = $y->{sign};
495 $x->{_m} += $add; # finally do add/sub
496 $x->{sign} = $x->{_m}->{sign}; # re-adjust signs
497 $x->{_m}->{sign} = '+'; # mantissa always positiv
498 # delete trailing zeros, then round
499 return $x->bnorm()->round($a,$p,$r,$y);
504 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
505 # subtract second arg from first, modify first
506 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
508 if ($y->is_zero()) # still round for not adding zero
510 return $x->round($a,$p,$r);
513 $y->{sign} =~ tr/+\-/-+/; # does nothing for NaN
514 $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
515 $y->{sign} =~ tr/+\-/-+/; # refix $y (does nothing for NaN)
516 $x; # already rounded by badd()
521 # increment arg by one
522 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
524 if ($x->{_e}->sign() eq '-')
526 return $x->badd($self->bone(),$a,$p,$r); # digits after dot
529 if (!$x->{_e}->is_zero())
531 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
535 if ($x->{sign} eq '+')
538 return $x->bnorm()->bround($a,$p,$r);
540 elsif ($x->{sign} eq '-')
543 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
544 return $x->bnorm()->bround($a,$p,$r);
546 # inf, nan handling etc
547 $x->badd($self->__one(),$a,$p,$r); # does round
552 # decrement arg by one
553 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
555 if ($x->{_e}->sign() eq '-')
557 return $x->badd($self->bone('-'),$a,$p,$r); # digits after dot
560 if (!$x->{_e}->is_zero())
562 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
566 my $zero = $x->is_zero();
568 if (($x->{sign} eq '-') || $zero)
571 $x->{sign} = '-' if $zero; # 0 => 1 => -1
572 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
573 return $x->bnorm()->round($a,$p,$r);
576 elsif ($x->{sign} eq '+')
579 return $x->bnorm()->round($a,$p,$r);
581 # inf, nan handling etc
582 $x->badd($self->bone('-'),$a,$p,$r); # does round
587 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(2,@_);
589 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
593 # Taylor: | u 1 u^3 1 u^5 |
594 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
595 # |_ v 3 v^3 5 v^5 _|
597 # This takes much more steps to calculate the result:
600 # Taylor: | u 1 u^2 1 u^3 |
601 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
602 # |_ x 2 x^2 3 x^3 _|
604 # we need to limit the accuracy to protect against overflow
607 my @params = $x->_find_round_parameters($a,$p,$r);
609 # no rounding at all, so must use fallback
610 if (scalar @params == 1)
612 # simulate old behaviour
613 $params[1] = $self->div_scale(); # and round to it as accuracy
614 $scale = $params[1]+4; # at least four more for proper round
615 $params[3] = $r; # round mode by caller or undef
616 $fallback = 1; # to clear a/p afterwards
620 # the 4 below is empirical, and there might be cases where it is not
622 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
625 return $x->bzero(@params) if $x->is_one();
626 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
627 #return $x->bone('+',@params) if $x->bcmp($base) == 0;
629 # when user set globals, they would interfere with our calculation, so
630 # disable then and later re-enable them
632 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
633 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
634 # we also need to disable any set A or P on $x (_find_round_parameters took
635 # them already into account), since these would interfere, too
636 delete $x->{_a}; delete $x->{_p};
637 # need to disable $upgrade in BigInt, to avoid deep recursion
638 local $Math::BigInt::upgrade = undef;
640 my ($case,$limit,$v,$u,$below,$factor,$two,$next,$over,$f);
643 #if ($x <= Math::BigFloat->new("0.5"))
646 # print "case $case $x < 0.5\n";
647 $v = $x->copy(); $v->binc(); # v = x+1
648 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
649 $x->bdiv($v,$scale); # first term: u/v
652 $u *= $u; $v *= $v; # u^2, v^2
653 $below->bmul($v); # u^3, v^3
655 $factor = $self->new(3); $f = $self->new(2);
660 # print "case 1 $x > 0.5\n";
661 # $v = $x->copy(); # v = x
662 # $u = $x->copy(); $u->bdec(); # u = x-1;
663 # $x->bdec(); $x->bdiv($v,$scale); # first term: x-1/x
664 # $below = $v->copy();
665 # $over = $u->copy();
666 # $below->bmul($v); # u^2, v^2
668 # $factor = $self->new(2); $f = $self->bone();
670 $limit = $self->new("1E-". ($scale-1));
674 # we calculate the next term, and add it to the last
675 # when the next term is below our limit, it won't affect the outcome
676 # anymore, so we stop
677 $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
678 last if $next->bcmp($limit) <= 0;
681 # calculate things for the next term
682 $over *= $u; $below *= $v; $factor->badd($f);
685 $x->bmul(2) if $case == 0;
686 #print "took $steps steps\n";
688 # shortcut to not run trough _find_round_parameters again
689 if (defined $params[1])
691 $x->bround($params[1],$params[3]); # then round accordingly
695 $x->bfround($params[2],$params[3]); # then round accordingly
699 # clear a/p after round, since user did not request it
700 $x->{_a} = undef; $x->{_p} = undef;
703 $$abr = $ab; $$pbr = $pb;
710 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
711 # does not modify arguments, but returns new object
712 # Lowest Common Multiplicator
714 my ($self,@arg) = objectify(0,@_);
715 my $x = $self->new(shift @arg);
716 while (@arg) { $x = _lcm($x,shift @arg); }
722 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
723 # does not modify arguments, but returns new object
724 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
726 my ($self,@arg) = objectify(0,@_);
727 my $x = $self->new(shift @arg);
728 while (@arg) { $x = _gcd($x,shift @arg); }
732 ###############################################################################
733 # is_foo methods (is_negative, is_positive are inherited from BigInt)
737 # return true if arg (BFLOAT or num_str) is an integer
738 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
740 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
741 $x->{_e}->{sign} eq '+'; # 1e-1 => no integer
747 # return true if arg (BFLOAT or num_str) is zero
748 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
750 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
756 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
757 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
759 my $sign = shift || ''; $sign = '+' if $sign ne '-';
761 if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
767 # return true if arg (BFLOAT or num_str) is odd or false if even
768 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
770 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
771 ($x->{_e}->is_zero() && $x->{_m}->is_odd());
777 # return true if arg (BINT or num_str) is even or false if odd
778 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
780 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
781 return 1 if ($x->{_e}->{sign} eq '+' # 123.45 is never
782 && $x->{_m}->is_even()); # but 1200 is
788 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
789 # (BINT or num_str, BINT or num_str) return BINT
790 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
792 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
795 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
797 return $x->bnan() if $x->is_zero() || $y->is_zero();
798 # result will always be +-inf:
799 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
800 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
801 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
802 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
803 return $x->binf('-');
806 return $x->bzero() if $x->is_zero() || $y->is_zero();
808 return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
809 ((!$x->isa($self)) || (!$y->isa($self)));
811 # aEb * cEd = (a*c)E(b+d)
812 $x->{_m}->bmul($y->{_m});
813 $x->{_e}->badd($y->{_e});
815 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
816 return $x->bnorm()->round($a,$p,$r,$y);
821 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
822 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
823 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
825 return $self->_div_inf($x,$y)
826 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
828 # x== 0 # also: or y == 1 or y == -1
829 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
832 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
834 # we need to limit the accuracy to protect against overflow
837 my @params = $x->_find_round_parameters($a,$p,$r,$y);
839 # no rounding at all, so must use fallback
840 if (scalar @params == 1)
842 # simulate old behaviour
843 $params[1] = $self->div_scale(); # and round to it as accuracy
844 $scale = $params[1]+4; # at least four more for proper round
845 $params[3] = $r; # round mode by caller or undef
846 $fallback = 1; # to clear a/p afterwards
850 # the 4 below is empirical, and there might be cases where it is not
852 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
854 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
855 $scale = $lx if $lx > $scale;
856 $scale = $ly if $ly > $scale;
857 my $diff = $ly - $lx;
858 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
860 # make copy of $x in case of list context for later reminder calculation
862 if (wantarray && !$y->is_one())
867 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
869 # check for / +-1 ( +/- 1E0)
872 # promote BigInts and it's subclasses (except when already a BigFloat)
873 $y = $self->new($y) unless $y->isa('Math::BigFloat');
875 #print "bdiv $y ",ref($y),"\n";
876 # need to disable $upgrade in BigInt, to avoid deep recursion
877 local $Math::BigInt::upgrade = undef; # should be parent class vs MBI
879 # calculate the result to $scale digits and then round it
880 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
881 $x->{_m}->blsft($scale,10);
882 $x->{_m}->bdiv( $y->{_m} ); # a/c
883 $x->{_e}->bsub( $y->{_e} ); # b-d
884 $x->{_e}->bsub($scale); # correct for 10**scale
885 $x->bnorm(); # remove trailing 0's
888 # shortcut to not run trough _find_round_parameters again
889 if (defined $params[1])
891 $x->bround($params[1],$params[3]); # then round accordingly
895 $x->bfround($params[2],$params[3]); # then round accordingly
899 # clear a/p after round, since user did not request it
900 $x->{_a} = undef; $x->{_p} = undef;
907 $rem->bmod($y,$params[1],$params[2],$params[3]); # copy already done
911 $rem = $self->bzero();
915 # clear a/p after round, since user did not request it
916 $rem->{_a} = undef; $rem->{_p} = undef;
925 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
926 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
928 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
930 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
931 return $re->round($a,$p,$r,$y);
933 return $x->bnan() if $x->is_zero() && $y->is_zero();
934 return $x if $y->is_zero();
935 return $x->bnan() if $x->is_nan() || $y->is_nan();
936 return $x->bzero() if $y->is_one() || $x->is_zero();
938 # inf handling is missing here
940 my $cmp = $x->bacmp($y); # equal or $x < $y?
941 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
943 # only $y of the operands negative?
944 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
946 $x->{sign} = $y->{sign}; # calc sign first
947 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
949 my $ym = $y->{_m}->copy();
952 $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
954 # if $y has digits after dot
955 my $shifty = 0; # correct _e of $x by this
956 if ($y->{_e}->{sign} eq '-') # has digits after dot
958 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
959 $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
960 $x->blsft($shifty,10); # 123 => 1230, $y->{_m} is already 25
962 # $ym is now mantissa of $y based on exponent 0
964 my $shiftx = 0; # correct _e of $x by this
965 if ($x->{_e}->{sign} eq '-') # has digits after dot
967 # 123.4 % 20 => 1234 % 200
968 $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
969 $ym->blsft($shiftx,10);
971 # 123e1 % 20 => 1230 % 20
972 if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
974 $x->{_m}->blsft($x->{_e},10);
976 $x->{_e} = $MBI->bzero() unless $x->{_e}->is_zero();
978 $x->{_e}->bsub($shiftx) if $shiftx != 0;
979 $x->{_e}->bsub($shifty) if $shifty != 0;
981 # now mantissas are equalized, exponent of $x is adjusted, so calc result
985 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
988 if ($neg != 0) # one of them negative => correct in place
993 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
997 $x->round($a,$p,$r,$y); # round and return
1002 # calculate square root; this should probably
1003 # use a different test to see whether the accuracy we want is...
1004 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1006 return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
1007 return $x if $x->{sign} eq '+inf'; # +inf
1008 return $x if $x->is_zero() || $x->is_one();
1010 # we need to limit the accuracy to protect against overflow
1013 my @params = $x->_find_round_parameters($a,$p,$r);
1015 # no rounding at all, so must use fallback
1016 if (scalar @params == 1)
1018 # simulate old behaviour
1019 $params[1] = $self->div_scale(); # and round to it as accuracy
1020 $scale = $params[1]+4; # at least four more for proper round
1021 $params[3] = $r; # round mode by caller or undef
1022 $fallback = 1; # to clear a/p afterwards
1026 # the 4 below is empirical, and there might be cases where it is not
1028 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1031 # when user set globals, they would interfere with our calculation, so
1032 # disable them and later re-enable them
1034 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1035 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1036 # we also need to disable any set A or P on $x (_find_round_parameters took
1037 # them already into account), since these would interfere, too
1038 delete $x->{_a}; delete $x->{_p};
1039 # need to disable $upgrade in BigInt, to avoid deep recursion
1040 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1042 my $xas = $x->as_number();
1043 my $gs = $xas->copy()->bsqrt(); # some guess
1045 # print "guess $gs\n";
1046 if (($x->{_e}->{sign} ne '-') # guess can't be accurate if there are
1047 # digits after the dot
1048 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
1051 $x->{_m} = $gs; $x->{_e} = $MBI->bzero(); $x->bnorm();
1052 # shortcut to not run trough _find_round_parameters again
1053 if (defined $params[1])
1055 $x->bround($params[1],$params[3]); # then round accordingly
1059 $x->bfround($params[2],$params[3]); # then round accordingly
1063 # clear a/p after round, since user did not request it
1064 $x->{_a} = undef; $x->{_p} = undef;
1066 # re-enable A and P, upgrade is taken care of by "local"
1067 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1070 $gs = $self->new( $gs ); # BigInt to BigFloat
1072 my $lx = $x->{_m}->length();
1073 $scale = $lx if $scale < $lx;
1074 my $e = $self->new("1E-$scale"); # make test variable
1077 my $two = $self->new(2);
1079 # promote BigInts and it's subclasses (except when already a BigFloat)
1080 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1083 while ($diff->bacmp($e) >= 0)
1085 $rem = $y->copy()->bdiv($gs,$scale);
1086 $rem = $y->copy()->bdiv($gs,$scale)->badd($gs)->bdiv($two,$scale);
1087 $diff = $rem->copy()->bsub($gs);
1090 # copy over to modify $x
1091 $x->{_m} = $rem->{_m}; $x->{_e} = $rem->{_e};
1093 # shortcut to not run trough _find_round_parameters again
1094 if (defined $params[1])
1096 $x->bround($params[1],$params[3]); # then round accordingly
1100 $x->bfround($params[2],$params[3]); # then round accordingly
1104 # clear a/p after round, since user did not request it
1105 $x->{_a} = undef; $x->{_p} = undef;
1108 $$abr = $ab; $$pbr = $pb;
1114 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1115 # compute factorial numbers
1116 # modifies first argument
1117 my ($self,$x,@r) = objectify(1,@_);
1120 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1121 ($x->{_e}->{sign} ne '+')); # digits after dot?
1123 return $x->bone(@r) if $x->is_zero() || $x->is_one(); # 0 or 1 => 1
1125 # use BigInt's bfac() for faster calc
1126 $x->{_m}->blsft($x->{_e},10); # un-norm m
1127 $x->{_e}->bzero(); # norm $x again
1128 $x->{_m}->bfac(); # factorial
1129 $x->bnorm()->round(@r);
1134 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1135 my ($x,$y,$a,$p,$r) = @_;
1138 # we need to limit the accuracy to protect against overflow
1141 my @params = $x->_find_round_parameters($a,$p,$r);
1143 # no rounding at all, so must use fallback
1144 if (scalar @params == 1)
1146 # simulate old behaviour
1147 $params[1] = $self->div_scale(); # and round to it as accuracy
1148 $scale = $params[1]+4; # at least four more for proper round
1149 $params[3] = $r; # round mode by caller or undef
1150 $fallback = 1; # to clear a/p afterwards
1154 # the 4 below is empirical, and there might be cases where it is not
1156 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1159 # when user set globals, they would interfere with our calculation, so
1160 # disable then and later re-enable them
1162 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1163 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1164 # we also need to disable any set A or P on $x (_find_round_parameters took
1165 # them already into account), since these would interfere, too
1166 delete $x->{_a}; delete $x->{_p};
1167 # need to disable $upgrade in BigInt, to avoid deep recursion
1168 local $Math::BigInt::upgrade = undef;
1170 # split the second argument into its integer and fraction part
1171 # we calculate the result then from these two parts, like in
1172 # 2 ** 2.4 == (2 ** 2) * (2 ** 0.4)
1173 my $c = $self->new($y->as_number()); # integer part
1174 my $d = $y-$c; # fractional part
1175 my $xc = $x->copy(); # a temp. copy
1177 # now calculate binary fraction from the decimal fraction on the fly
1179 # 0.654 * 2 = 1.308 > 1 => 0.1 ( 1.308 - 1 = 0.308)
1180 # 0.308 * 2 = 0.616 < 1 => 0.10
1181 # 0.616 * 2 = 1.232 > 1 => 0.101 ( 1.232 - 1 = 0.232)
1183 # The process stops when the result is exactly one, or when we have
1186 # From the binary fraction we calculate the result as follows:
1187 # we assume the fraction ends in 1, and we remove this one first.
1188 # For each digit after the dot, assume 1 eq R and 0 eq XR, where R means
1189 # take square root and X multiply with the original X.
1195 last if $d->is_one(); # == 1
1199 $x->bsqrt(); $x->bmul($xc); $d->bdec(); # 1
1202 # assume fraction ends in 1
1206 $x->bmul( $xc->bpow($c) );
1208 elsif (!$c->is_zero())
1214 # shortcut to not run trough _find_round_parameters again
1215 if (defined $params[1])
1217 $x->bround($params[1],$params[3]); # then round accordingly
1221 $x->bfround($params[2],$params[3]); # then round accordingly
1225 # clear a/p after round, since user did not request it
1226 $x->{_a} = undef; $x->{_p} = undef;
1229 $$abr = $ab; $$pbr = $pb;
1235 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1236 my ($x,$y,$a,$p,$r) = @_;
1239 # if $y == 0.5, it is sqrt($x)
1240 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1244 # Taylor: | u u^2 u^3 |
1245 # x ** y = 1 + | --- + --- + * ----- + ... |
1248 # we need to limit the accuracy to protect against overflow
1251 my @params = $x->_find_round_parameters($a,$p,$r);
1253 # no rounding at all, so must use fallback
1254 if (scalar @params == 1)
1256 # simulate old behaviour
1257 $params[1] = $self->div_scale(); # and round to it as accuracy
1258 $scale = $params[1]+4; # at least four more for proper round
1259 $params[3] = $r; # round mode by caller or undef
1260 $fallback = 1; # to clear a/p afterwards
1264 # the 4 below is empirical, and there might be cases where it is not
1266 $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1269 # when user set globals, they would interfere with our calculation, so
1270 # disable then and later re-enable them
1272 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1273 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1274 # we also need to disable any set A or P on $x (_find_round_parameters took
1275 # them already into account), since these would interfere, too
1276 delete $x->{_a}; delete $x->{_p};
1277 # need to disable $upgrade in BigInt, to avoid deep recursion
1278 local $Math::BigInt::upgrade = undef;
1280 my ($limit,$v,$u,$below,$factor,$next,$over);
1282 $u = $x->copy()->blog($scale)->bmul($y);
1283 $v = $self->bone(); # 1
1284 $factor = $self->new(2); # 2
1285 $x->bone(); # first term: 1
1287 $below = $v->copy();
1290 $limit = $self->new("1E-". ($scale-1));
1294 # we calculate the next term, and add it to the last
1295 # when the next term is below our limit, it won't affect the outcome
1296 # anymore, so we stop
1297 $next = $over->copy()->bdiv($below,$scale);
1298 last if $next->bcmp($limit) <= 0;
1301 # calculate things for the next term
1302 $over *= $u; $below *= $factor; $factor->binc();
1306 # shortcut to not run trough _find_round_parameters again
1307 if (defined $params[1])
1309 $x->bround($params[1],$params[3]); # then round accordingly
1313 $x->bfround($params[2],$params[3]); # then round accordingly
1317 # clear a/p after round, since user did not request it
1318 $x->{_a} = undef; $x->{_p} = undef;
1321 $$abr = $ab; $$pbr = $pb;
1327 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1328 # compute power of two numbers, second arg is used as integer
1329 # modifies first argument
1331 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1333 return $x if $x->{sign} =~ /^[+-]inf$/;
1334 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1335 return $x->bone() if $y->is_zero();
1336 return $x if $x->is_one() || $y->is_one();
1338 return $x->_pow($y,$a,$p,$r) if !$y->is_int(); # non-integer power
1340 my $y1 = $y->as_number(); # make bigint
1342 if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1344 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
1345 return $y1->is_odd() ? $x : $x->babs(1);
1349 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1350 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1354 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1356 $x->{_m}->bpow($y1);
1357 $x->{_e}->bmul($y1);
1358 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1360 if ($y->{sign} eq '-')
1362 # modify $x in place!
1363 my $z = $x->copy(); $x->bzero()->binc();
1364 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
1366 $x->round($a,$p,$r,$y);
1369 ###############################################################################
1370 # rounding functions
1374 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1375 # $n == 0 means round to integer
1376 # expects and returns normalized numbers!
1377 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1379 return $x if $x->modify('bfround');
1381 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1382 return $x if !defined $scale; # no-op
1384 # never round a 0, +-inf, NaN
1387 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1390 return $x if $x->{sign} !~ /^[+-]$/;
1391 # print "MBF bfround $x to scale $scale mode $mode\n";
1393 # don't round if x already has lower precision
1394 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1396 $x->{_p} = $scale; # remember round in any case
1397 $x->{_a} = undef; # and clear A
1400 # print "bfround scale $scale e $x->{_e}\n";
1401 # round right from the '.'
1402 return $x if $x->{_e} >= 0; # nothing to round
1403 $scale = -$scale; # positive for simplicity
1404 my $len = $x->{_m}->length(); # length of mantissa
1405 my $dad = -$x->{_e}; # digits after dot
1406 my $zad = 0; # zeros after dot
1407 $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
1408 #print "scale $scale dad $dad zad $zad len $len\n";
1410 # number bsstr len zad dad
1411 # 0.123 123e-3 3 0 3
1412 # 0.0123 123e-4 3 1 4
1415 # 1.2345 12345e-4 5 0 4
1417 # do not round after/right of the $dad
1418 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1420 # round to zero if rounding inside the $zad, but not for last zero like:
1421 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1422 return $x->bzero() if $scale < $zad;
1423 if ($scale == $zad) # for 0.006, scale -3 and trunc
1429 # adjust round-point to be inside mantissa
1432 $scale = $scale-$zad;
1436 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
1437 $scale = $dbd+$scale;
1440 # print "round to $x->{_m} to $scale\n";
1444 # 123 => 100 means length(123) = 3 - $scale (2) => 1
1446 my $dbt = $x->{_m}->length();
1448 my $dbd = $dbt + $x->{_e};
1449 # should be the same, so treat it as this
1450 $scale = 1 if $scale == 0;
1451 # shortcut if already integer
1452 return $x if $scale == 1 && $dbt <= $dbd;
1453 # maximum digits before dot
1458 # not enough digits before dot, so round to zero
1461 elsif ( $scale == $dbd )
1468 $scale = $dbd - $scale;
1472 # print "using $scale for $x->{_m} with '$mode'\n";
1473 # pass sign to bround for rounding modes '+inf' and '-inf'
1474 $x->{_m}->{sign} = $x->{sign};
1475 $x->{_m}->bround($scale,$mode);
1476 $x->{_m}->{sign} = '+'; # fix sign back
1482 # accuracy: preserve $N digits, and overwrite the rest with 0's
1483 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1485 die ('bround() needs positive accuracy') if ($_[0] || 0) < 0;
1487 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1488 return $x if !defined $scale; # no-op
1490 return $x if $x->modify('bround');
1492 # scale is now either $x->{_a}, $accuracy, or the user parameter
1493 # test whether $x already has lower accuracy, do nothing in this case
1494 # but do round if the accuracy is the same, since a math operation might
1495 # want to round a number with A=5 to 5 digits afterwards again
1496 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1498 # scale < 0 makes no sense
1499 # never round a +-inf, NaN
1500 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
1502 # 1: $scale == 0 => keep all digits
1503 # 2: never round a 0
1504 # 3: if we should keep more digits than the mantissa has, do nothing
1505 if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1507 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1511 # pass sign to bround for '+inf' and '-inf' rounding modes
1512 $x->{_m}->{sign} = $x->{sign};
1513 $x->{_m}->bround($scale,$mode); # round mantissa
1514 $x->{_m}->{sign} = '+'; # fix sign back
1515 # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1516 $x->{_a} = $scale; # remember rounding
1517 $x->{_p} = undef; # and clear P
1518 $x->bnorm(); # del trailing zeros gen. by bround()
1523 # return integer less or equal then $x
1524 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1526 return $x if $x->modify('bfloor');
1528 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1530 # if $x has digits after dot
1531 if ($x->{_e}->{sign} eq '-')
1533 #$x->{_m}->brsft(-$x->{_e},10);
1535 #$x-- if $x->{sign} eq '-';
1537 $x->{_e}->{sign} = '+'; # negate e
1538 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1539 $x->{_e}->bzero(); # trunc/norm
1540 $x->{_m}->binc() if $x->{sign} eq '-'; # decrement if negative
1542 $x->round($a,$p,$r);
1547 # return integer greater or equal then $x
1548 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1550 return $x if $x->modify('bceil');
1551 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1553 # if $x has digits after dot
1554 if ($x->{_e}->{sign} eq '-')
1556 #$x->{_m}->brsft(-$x->{_e},10);
1558 #$x++ if $x->{sign} eq '+';
1560 $x->{_e}->{sign} = '+'; # negate e
1561 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1562 $x->{_e}->bzero(); # trunc/norm
1563 $x->{_m}->binc() if $x->{sign} eq '+'; # decrement if negative
1565 $x->round($a,$p,$r);
1570 # shift right by $y (divide by power of 2)
1571 my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1573 return $x if $x->modify('brsft');
1574 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1576 $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1577 $x->bdiv($n ** $y,$a,$p,$r,$y);
1582 # shift right by $y (divide by power of 2)
1583 my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1585 return $x if $x->modify('brsft');
1586 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1588 $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1589 $x->bmul($n ** $y,$a,$p,$r,$y);
1592 ###############################################################################
1596 # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1601 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
1602 # or falling back to MBI::bxxx()
1603 my $name = $AUTOLOAD;
1605 $name =~ s/.*:://; # split package
1607 if (!method_alias($name))
1611 # delayed load of Carp and avoid recursion
1613 Carp::croak ("Can't call a method without name");
1615 if (!method_hand_up($name))
1617 # delayed load of Carp and avoid recursion
1619 Carp::croak ("Can't call $class\-\>$name, not a valid method");
1621 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1623 return &{"$MBI"."::$name"}(@_);
1625 my $bname = $name; $bname =~ s/^f/b/;
1626 *{$class."::$name"} = \&$bname;
1632 # return a copy of the exponent
1633 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1635 if ($x->{sign} !~ /^[+-]$/)
1637 my $s = $x->{sign}; $s =~ s/^[+-]//;
1638 return $self->new($s); # -inf, +inf => +inf
1640 return $x->{_e}->copy();
1645 # return a copy of the mantissa
1646 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1648 if ($x->{sign} !~ /^[+-]$/)
1650 my $s = $x->{sign}; $s =~ s/^[+]//;
1651 return $self->new($s); # -inf, +inf => +inf
1653 my $m = $x->{_m}->copy(); # faster than going via bstr()
1654 $m->bneg() if $x->{sign} eq '-';
1661 # return a copy of both the exponent and the mantissa
1662 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1664 if ($x->{sign} !~ /^[+-]$/)
1666 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
1667 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
1669 my $m = $x->{_m}->copy(); # faster than going via bstr()
1670 $m->bneg() if $x->{sign} eq '-';
1671 return ($m,$x->{_e}->copy());
1674 ##############################################################################
1675 # private stuff (internal use only)
1681 my $lib = ''; my @a;
1682 for ( my $i = 0; $i < $l ; $i++)
1684 # print "at $_[$i] (",$_[$i+1]||'undef',")\n";
1685 if ( $_[$i] eq ':constant' )
1687 # this rest causes overlord er load to step in
1688 # print "overload @_\n";
1689 overload::constant float => sub { $self->new(shift); };
1691 elsif ($_[$i] eq 'upgrade')
1693 # this causes upgrading
1694 $upgrade = $_[$i+1]; # or undef to disable
1697 elsif ($_[$i] eq 'downgrade')
1699 # this causes downgrading
1700 $downgrade = $_[$i+1]; # or undef to disable
1703 elsif ($_[$i] eq 'lib')
1705 $lib = $_[$i+1] || ''; # default Calc
1708 elsif ($_[$i] eq 'with')
1710 $MBI = $_[$i+1] || 'Math::BigInt'; # default Math::BigInt
1720 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
1721 my $mbilib = eval { Math::BigInt->config()->{lib} };
1722 if ((defined $mbilib) && ($MBI eq 'Math::BigInt'))
1724 # MBI already loaded
1725 $MBI->import('lib',"$lib,$mbilib", 'objectify');
1729 # MBI not loaded, or with ne "Math::BigInt"
1730 $lib .= ",$mbilib" if defined $mbilib;
1732 # my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
1733 # my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
1734 # $file = File::Spec->catfile (@parts, $file);
1738 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
1739 # used in the same script, or eval inside import().
1740 my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
1741 my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
1742 $file = File::Spec->catfile (@parts, $file);
1743 eval { require $file; $MBI->import( lib => '$lib', 'objectify' ); }
1747 my $rc = "use $MBI lib => '$lib', 'objectify';";
1751 die ("Couldn't load $MBI: $! $@") if $@;
1753 # any non :constant stuff is handled by our parent, Exporter
1754 # even if @_ is empty, to give it a chance
1755 $self->SUPER::import(@a); # for subclasses
1756 $self->export_to_level(1,$self,@a); # need this, too
1761 # adjust m and e so that m is smallest possible
1762 # round number according to accuracy and precision settings
1763 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1765 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
1767 # if (!$x->{_m}->is_odd())
1769 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
1772 $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros);
1774 # for something like 0Ey, set y to 1, and -0 => +0
1775 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
1777 # this is to prevent automatically rounding when MBI's globals are set
1778 $x->{_m}->{_f} = MB_NEVER_ROUND;
1779 $x->{_e}->{_f} = MB_NEVER_ROUND;
1780 # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
1781 $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
1782 $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
1783 $x; # MBI bnorm is no-op, so dont call it
1786 ##############################################################################
1787 # internal calculation routines
1791 # return copy as a bigint representation of this BigFloat number
1792 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1794 my $z = $x->{_m}->copy();
1795 if ($x->{_e}->{sign} eq '-') # < 0
1797 $x->{_e}->{sign} = '+'; # flip
1798 $z->brsft($x->{_e},10);
1799 $x->{_e}->{sign} = '-'; # flip back
1801 elsif (!$x->{_e}->is_zero()) # > 0
1803 $z->blsft($x->{_e},10);
1805 $z->{sign} = $x->{sign};
1812 my $class = ref($x) || $x;
1813 $x = $class->new(shift) unless ref($x);
1815 return 1 if $x->{_m}->is_zero();
1816 my $len = $x->{_m}->length();
1817 $len += $x->{_e} if $x->{_e}->sign() eq '+';
1820 my $t = $MBI->bzero();
1821 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1832 Math::BigFloat - Arbitrary size floating point math package
1839 $x = Math::BigFloat->new($str); # defaults to 0
1840 $nan = Math::BigFloat->bnan(); # create a NotANumber
1841 $zero = Math::BigFloat->bzero(); # create a +0
1842 $inf = Math::BigFloat->binf(); # create a +inf
1843 $inf = Math::BigFloat->binf('-'); # create a -inf
1844 $one = Math::BigFloat->bone(); # create a +1
1845 $one = Math::BigFloat->bone('-'); # create a -1
1848 $x->is_zero(); # true if arg is +0
1849 $x->is_nan(); # true if arg is NaN
1850 $x->is_one(); # true if arg is +1
1851 $x->is_one('-'); # true if arg is -1
1852 $x->is_odd(); # true if odd, false for even
1853 $x->is_even(); # true if even, false for odd
1854 $x->is_positive(); # true if >= 0
1855 $x->is_negative(); # true if < 0
1856 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
1858 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
1859 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
1860 $x->sign(); # return the sign, either +,- or NaN
1861 $x->digit($n); # return the nth digit, counting from right
1862 $x->digit(-$n); # return the nth digit, counting from left
1864 # The following all modify their first argument:
1867 $x->bzero(); # set $i to 0
1868 $x->bnan(); # set $i to NaN
1869 $x->bone(); # set $x to +1
1870 $x->bone('-'); # set $x to -1
1871 $x->binf(); # set $x to inf
1872 $x->binf('-'); # set $x to -inf
1874 $x->bneg(); # negation
1875 $x->babs(); # absolute value
1876 $x->bnorm(); # normalize (no-op)
1877 $x->bnot(); # two's complement (bit wise not)
1878 $x->binc(); # increment x by 1
1879 $x->bdec(); # decrement x by 1
1881 $x->badd($y); # addition (add $y to $x)
1882 $x->bsub($y); # subtraction (subtract $y from $x)
1883 $x->bmul($y); # multiplication (multiply $x by $y)
1884 $x->bdiv($y); # divide, set $i to quotient
1885 # return (quo,rem) or quo if scalar
1887 $x->bmod($y); # modulus
1888 $x->bpow($y); # power of arguments (a**b)
1889 $x->blsft($y); # left shift
1890 $x->brsft($y); # right shift
1891 # return (quo,rem) or quo if scalar
1893 $x->blog($base); # logarithm of $x, base defaults to e
1894 # (other bases than e not supported yet)
1896 $x->band($y); # bit-wise and
1897 $x->bior($y); # bit-wise inclusive or
1898 $x->bxor($y); # bit-wise exclusive or
1899 $x->bnot(); # bit-wise not (two's complement)
1901 $x->bsqrt(); # calculate square-root
1902 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
1904 $x->bround($N); # accuracy: preserver $N digits
1905 $x->bfround($N); # precision: round to the $Nth digit
1907 # The following do not modify their arguments:
1908 bgcd(@values); # greatest common divisor
1909 blcm(@values); # lowest common multiplicator
1911 $x->bstr(); # return string
1912 $x->bsstr(); # return string in scientific notation
1914 $x->bfloor(); # return integer less or equal than $x
1915 $x->bceil(); # return integer greater or equal than $x
1917 $x->exponent(); # return exponent as BigInt
1918 $x->mantissa(); # return mantissa as BigInt
1919 $x->parts(); # return (mantissa,exponent) as BigInt
1921 $x->length(); # number of digits (w/o sign and '.')
1922 ($l,$f) = $x->length(); # number of digits, and length of fraction
1926 All operators (inlcuding basic math operations) are overloaded if you
1927 declare your big floating point numbers as
1929 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1931 Operations with overloaded operators preserve the arguments, which is
1932 exactly what you expect.
1934 =head2 Canonical notation
1936 Input to these routines are either BigFloat objects, or strings of the
1937 following four forms:
1951 C</^[+-]\d+E[+-]?\d+$/>
1955 C</^[+-]\d*\.\d+E[+-]?\d+$/>
1959 all with optional leading and trailing zeros and/or spaces. Additonally,
1960 numbers are allowed to have an underscore between any two digits.
1962 Empty strings as well as other illegal numbers results in 'NaN'.
1964 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1965 are always stored in normalized form. On a string, it creates a BigFloat
1970 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1972 The string output will always have leading and trailing zeros stripped and drop
1973 a plus sign. C<bstr()> will give you always the form with a decimal point,
1974 while C<bsstr()> (for scientific) gives you the scientific notation.
1976 Input bstr() bsstr()
1978 ' -123 123 123' '-123123123' '-123123123E0'
1979 '00.0123' '0.0123' '123E-4'
1980 '123.45E-2' '1.2345' '12345E-4'
1981 '10E+3' '10000' '1E4'
1983 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1984 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1985 return either undef, <0, 0 or >0 and are suited for sort.
1987 Actual math is done by using BigInts to represent the mantissa and exponent.
1988 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1989 represent the result when input arguments are not numbers, as well as
1990 the result of dividing by zero.
1992 =head2 C<mantissa()>, C<exponent()> and C<parts()>
1994 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1995 as BigInts such that:
1997 $m = $x->mantissa();
1998 $e = $x->exponent();
1999 $y = $m * ( 10 ** $e );
2000 print "ok\n" if $x == $y;
2002 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2004 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2006 Currently the mantissa is reduced as much as possible, favouring higher
2007 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2008 This might change in the future, so do not depend on it.
2010 =head2 Accuracy vs. Precision
2012 See also: L<Rounding|Rounding>.
2014 Math::BigFloat supports both precision and accuracy. For a full documentation,
2015 examples and tips on these topics please see the large section in
2018 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2019 a operation consumes all resources, each operation produces no more than
2020 C<Math::BigFloat::precision()> digits.
2022 In case the result of one operation has more precision than specified,
2023 it is rounded. The rounding mode taken is either the default mode, or the one
2024 supplied to the operation after the I<scale>:
2026 $x = Math::BigFloat->new(2);
2027 Math::BigFloat::precision(5); # 5 digits max
2028 $y = $x->copy()->bdiv(3); # will give 0.66666
2029 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2030 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
2031 Math::BigFloat::round_mode('zero');
2032 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2038 =item ffround ( +$scale )
2040 Rounds to the $scale'th place left from the '.', counting from the dot.
2041 The first digit is numbered 1.
2043 =item ffround ( -$scale )
2045 Rounds to the $scale'th place right from the '.', counting from the dot.
2049 Rounds to an integer.
2051 =item fround ( +$scale )
2053 Preserves accuracy to $scale digits from the left (aka significant digits)
2054 and pads the rest with zeros. If the number is between 1 and -1, the
2055 significant digits count from the first non-zero after the '.'
2057 =item fround ( -$scale ) and fround ( 0 )
2059 These are effetively no-ops.
2063 All rounding functions take as a second parameter a rounding mode from one of
2064 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2066 The default rounding mode is 'even'. By using
2067 C<< Math::BigFloat::round_mode($round_mode); >> you can get and set the default
2068 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2069 no longer supported.
2070 The second parameter to the round functions then overrides the default
2073 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
2074 'trunc' as rounding mode to make it equivalent to:
2079 You can override this by passing the desired rounding mode as parameter to
2082 $x = Math::BigFloat->new(2.5);
2083 $y = $x->as_number('odd'); # $y = 3
2089 =head1 Autocreating constants
2091 After C<use Math::BigFloat ':constant'> all the floating point constants
2092 in the given scope are converted to C<Math::BigFloat>. This conversion
2093 happens at compile time.
2097 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2099 prints the value of C<2E-100>. Note that without conversion of
2100 constants the expression 2E-100 will be calculated as normal floating point
2103 Please note that ':constant' does not affect integer constants, nor binary
2104 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2109 Math with the numbers is done (by default) by a module called
2110 Math::BigInt::Calc. This is equivalent to saying:
2112 use Math::BigFloat lib => 'Calc';
2114 You can change this by using:
2116 use Math::BigFloat lib => 'BitVect';
2118 The following would first try to find Math::BigInt::Foo, then
2119 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2121 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2123 Calc.pm uses as internal format an array of elements of some decimal base
2124 (usually 1e7, but this might be differen for some systems) with the least
2125 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2126 significant bit first. Other modules might use even different means of
2127 representing the numbers. See the respective module documentation for further
2130 Please note that Math::BigFloat does B<not> use the denoted library itself,
2131 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2134 use Math::BigInt lib => 'GMP';
2137 you can roll it all into one line:
2139 use Math::BigFloat lib => 'GMP';
2141 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details.
2143 =head2 Using Math::BigInt::Lite
2145 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2148 use Math::BigFloat with => 'Math::BigInt::Lite';
2150 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2151 can combine these if you want. For instance, you may want to use
2152 Math::BigInt objects in your main script, too.
2156 use Math::BigFloat with => 'Math::BigInt::Lite';
2158 Of course, you can combine this with the C<lib> parameter.
2161 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2163 If you want to use Math::BigInt's, too, simple add a Math::BigInt B<before>:
2167 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2169 Notice that the module with the last C<lib> will "win" and thus
2170 it's lib will be used if the lib is available:
2173 use Math::BigInt lib => 'Bar,Baz';
2174 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2176 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2177 words, Math::BigFloat will try to retain previously loaded libs when you
2178 don't specify it one.
2180 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2181 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2182 same as trying the latter load alone, except for the fact that Bar or Baz
2183 might be loaded needlessly in an intermidiate step
2185 The old way still works though:
2188 use Math::BigInt lib => 'Bar,Baz';
2191 But B<examples #3 and #4 are recommended> for usage.
2199 The following does not work yet:
2201 $m = $x->mantissa();
2202 $e = $x->exponent();
2203 $y = $m * ( 10 ** $e );
2204 print "ok\n" if $x == $y;
2208 There is no fmod() function yet.
2216 =item stringify, bstr()
2218 Both stringify and bstr() now drop the leading '+'. The old code would return
2219 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2220 reasoning and details.
2224 The following will probably not do what you expect:
2226 print $c->bdiv(123.456),"\n";
2228 It prints both quotient and reminder since print works in list context. Also,
2229 bdiv() will modify $c, so be carefull. You probably want to use
2231 print $c / 123.456,"\n";
2232 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2236 =item Modifying and =
2240 $x = Math::BigFloat->new(5);
2243 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2244 a second reference to the B<same> object and stores it in $y. Thus anything
2245 that modifies $x will modify $y, and vice versa.
2248 print "$x, $y\n"; # prints '10, 10'
2250 If you want a true copy of $x, use:
2254 See also the documentation in L<overload> regarding C<=>.
2258 C<bpow()> now modifies the first argument, unlike the old code which left
2259 it alone and only returned the result. This is to be consistent with
2260 C<badd()> etc. The first will modify $x, the second one won't:
2262 print bpow($x,$i),"\n"; # modify $x
2263 print $x->bpow($i),"\n"; # ditto
2264 print $x ** $i,"\n"; # leave $x alone
2270 This program is free software; you may redistribute it and/or modify it under
2271 the same terms as Perl itself.
2275 Mark Biggar, overloaded interface by Ilya Zakharevich.
2276 Completely rewritten by Tels http://bloodgate.com in 2001.