3 # The following hash values are internally used:
4 # _e: exponent (BigInt)
5 # _m: mantissa (absolute BigInt)
6 # sign: +,-,"NaN" if not a number
9 # _f: flags, used to signal MBI not to touch our private parts
10 # _cow: Copy-On-Write (NRY)
12 package Math::BigFloat;
17 use Math::BigInt qw/objectify/;
18 @ISA = qw( Exporter Math::BigInt);
19 # can not export bneg/babs since the are only in MBI
22 badd bmul bdiv bmod bnorm bsub
23 bgcd blcm bround bfround
24 bpow bnan bzero bfloor bceil
25 bacmp bstr binc bdec binf
26 is_odd is_even is_nan is_inf is_positive is_negative
32 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $rnd_mode/;
33 my $class = "Math::BigFloat";
36 '<=>' => sub { $_[2] ?
37 ref($_[0])->bcmp($_[1],$_[0]) :
38 ref($_[0])->bcmp($_[0],$_[1])},
39 'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
42 ##############################################################################
43 # global constants, flags and accessory
45 use constant MB_NEVER_ROUND => 0x0001;
49 # constant for easier life
52 # Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
58 # in case we call SUPER::->foo() and this wants to call modify()
59 # sub modify () { 0; }
63 my %methods = map { $_ => 1 }
64 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
65 fabs fneg fint fcmp fzero fnan finc fdec
68 sub method_valid { return exists $methods{$_[0]||''}; }
71 ##############################################################################
76 # create a new BigFloat object from a string or another bigfloat object.
79 # sign => sign (+/-), or "NaN"
83 my $wanted = shift; # avoid numify call by not using || here
84 return $class->bzero() if !defined $wanted; # default to 0
85 return $wanted->copy() if ref($wanted) eq $class;
87 my $round = shift; $round = 0 if !defined $round; # no rounding as default
88 my $self = {}; bless $self, $class;
89 # shortcut for bigints and its subclasses
90 if ((ref($wanted)) && (ref($wanted) ne $class))
92 $self->{_m} = $wanted->as_number(); # get us a bigint copy
93 $self->{_e} = Math::BigInt->new(0);
95 $self->{sign} = $wanted->sign();
96 return $self->bnorm();
99 # handle '+inf', '-inf' first
100 if ($wanted =~ /^[+-]inf$/)
102 $self->{_e} = Math::BigInt->new(0);
103 $self->{_m} = Math::BigInt->new(0);
104 $self->{sign} = $wanted;
105 return $self->bnorm();
107 #print "new string '$wanted'\n";
108 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
111 die "$wanted is not a number initialized to $class" if !$NaNOK;
112 $self->{_e} = Math::BigInt->new(0);
113 $self->{_m} = Math::BigInt->new(0);
114 $self->{sign} = $nan;
118 # make integer from mantissa by adjusting exp, then convert to bigint
119 $self->{_e} = Math::BigInt->new("$$es$$ev"); # exponent
120 $self->{_m} = Math::BigInt->new("$$mis$$miv$$mfv"); # create mantissa
121 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
122 $self->{_e} -= CORE::length($$mfv);
123 $self->{sign} = $self->{_m}->sign(); $self->{_m}->babs();
125 #print "$wanted => $self->{sign} $self->{value}\n";
126 $self->bnorm(); # first normalize
127 # if any of the globals is set, round to them and thus store them insid $self
128 $self->round($accuracy,$precision,$rnd_mode)
129 if defined $accuracy || defined $precision;
135 # create a bigfloat 'NaN', if given a BigFloat, set it to 'NaN'
137 $self = $class if !defined $self;
140 my $c = $self; $self = {}; bless $self, $c;
142 $self->{_m} = Math::BigInt->bzero();
143 $self->{_e} = Math::BigInt->bzero();
144 $self->{sign} = $nan;
150 # create a bigfloat '+-inf', if given a BigFloat, set it to '+-inf'
152 my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
154 $self = $class if !defined $self;
157 my $c = $self; $self = {}; bless $self, $c;
159 $self->{_m} = Math::BigInt->bzero();
160 $self->{_e} = Math::BigInt->bzero();
161 $self->{sign} = $sign.'inf';
167 # create a bigfloat '+-1', if given a BigFloat, set it to '+-1'
169 my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
171 $self = $class if !defined $self;
174 my $c = $self; $self = {}; bless $self, $c;
176 $self->{_m} = Math::BigInt->bone();
177 $self->{_e} = Math::BigInt->bzero();
178 $self->{sign} = $sign;
184 # create a bigfloat '+0', if given a BigFloat, set it to 0
186 $self = $class if !defined $self;
189 my $c = $self; $self = {}; bless $self, $c;
191 $self->{_m} = Math::BigInt->bzero();
192 $self->{_e} = Math::BigInt->bone();
197 ##############################################################################
198 # string conversation
202 # (ref to BFLOAT or num_str ) return num_str
203 # Convert number from internal format to (non-scientific) string format.
204 # internal format is always normalized (no leading zeros, "-0" => "+0")
205 my ($self,$x) = objectify(1,@_);
207 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
208 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
209 if ($x->{sign} !~ /^[+-]$/)
211 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
215 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
217 my $not_zero = !$x->is_zero();
220 $es = $x->{_m}->bstr();
221 $len = CORE::length($es);
222 if (!$x->{_e}->is_zero())
224 # $es = $x->{sign}.$es if $x->{sign} eq '-';
228 if ($x->{_e}->sign() eq '-')
231 if ($x->{_e} <= -$len)
233 # print "style: 0.xxxx\n";
234 my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
235 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
239 # print "insert '.' at $x->{_e} in '$es'\n";
240 substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
246 $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
250 $es = $x->{sign}.$es if $x->{sign} eq '-';
251 # if set accuracy or precision, pad with zeros
252 if ((defined $x->{_a}) && ($not_zero))
254 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
255 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
256 $zeros = $x->{_a} - $len if $cad != $len;
257 #print "acc padd $x->{_a} $zeros (len $len cad $cad)\n";
258 $es .= $dot.'0' x $zeros if $zeros > 0;
260 elsif ($x->{_p} || 0 < 0)
262 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
263 my $zeros = -$x->{_p} + $cad;
264 #print "pre padd $x->{_p} $zeros (len $len cad $cad)\n";
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) = objectify(1,@_);
277 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
278 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
279 if ($x->{sign} !~ /^[+-]$/)
281 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
284 my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
286 return $x->{_m}->bstr().$sep.$x->{_e}->bstr();
291 # Make a number from a BigFloat object
292 # simple return string and let Perl's atoi()/atof() handle the rest
293 my ($self,$x) = objectify(1,@_);
297 ##############################################################################
298 # public stuff (usually prefixed with "b")
300 # really? Just for exporting them is not what I had in mind
303 # $class->SUPER::babs($class,@_);
307 # $class->SUPER::bneg($class,@_);
311 # todo: this must be overwritten and return NaN for non-integer values
312 # band(), bior(), bxor(), too
315 # $class->SUPER::bnot($class,@_);
320 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
321 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
322 my ($self,$x,$y) = objectify(2,@_);
324 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
326 # handle +-inf and NaN
327 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
328 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
329 return +1 if $x->{sign} eq '+inf';
330 return -1 if $x->{sign} eq '-inf';
331 return -1 if $y->{sign} eq '+inf';
332 return +1 if $y->{sign} eq '-inf';
335 # check sign for speed first
336 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
337 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
340 my $xz = $x->is_zero();
341 my $yz = $y->is_zero();
342 return 0 if $xz && $yz; # 0 <=> 0
343 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
344 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
346 # adjust so that exponents are equal
347 my $lxm = $x->{_m}->length();
348 my $lym = $y->{_m}->length();
349 my $lx = $lxm + $x->{_e};
350 my $ly = $lym + $y->{_e};
351 # print "x $x y $y lx $lx ly $ly\n";
352 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
353 # print "$l $x->{sign}\n";
354 return $l <=> 0 if $l != 0;
356 # lengths (corrected by exponent) are equal
357 # so make mantissa euqal length by padding with zero (shift left)
358 my $diff = $lxm - $lym;
359 my $xm = $x->{_m}; # not yet copy it
363 $ym = $y->{_m}->copy()->blsft($diff,10);
367 $xm = $x->{_m}->copy()->blsft(-$diff,10);
369 my $rc = $xm->bcmp($ym);
370 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
376 # Compares 2 values, ignoring their signs.
377 # Returns one of undef, <0, =0, >0. (suitable for sort)
378 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
379 my ($self,$x,$y) = objectify(2,@_);
380 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
382 # signs are ignored, so check length
383 # length(x) is length(m)+e aka length of non-fraction part
384 # the longer one is bigger
385 my $l = $x->length() - $y->length();
387 return $l if $l != 0;
388 #print "equal lengths\n";
390 # if both are equal long, make full compare
391 # first compare only the mantissa
392 # if mantissa are equal, compare fractions
394 return $x->{_m} <=> $y->{_m} || $x->{_e} <=> $y->{_e};
399 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
400 # return result as BFLOAT
401 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
403 # inf and NaN handling
404 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
407 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
409 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
411 # + and + => +, - and - => -, + and - => 0, - and + => 0
412 return $x->bzero() if $x->{sign} ne $y->{sign};
415 # +-inf + something => +inf
416 # something +-inf => +-inf
417 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
421 # speed: no add for 0+y or x+0
422 return $x if $y->is_zero(); # x+0
423 if ($x->is_zero()) # 0+y
425 # make copy, clobbering up x (modify in place!)
426 $x->{_e} = $y->{_e}->copy();
427 $x->{_m} = $y->{_m}->copy();
428 $x->{sign} = $y->{sign} || $nan;
429 return $x->round($a,$p,$r,$y);
432 # take lower of the two e's and adapt m1 to it to match m2
433 my $e = $y->{_e}; $e = Math::BigInt::bzero() if !defined $e; # if no BFLOAT
435 my $add = $y->{_m}->copy();
439 #print "\$x->{_m}: $x->{_m} ";
440 #print "\$x->{_e}: $x->{_e}\n";
441 my $e1 = $e->copy()->babs();
442 $x->{_m} *= (10 ** $e1);
443 $x->{_e} += $e; # need the sign of e
444 #$x->{_m} += $y->{_m};
445 #print "\$x->{_m}: $x->{_m} ";
446 #print "\$x->{_e}: $x->{_e}\n";
451 #print "\$x->{_m}: $x->{_m} \$y->{_m}: $y->{_m} \$e: $e ",ref($e),"\n";
453 #$x->{_m} += $y->{_m} * (10 ** $e);
454 #print "\$x->{_m}: $x->{_m}\n";
456 # else: both e are same, so leave them
457 #print "badd $x->{sign}$x->{_m} + $y->{sign}$add\n";
459 $x->{_m}->{sign} = $x->{sign};
460 $add->{sign} = $y->{sign};
464 $x->{sign} = $x->{_m}->{sign};
465 $x->{_m}->{sign} = '+';
466 #$x->bnorm(); # delete trailing zeros
467 return $x->round($a,$p,$r,$y);
472 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
473 # subtract second arg from first, modify first
474 my ($self,$x,$y) = objectify(2,@_);
476 $x->badd($y->bneg()); # badd does not leave internal zeros
477 $y->bneg(); # refix y, assumes no one reads $y in between
478 return $x; # badd() already normalized and rounded
483 # increment arg by one
484 my ($self,$x,$a,$p,$r) = objectify(1,@_);
485 $x->badd($self->_one())->round($a,$p,$r);
490 # decrement arg by one
491 my ($self,$x,$a,$p,$r) = objectify(1,@_);
492 $x->badd($self->_one('-'))->round($a,$p,$r);
497 # (BINT or num_str, BINT or num_str) return BINT
498 # does not modify arguments, but returns new object
499 # Lowest Common Multiplicator
501 my ($self,@arg) = objectify(0,@_);
502 my $x = $self->new(shift @arg);
503 while (@arg) { $x = _lcm($x,shift @arg); }
509 # (BINT or num_str, BINT or num_str) return BINT
510 # does not modify arguments, but returns new object
511 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
513 my ($self,@arg) = objectify(0,@_);
514 my $x = $self->new(shift @arg);
515 while (@arg) { $x = _gcd($x,shift @arg); }
521 # return true if arg (BINT or num_str) is zero (array '+', '0')
522 my $x = shift; $x = $class->new($x) unless ref $x;
524 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
530 # return true if arg (BINT or num_str) is +1 (array '+', '1')
531 # or -1 if signis given
532 my $x = shift; $x = $class->new($x) unless ref $x;
533 #my ($self,$x) = objectify(1,@_);
534 my $sign = $_[2] || '+';
535 return ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
540 # return true if arg (BINT or num_str) is odd or false if even
541 my $x = shift; $x = $class->new($x) unless ref $x;
542 #my ($self,$x) = objectify(1,@_);
544 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
545 return ($x->{_e}->is_zero() && $x->{_m}->is_odd());
550 # return true if arg (BINT or num_str) is even or false if odd
551 my $x = shift; $x = $class->new($x) unless ref $x;
552 #my ($self,$x) = objectify(1,@_);
554 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
555 return 1 if $x->{_m}->is_zero(); # 0e1 is even
556 return ($x->{_e}->is_zero() && $x->{_m}->is_even()); # 123.45 is never
561 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
562 # (BINT or num_str, BINT or num_str) return BINT
563 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
565 # print "mbf bmul $x->{_m}e$x->{_e} $y->{_m}e$y->{_e}\n";
566 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
569 return $x->bzero() if $x->is_zero() || $y->is_zero();
571 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
573 # result will always be +-inf:
574 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
575 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
576 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
577 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
578 return $x->binf('-');
581 # aEb * cEd = (a*c)E(b+d)
582 $x->{_m} = $x->{_m} * $y->{_m};
583 #print "m: $x->{_m}\n";
584 $x->{_e} = $x->{_e} + $y->{_e};
585 #print "e: $x->{_m}\n";
587 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
588 #print "s: $x->{sign}\n";
590 return $x->round($a,$p,$r,$y);
595 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
596 # (BFLOAT,BFLOAT) (quo,rem) or BINT (only rem)
597 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
599 # x / +-inf => 0, reminder x
600 return wantarray ? ($x->bzero(),$x->copy()) : $x->bzero()
601 if $y->{sign} =~ /^[+-]inf$/;
603 # NaN if x == NaN or y == NaN or x==y==0
604 return wantarray ? ($x->bnan(),bnan()) : $x->bnan()
605 if (($x->is_nan() || $y->is_nan()) ||
606 ($x->is_zero() && $y->is_zero()));
608 # 5 / 0 => +inf, -6 / 0 => -inf
610 ? ($x->binf($x->{sign}),$self->bnan()) : $x->binf($x->{sign})
611 if ($x->{sign} =~ /^[+-]$/ && $y->is_zero());
613 $y = $class->new($y) if ref($y) ne $class; # promote bigints
615 # print "mbf bdiv $x ",ref($x)," ",$y," ",ref($y),"\n";
616 # we need to limit the accuracy to protect against overflow
617 my ($scale) = $x->_scale_a($accuracy,$rnd_mode,$a,$r); # ignore $p
621 # simulate old behaviour
622 $scale = $div_scale+1; # one more for proper riund
623 $a = $div_scale; # and round to it
624 $fallback = 1; # to clear a/p afterwards
626 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
627 $scale = $lx if $lx > $scale;
628 $scale = $ly if $ly > $scale;
629 #print "scale $scale $lx $ly\n";
630 my $diff = $ly - $lx;
631 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
633 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
635 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
637 # check for / +-1 ( +/- 1E0)
640 return wantarray ? ($x,$self->bzero()) : $x;
643 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
644 #print "self: $self x: $x ref(x) ", ref($x)," m: $x->{_m}\n";
645 # my $scale_10 = 10 ** $scale; $x->{_m}->bmul($scale_10);
646 $x->{_m}->blsft($scale,10);
647 #print "m: $x->{_m} $y->{_m}\n";
648 $x->{_m}->bdiv( $y->{_m} ); # a/c
649 #print "m: $x->{_m}\n";
650 #print "e: $x->{_e} $y->{_e}",$scale,"\n";
651 $x->{_e}->bsub($y->{_e}); # b-d
652 #print "e: $x->{_e}\n";
653 $x->{_e}->bsub($scale); # correct for 10**scale
654 #print "after div: m: $x->{_m} e: $x->{_e}\n";
655 $x->bnorm(); # remove trailing 0's
656 #print "after div: m: $x->{_m} e: $x->{_e}\n";
657 $x->round($a,$p,$r); # then round accordingly
660 # clear a/p after round, since user did not request it
667 my $rem = $x->copy();
668 $rem->bmod($y,$a,$p,$r);
671 # clear a/p after round, since user did not request it
682 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
683 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
685 return $x->bnan() if ($x->{sign} eq $nan || $y->is_nan() || $y->is_zero());
686 return $x->bzero() if $y->is_one();
688 # XXX tels: not done yet
689 return $x->round($a,$p,$r,$y);
694 # calculate square root; this should probably
695 # use a different test to see whether the accuracy we want is...
696 my ($self,$x,$a,$p,$r) = objectify(1,@_);
698 return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
699 return $x if $x->{sign} eq '+inf'; # +inf
700 return $x if $x->is_zero() || $x == 1;
702 # we need to limit the accuracy to protect against overflow
703 my ($scale) = $x->_scale_a($accuracy,$rnd_mode,$a,$r); # ignore $p
707 # simulate old behaviour
708 $scale = $div_scale+1; # one more for proper riund
709 $a = $div_scale; # and round to it
710 $fallback = 1; # to clear a/p afterwards
712 my $lx = $x->{_m}->length();
713 $scale = $lx if $scale < $lx;
714 my $e = Math::BigFloat->new("1E-$scale"); # make test variable
715 return $x->bnan() if $e->sign() eq 'NaN';
717 # start with some reasonable guess
718 #$x *= 10 ** ($len - $org->{_e}); $x /= 2; # !?!?
721 my $gs = Math::BigFloat->new('1'. ('0' x $lx));
723 # print "first guess: $gs (x $x) scale $scale\n";
727 my $two = Math::BigFloat->new(2);
728 $x = Math::BigFloat->new($x) if ref($x) ne $class; # promote BigInts
732 return $x->bnan() if $gs->is_zero();
733 $r = $y->copy(); $r->bdiv($gs,$scale);
735 $x->bdiv($two,$scale);
736 $diff = $x->copy()->bsub($gs)->babs();
742 # clear a/p after round, since user did not request it
751 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
752 # compute power of two numbers, second arg is used as integer
753 # modifies first argument
755 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
757 return $x if $x->{sign} =~ /^[+-]inf$/;
758 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
759 return $x->bone() if $y->is_zero();
760 return $x if $x->is_one() || $y->is_one();
761 my $y1 = $y->as_number(); # make bigint
764 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
765 return $y1->is_odd() ? $x : $x->babs(1);
767 return $x if $x->is_zero() && $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
768 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
769 return $x->binf() if $x->is_zero() && $y->{sign} eq '-';
771 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
775 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
777 if ($y->{sign} eq '-')
779 # modify $x in place!
780 my $z = $x->copy(); $x->bzero()->binc();
781 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
783 return $x->round($a,$p,$r,$y);
786 ###############################################################################
791 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
792 # $n == 0 means round to integer
793 # expects and returns normalized numbers!
794 my $x = shift; $x = $class->new($x) unless ref $x;
796 return $x if $x->modify('bfround');
798 my ($scale,$mode) = $x->_scale_p($precision,$rnd_mode,@_);
799 return $x if !defined $scale; # no-op
801 # never round a 0, +-inf, NaN
802 return $x if $x->{sign} !~ /^[+-]$/ || $x->is_zero();
803 # print "MBF bfround $x to scale $scale mode $mode\n";
807 # print "bfround scale $scale e $x->{_e}\n";
808 # round right from the '.'
809 return $x if $x->{_e} >= 0; # nothing to round
810 $scale = -$scale; # positive for simplicity
811 my $len = $x->{_m}->length(); # length of mantissa
812 my $dad = -$x->{_e}; # digits after dot
813 my $zad = 0; # zeros after dot
814 $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
815 # print "scale $scale dad $dad zad $zad len $len\n";
817 # number bsstr len zad dad
819 # 0.0123 123e-4 3 1 4
822 # 1.2345 12345e-4 5 0 4
824 # do not round after/right of the $dad
825 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
827 # round to zero if rounding inside the $zad, but not for last zero like:
828 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
833 if ($scale == $zad) # for 0.006, scale -2 and trunc
839 # adjust round-point to be inside mantissa
842 $scale = $scale-$zad;
846 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
847 $scale = $dbd+$scale;
850 # print "round to $x->{_m} to $scale\n";
854 # 123 => 100 means length(123) = 3 - $scale (2) => 1
856 # calculate digits before dot
857 my $dbt = $x->{_m}->length(); $dbt += $x->{_e} if $x->{_e}->sign() eq '-';
858 if (($scale > $dbt) && ($dbt < 0))
860 # if not enough digits before dot, round to zero
863 if (($scale >= 0) && ($dbt == 0))
865 # 0.49->bfround(1): scale == 1, dbt == 0: => 0.0
866 # 0.51->bfround(0): scale == 0, dbt == 0: => 1.0
867 # 0.5->bfround(0): scale == 0, dbt == 0: => 0
868 # 0.05->bfround(0): scale == 0, dbt == 0: => 0
869 # print "$scale $dbt $x->{_m}\n";
870 $scale = -$x->{_m}->length();
874 # correct by subtracting scale
875 $scale = $dbt - $scale;
879 $scale = $x->{_m}->length() - $scale;
882 # print "using $scale for $x->{_m} with '$mode'\n";
883 # pass sign to bround for rounding modes '+inf' and '-inf'
884 $x->{_m}->{sign} = $x->{sign};
885 $x->{_m}->bround($scale,$mode);
886 $x->{_m}->{sign} = '+'; # fix sign back
892 # accuracy: preserve $N digits, and overwrite the rest with 0's
893 my $x = shift; $x = $class->new($x) unless ref $x;
894 my ($scale,$mode) = $x->_scale_a($accuracy,$rnd_mode,@_);
895 return $x if !defined $scale; # no-op
897 return $x if $x->modify('bround');
899 # print "bround $scale $mode\n";
900 # 0 => return all digits, scale < 0 makes no sense
901 return $x if ($scale <= 0);
902 # never round a 0, +-inf, NaN
903 return $x if $x->{sign} !~ /^[+-]$/ || $x->is_zero();
905 # if $e longer than $m, we have 0.0000xxxyyy style number, and must
906 # subtract the delta from scale, to simulate keeping the zeros
907 # -5 +5 => 1; -10 +5 => -4
908 my $delta = $x->{_e} + $x->{_m}->length() + 1;
909 # removed by tlr, since causes problems with fraction tests:
910 # $scale += $delta if $delta < 0;
912 # if we should keep more digits than the mantissa has, do nothing
913 return $x if $x->{_m}->length() <= $scale;
915 # pass sign to bround for '+inf' and '-inf' rounding modes
916 $x->{_m}->{sign} = $x->{sign};
917 $x->{_m}->bround($scale,$mode); # round mantissa
918 $x->{_m}->{sign} = '+'; # fix sign back
919 $x->bnorm(); # del trailing zeros gen. by bround()
924 # return integer less or equal then $x
925 my ($self,$x,$a,$p,$r) = objectify(1,@_);
927 return $x if $x->modify('bfloor');
929 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
931 # if $x has digits after dot
932 if ($x->{_e}->{sign} eq '-')
934 $x->{_m}->brsft(-$x->{_e},10);
936 $x-- if $x->{sign} eq '-';
938 return $x->round($a,$p,$r);
943 # return integer greater or equal then $x
944 my ($self,$x,$a,$p,$r) = objectify(1,@_);
946 return $x if $x->modify('bceil');
947 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
949 # if $x has digits after dot
950 if ($x->{_e}->{sign} eq '-')
952 $x->{_m}->brsft(-$x->{_e},10);
954 $x++ if $x->{sign} eq '+';
956 return $x->round($a,$p,$r);
959 ###############################################################################
963 # going trough AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
968 # make fxxx and bxxx work
970 my $name = $AUTOLOAD;
972 $name =~ s/.*:://; # split package
974 if (!method_valid($name))
978 #&{$class."::SUPER->$name"}(@_);
979 # delayed load of Carp and avoid recursion
981 Carp::croak ("Can't call $class\-\>$name, not a valid method");
984 my $bname = $name; $bname =~ s/^f/b/;
985 *{$class."\:\:$name"} = \&$bname;
991 # return a copy of the exponent
993 $self = $class->new($self) unless ref $self;
995 return bnan() if $self->is_nan();
996 return $self->{_e}->copy();
1001 # return a copy of the mantissa
1003 $self = $class->new($self) unless ref $self;
1005 return bnan() if $self->is_nan();
1006 my $m = $self->{_m}->copy(); # faster than going via bstr()
1007 $m->bneg() if $self->{sign} eq '-';
1014 # return a copy of both the exponent and the mantissa
1016 $self = $class->new($self) unless ref $self;
1018 return (bnan(),bnan()) if $self->is_nan();
1019 my $m = $self->{_m}->copy(); # faster than going via bstr()
1020 $m->bneg() if $self->{sign} eq '-';
1021 return ($m,$self->{_e}->copy());
1024 ##############################################################################
1025 # private stuff (internal use only)
1029 # internal speedup, set argument to 1, or create a +/- 1
1030 my $self = shift; $self = ref($self) if ref($self);
1031 my $x = {}; bless $x, $self;
1032 $x->{_m} = Math::BigInt->new(1);
1033 $x->{_e} = Math::BigInt->new(0);
1034 $x->{sign} = shift || '+';
1041 #print "import $self\n";
1042 for ( my $i = 0; $i < @_ ; $i++ )
1044 if ( $_[$i] eq ':constant' )
1046 # this rest causes overlord er load to step in
1047 # print "overload @_\n";
1048 overload::constant float => sub { $self->new(shift); };
1049 splice @_, $i, 1; last;
1052 # any non :constant stuff is handled by our parent, Exporter
1053 # even if @_ is empty, to give it a chance
1054 #$self->SUPER::import(@_); # does not work (would call MBI)
1055 $self->export_to_level(1,$self,@_); # need this instead
1060 # adjust m and e so that m is smallest possible
1061 # round number according to accuracy and precision settings
1064 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
1066 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
1069 $x->{_m}->brsft($zeros,10); $x->{_e} += $zeros;
1071 # for something like 0Ey, set y to 1
1072 $x->{sign} = '+', $x->{_e}->bzero()->binc() if $x->{_m}->is_zero();
1073 $x->{_m}->{_f} = MB_NEVER_ROUND;
1074 $x->{_e}->{_f} = MB_NEVER_ROUND;
1075 return $x; # MBI bnorm is no-op
1078 ##############################################################################
1079 # internal calculation routines
1083 # return a bigint representation of this BigFloat number
1084 my ($self,$x) = objectify(1,@_);
1087 if ($x->{_e}->is_zero())
1089 $z = $x->{_m}->copy();
1090 $z->{sign} = $x->{sign};
1093 $z = $x->{_m}->copy();
1096 $z->brsft(-$x->{_e},10);
1100 $z->blsft($x->{_e},10);
1102 $z->{sign} = $x->{sign};
1108 my $x = shift; $x = $class->new($x) unless ref $x;
1110 my $len = $x->{_m}->length();
1111 $len += $x->{_e} if $x->{_e}->sign() eq '+';
1114 my $t = Math::BigInt::bzero();
1115 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1126 Math::BigFloat - Arbitrary size floating point math package
1133 $x = Math::BigInt->new($str); # defaults to 0
1134 $nan = Math::BigInt->bnan(); # create a NotANumber
1135 $zero = Math::BigInt->bzero();# create a "+0"
1138 $x->is_zero(); # return whether arg is zero or not
1139 $x->is_nan(); # return whether arg is NaN or not
1140 $x->is_one(); # true if arg is +1
1141 $x->is_one('-'); # true if arg is -1
1142 $x->is_odd(); # true if odd, false for even
1143 $x->is_even(); # true if even, false for odd
1144 $x->is_positive(); # true if >= 0
1145 $x->is_negative(); # true if < 0
1146 $x->is_inf(sign) # true if +inf or -inf (sign default '+')
1147 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
1148 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
1149 $x->sign(); # return the sign, either +,- or NaN
1151 # The following all modify their first argument:
1154 $x->bzero(); # set $i to 0
1155 $x->bnan(); # set $i to NaN
1157 $x->bneg(); # negation
1158 $x->babs(); # absolute value
1159 $x->bnorm(); # normalize (no-op)
1160 $x->bnot(); # two's complement (bit wise not)
1161 $x->binc(); # increment x by 1
1162 $x->bdec(); # decrement x by 1
1164 $x->badd($y); # addition (add $y to $x)
1165 $x->bsub($y); # subtraction (subtract $y from $x)
1166 $x->bmul($y); # multiplication (multiply $x by $y)
1167 $x->bdiv($y); # divide, set $i to quotient
1168 # return (quo,rem) or quo if scalar
1170 $x->bmod($y); # modulus
1171 $x->bpow($y); # power of arguments (a**b)
1172 $x->blsft($y); # left shift
1173 $x->brsft($y); # right shift
1174 # return (quo,rem) or quo if scalar
1176 $x->band($y); # bit-wise and
1177 $x->bior($y); # bit-wise inclusive or
1178 $x->bxor($y); # bit-wise exclusive or
1179 $x->bnot(); # bit-wise not (two's complement)
1181 $x->bround($N); # accuracy: preserver $N digits
1182 $x->bfround($N); # precision: round to the $Nth digit
1184 # The following do not modify their arguments:
1186 bgcd(@values); # greatest common divisor
1187 blcm(@values); # lowest common multiplicator
1189 $x->bstr(); # return string
1190 $x->bsstr(); # return string in scientific notation
1192 $x->exponent(); # return exponent as BigInt
1193 $x->mantissa(); # return mantissa as BigInt
1194 $x->parts(); # return (mantissa,exponent) as BigInt
1196 $x->length(); # number of digits (w/o sign and '.')
1197 ($l,$f) = $x->length(); # number of digits, and length of fraction
1201 All operators (inlcuding basic math operations) are overloaded if you
1202 declare your big floating point numbers as
1204 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1206 Operations with overloaded operators preserve the arguments, which is
1207 exactly what you expect.
1209 =head2 Canonical notation
1211 Input to these routines are either BigFloat objects, or strings of the
1212 following four forms:
1226 C</^[+-]\d+E[+-]?\d+$/>
1230 C</^[+-]\d*\.\d+E[+-]?\d+$/>
1234 all with optional leading and trailing zeros and/or spaces. Additonally,
1235 numbers are allowed to have an underscore between any two digits.
1237 Empty strings as well as other illegal numbers results in 'NaN'.
1239 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1240 are always stored in normalized form. On a string, it creates a BigFloat
1245 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1247 The string output will always have leading and trailing zeros stripped and drop
1248 a plus sign. C<bstr()> will give you always the form with a decimal point,
1249 while C<bsstr()> (for scientific) gives you the scientific notation.
1251 Input bstr() bsstr()
1253 ' -123 123 123' '-123123123' '-123123123E0'
1254 '00.0123' '0.0123' '123E-4'
1255 '123.45E-2' '1.2345' '12345E-4'
1256 '10E+3' '10000' '1E4'
1258 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1259 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1260 return either undef, <0, 0 or >0 and are suited for sort.
1262 Actual math is done by using BigInts to represent the mantissa and exponent.
1263 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1264 represent the result when input arguments are not numbers, as well as
1265 the result of dividing by zero.
1267 =head2 C<mantissa()>, C<exponent()> and C<parts()>
1269 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1270 as BigInts such that:
1272 $m = $x->mantissa();
1273 $e = $x->exponent();
1274 $y = $m * ( 10 ** $e );
1275 print "ok\n" if $x == $y;
1277 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
1279 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
1281 Currently the mantissa is reduced as much as possible, favouring higher
1282 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
1283 This might change in the future, so do not depend on it.
1285 =head2 Accuracy vs. Precision
1287 See also: L<Rounding|Rounding>.
1289 Math::BigFloat supports both precision and accuracy. (here should follow
1290 a short description of both).
1292 Precision: digits after the '.', laber, schwad
1293 Accuracy: Significant digits blah blah
1295 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
1296 a operation consumes all resources, each operation produces no more than
1297 C<Math::BigFloat::precision()> digits.
1299 In case the result of one operation has more precision than specified,
1300 it is rounded. The rounding mode taken is either the default mode, or the one
1301 supplied to the operation after the I<scale>:
1303 $x = Math::BigFloat->new(2);
1304 Math::BigFloat::precision(5); # 5 digits max
1305 $y = $x->copy()->bdiv(3); # will give 0.66666
1306 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1307 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
1308 Math::BigFloat::round_mode('zero');
1309 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1315 =item ffround ( +$scale )
1317 Rounds to the $scale'th place left from the '.', counting from the dot.
1318 The first digit is numbered 1.
1320 =item ffround ( -$scale )
1322 Rounds to the $scale'th place right from the '.', counting from the dot.
1326 Rounds to an integer.
1328 =item fround ( +$scale )
1330 Preserves accuracy to $scale digits from the left (aka significant digits)
1331 and pads the rest with zeros. If the number is between 1 and -1, the
1332 significant digits count from the first non-zero after the '.'
1334 =item fround ( -$scale ) and fround ( 0 )
1336 These are effetively no-ops.
1340 All rounding functions take as a second parameter a rounding mode from one of
1341 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
1343 The default rounding mode is 'even'. By using
1344 C<< Math::BigFloat::round_mode($rnd_mode); >> you can get and set the default
1345 mode for subsequent rounding. The usage of C<$Math::BigFloat::$rnd_mode> is
1346 no longer supported.
1347 The second parameter to the round functions then overrides the default
1350 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
1351 'trunc' as rounding mode to make it equivalent to:
1356 You can override this by passing the desired rounding mode as parameter to
1359 $x = Math::BigFloat->new(2.5);
1360 $y = $x->as_number('odd'); # $y = 3
1366 =head1 Autocreating constants
1368 After C<use Math::BigFloat ':constant'> all the floating point constants
1369 in the given scope are converted to C<Math::BigFloat>. This conversion
1370 happens at compile time.
1374 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
1376 prints the value of C<2E-100>. Note that without conversion of
1377 constants the expression 2E-100 will be calculated as normal floating point
1386 The following does not work yet:
1388 $m = $x->mantissa();
1389 $e = $x->exponent();
1390 $y = $m * ( 10 ** $e );
1391 print "ok\n" if $x == $y;
1395 There is no fmod() function yet.
1403 =item stringify, bstr()
1405 Both stringify and bstr() now drop the leading '+'. The old code would return
1406 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
1407 reasoning and details.
1411 The following will probably not do what you expect:
1413 print $c->bdiv(123.456),"\n";
1415 It prints both quotient and reminder since print works in list context. Also,
1416 bdiv() will modify $c, so be carefull. You probably want to use
1418 print $c / 123.456,"\n";
1419 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
1423 =item Modifying and =
1427 $x = Math::BigFloat->new(5);
1430 It will not do what you think, e.g. making a copy of $x. Instead it just makes
1431 a second reference to the B<same> object and stores it in $y. Thus anything
1432 that modifies $x will modify $y, and vice versa.
1435 print "$x, $y\n"; # prints '10, 10'
1437 If you want a true copy of $x, use:
1441 See also the documentation in L<overload> regarding C<=>.
1445 C<bpow()> now modifies the first argument, unlike the old code which left
1446 it alone and only returned the result. This is to be consistent with
1447 C<badd()> etc. The first will modify $x, the second one won't:
1449 print bpow($x,$i),"\n"; # modify $x
1450 print $x->bpow($i),"\n"; # ditto
1451 print $x ** $i,"\n"; # leave $x alone
1457 This program is free software; you may redistribute it and/or modify it under
1458 the same terms as Perl itself.
1462 Mark Biggar, overloaded interface by Ilya Zakharevich.
1463 Completely rewritten by Tels http://bloodgate.com in 2001.