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";
38 $class->bcmp($_[1],$_[0]) :
39 $class->bcmp($_[0],$_[1])},
40 'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
43 ##############################################################################
44 # global constants, flags and accessory
46 use constant MB_NEVER_ROUND => 0x0001;
50 # constant for easier life
53 # Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
59 # in case we call SUPER::->foo() and this wants to call modify()
60 # sub modify () { 0; }
64 my %methods = map { $_ => 1 }
65 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
66 fabs fneg fint fcmp fzero fnan finc fdec
69 sub method_valid { return exists $methods{$_[0]||''}; }
72 ##############################################################################
77 # create a new BigFloat object from a string or another bigfloat object.
80 # sign => sign (+/-), or "NaN"
84 my $wanted = shift; # avoid numify call by not using || here
85 return $class->bzero() if !defined $wanted; # default to 0
86 return $wanted->copy() if ref($wanted) eq $class;
88 my $round = shift; $round = 0 if !defined $round; # no rounding as default
89 my $self = {}; bless $self, $class;
90 # shortcut for bigints and its subclasses
91 if ((ref($wanted)) && (ref($wanted) ne $class))
93 $self->{_m} = $wanted->as_number(); # get us a bigint copy
94 $self->{_e} = Math::BigInt->new(0);
96 $self->{sign} = $wanted->sign();
97 return $self->bnorm();
100 # handle '+inf', '-inf' first
101 if ($wanted =~ /^[+-]inf$/)
103 $self->{_e} = Math::BigInt->new(0);
104 $self->{_m} = Math::BigInt->new(0);
105 $self->{sign} = $wanted;
106 return $self->bnorm();
108 #print "new string '$wanted'\n";
109 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
112 die "$wanted is not a number initialized to $class" if !$NaNOK;
113 $self->{_e} = Math::BigInt->new(0);
114 $self->{_m} = Math::BigInt->new(0);
115 $self->{sign} = $nan;
119 # make integer from mantissa by adjusting exp, then convert to bigint
120 $self->{_e} = Math::BigInt->new("$$es$$ev"); # exponent
121 $self->{_m} = Math::BigInt->new("$$mis$$miv$$mfv"); # create mantissa
122 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
123 $self->{_e} -= CORE::length($$mfv);
124 $self->{sign} = $self->{_m}->sign(); $self->{_m}->babs();
126 #print "$wanted => $self->{sign} $self->{value}\n";
127 $self->bnorm(); # first normalize
128 # if any of the globals is set, round to them and thus store them insid $self
129 $self->round($accuracy,$precision,$rnd_mode)
130 if defined $accuracy || defined $precision;
136 # create a bigfloat 'NaN', if given a BigFloat, set it to 'NaN'
138 $self = $class if !defined $self;
141 my $c = $self; $self = {}; bless $self, $c;
143 $self->{_m} = Math::BigInt->bzero();
144 $self->{_e} = Math::BigInt->bzero();
145 $self->{sign} = $nan;
151 # create a bigfloat '+-inf', if given a BigFloat, set it to '+-inf'
153 my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
155 $self = $class if !defined $self;
158 my $c = $self; $self = {}; bless $self, $c;
160 $self->{_m} = Math::BigInt->bzero();
161 $self->{_e} = Math::BigInt->bzero();
162 $self->{sign} = $sign.'inf';
168 # create a bigfloat '+-1', if given a BigFloat, set it to '+-1'
170 my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
172 $self = $class if !defined $self;
175 my $c = $self; $self = {}; bless $self, $c;
177 $self->{_m} = Math::BigInt->bone();
178 $self->{_e} = Math::BigInt->bzero();
179 $self->{sign} = $sign;
185 # create a bigfloat '+0', if given a BigFloat, set it to 0
187 $self = $class if !defined $self;
190 my $c = $self; $self = {}; bless $self, $c;
192 $self->{_m} = Math::BigInt->bzero();
193 $self->{_e} = Math::BigInt->bone();
198 ##############################################################################
199 # string conversation
203 # (ref to BFLOAT or num_str ) return num_str
204 # Convert number from internal format to (non-scientific) string format.
205 # internal format is always normalized (no leading zeros, "-0" => "+0")
206 my ($self,$x) = objectify(1,@_);
208 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
209 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
210 if ($x->{sign} !~ /^[+-]$/)
212 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
216 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
218 my $not_zero = !$x->is_zero();
221 $es = $x->{_m}->bstr();
222 $len = CORE::length($es);
223 if (!$x->{_e}->is_zero())
225 # $es = $x->{sign}.$es if $x->{sign} eq '-';
229 if ($x->{_e}->sign() eq '-')
232 if ($x->{_e} <= -$len)
234 # print "style: 0.xxxx\n";
235 my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
236 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
240 # print "insert '.' at $x->{_e} in '$es'\n";
241 substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
247 $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
251 $es = $x->{sign}.$es if $x->{sign} eq '-';
252 # if set accuracy or precision, pad with zeros
253 if ((defined $x->{_a}) && ($not_zero))
255 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
256 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
257 $zeros = $x->{_a} - $len if $cad != $len;
258 #print "acc padd $x->{_a} $zeros (len $len cad $cad)\n";
259 $es .= $dot.'0' x $zeros if $zeros > 0;
261 elsif ($x->{_p} || 0 < 0)
263 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
264 my $zeros = -$x->{_p} + $cad;
265 #print "pre padd $x->{_p} $zeros (len $len cad $cad)\n";
266 $es .= $dot.'0' x $zeros if $zeros > 0;
273 # (ref to BFLOAT or num_str ) return num_str
274 # Convert number from internal format to scientific string format.
275 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
276 my ($self,$x) = objectify(1,@_);
278 #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
279 #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
280 if ($x->{sign} !~ /^[+-]$/)
282 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
285 my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
287 return $x->{_m}->bstr().$sep.$x->{_e}->bstr();
292 # Make a number from a BigFloat object
293 # simple return string and let Perl's atoi()/atof() handle the rest
294 my ($self,$x) = objectify(1,@_);
298 ##############################################################################
299 # public stuff (usually prefixed with "b")
301 # really? Just for exporting them is not what I had in mind
304 # $class->SUPER::babs($class,@_);
308 # $class->SUPER::bneg($class,@_);
312 # todo: this must be overwritten and return NaN for non-integer values
313 # band(), bior(), bxor(), too
316 # $class->SUPER::bnot($class,@_);
321 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
322 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
323 my ($self,$x,$y) = objectify(2,@_);
325 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
327 # handle +-inf and NaN
328 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
329 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
330 return +1 if $x->{sign} eq '+inf';
331 return -1 if $x->{sign} eq '-inf';
332 return -1 if $y->{sign} eq '+inf';
333 return +1 if $y->{sign} eq '-inf';
336 # check sign for speed first
337 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
338 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
341 my $xz = $x->is_zero();
342 my $yz = $y->is_zero();
343 return 0 if $xz && $yz; # 0 <=> 0
344 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
345 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
347 # adjust so that exponents are equal
348 my $lx = $x->{_m}->length() + $x->{_e};
349 my $ly = $y->{_m}->length() + $y->{_e};
350 # print "x $x y $y lx $lx ly $ly\n";
351 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
352 # print "$l $x->{sign}\n";
353 return $l if $l != 0;
355 # lengths are equal, so compare mantissa, if equal, compare exponents
356 # this assumes normalized numbers (no trailing zeros etc!)
357 my $rc = $x->{_m} <=> $y->{_m} || $x->{_e} <=> $y->{_e};
358 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
364 # Compares 2 values, ignoring their signs.
365 # Returns one of undef, <0, =0, >0. (suitable for sort)
366 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
367 my ($self,$x,$y) = objectify(2,@_);
368 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
370 # signs are ignored, so check length
371 # length(x) is length(m)+e aka length of non-fraction part
372 # the longer one is bigger
373 my $l = $x->length() - $y->length();
375 return $l if $l != 0;
376 #print "equal lengths\n";
378 # if both are equal long, make full compare
379 # first compare only the mantissa
380 # if mantissa are equal, compare fractions
382 return $x->{_m} <=> $y->{_m} || $x->{_e} <=> $y->{_e};
387 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
388 # return result as BFLOAT
389 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
391 # inf and NaN handling
392 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
395 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
397 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
399 # + and + => +, - and - => -, + and - => 0, - and + => 0
400 return $x->bzero() if $x->{sign} ne $y->{sign};
403 # +-inf + something => +inf
404 # something +-inf => +-inf
405 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
409 # speed: no add for 0+y or x+0
410 return $x if $y->is_zero(); # x+0
411 if ($x->is_zero()) # 0+y
413 # make copy, clobbering up x (modify in place!)
414 $x->{_e} = $y->{_e}->copy();
415 $x->{_m} = $y->{_m}->copy();
416 $x->{sign} = $y->{sign} || $nan;
417 return $x->round($a,$p,$r,$y);
420 # take lower of the two e's and adapt m1 to it to match m2
421 my $e = $y->{_e}; $e = Math::BigInt::bzero() if !defined $e; # if no BFLOAT
423 my $add = $y->{_m}->copy();
427 #print "\$x->{_m}: $x->{_m} ";
428 #print "\$x->{_e}: $x->{_e}\n";
429 my $e1 = $e->copy()->babs();
430 $x->{_m} *= (10 ** $e1);
431 $x->{_e} += $e; # need the sign of e
432 #$x->{_m} += $y->{_m};
433 #print "\$x->{_m}: $x->{_m} ";
434 #print "\$x->{_e}: $x->{_e}\n";
439 #print "\$x->{_m}: $x->{_m} \$y->{_m}: $y->{_m} \$e: $e ",ref($e),"\n";
441 #$x->{_m} += $y->{_m} * (10 ** $e);
442 #print "\$x->{_m}: $x->{_m}\n";
444 # else: both e are same, so leave them
445 #print "badd $x->{sign}$x->{_m} + $y->{sign}$add\n";
447 $x->{_m}->{sign} = $x->{sign};
448 $add->{sign} = $y->{sign};
452 $x->{sign} = $x->{_m}->{sign};
453 $x->{_m}->{sign} = '+';
454 #$x->bnorm(); # delete trailing zeros
455 return $x->round($a,$p,$r,$y);
460 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
461 # subtract second arg from first, modify first
462 my ($self,$x,$y) = objectify(2,@_);
464 $x->badd($y->bneg()); # badd does not leave internal zeros
465 $y->bneg(); # refix y, assumes no one reads $y in between
466 return $x; # badd() already normalized and rounded
471 # increment arg by one
472 my ($self,$x,$a,$p,$r) = objectify(1,@_);
473 $x->badd($self->_one())->round($a,$p,$r);
478 # decrement arg by one
479 my ($self,$x,$a,$p,$r) = objectify(1,@_);
480 $x->badd($self->_one('-'))->round($a,$p,$r);
485 # (BINT or num_str, BINT or num_str) return BINT
486 # does not modify arguments, but returns new object
487 # Lowest Common Multiplicator
489 my ($self,@arg) = objectify(0,@_);
490 my $x = $self->new(shift @arg);
491 while (@arg) { $x = _lcm($x,shift @arg); }
497 # (BINT or num_str, BINT or num_str) return BINT
498 # does not modify arguments, but returns new object
499 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
501 my ($self,@arg) = objectify(0,@_);
502 my $x = $self->new(shift @arg);
503 while (@arg) { $x = _gcd($x,shift @arg); }
509 # return true if arg (BINT or num_str) is zero (array '+', '0')
510 my $x = shift; $x = $class->new($x) unless ref $x;
512 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
518 # return true if arg (BINT or num_str) is +1 (array '+', '1')
519 # or -1 if signis given
520 my $x = shift; $x = $class->new($x) unless ref $x;
521 #my ($self,$x) = objectify(1,@_);
522 my $sign = $_[2] || '+';
523 return ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
528 # return true if arg (BINT or num_str) is odd or false if even
529 my $x = shift; $x = $class->new($x) unless ref $x;
530 #my ($self,$x) = objectify(1,@_);
532 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
533 return ($x->{_e}->is_zero() && $x->{_m}->is_odd());
538 # return true if arg (BINT or num_str) is even or false if odd
539 my $x = shift; $x = $class->new($x) unless ref $x;
540 #my ($self,$x) = objectify(1,@_);
542 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
543 return 1 if $x->{_m}->is_zero(); # 0e1 is even
544 return ($x->{_e}->is_zero() && $x->{_m}->is_even()); # 123.45 is never
549 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
550 # (BINT or num_str, BINT or num_str) return BINT
551 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
553 # print "mbf bmul $x->{_m}e$x->{_e} $y->{_m}e$y->{_e}\n";
554 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
557 return $x->bzero() if $x->is_zero() || $y->is_zero();
559 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
561 # result will always be +-inf:
562 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
563 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
564 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
565 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
566 return $x->binf('-');
569 # aEb * cEd = (a*c)E(b+d)
570 $x->{_m} = $x->{_m} * $y->{_m};
571 #print "m: $x->{_m}\n";
572 $x->{_e} = $x->{_e} + $y->{_e};
573 #print "e: $x->{_m}\n";
575 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
576 #print "s: $x->{sign}\n";
578 return $x->round($a,$p,$r,$y);
583 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
584 # (BFLOAT,BFLOAT) (quo,rem) or BINT (only rem)
585 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
587 # x / +-inf => 0, reminder x
588 return wantarray ? ($x->bzero(),$x->copy()) : $x->bzero()
589 if $y->{sign} =~ /^[+-]inf$/;
591 # NaN if x == NaN or y == NaN or x==y==0
592 return wantarray ? ($x->bnan(),bnan()) : $x->bnan()
593 if (($x->is_nan() || $y->is_nan()) ||
594 ($x->is_zero() && $y->is_zero()));
596 # 5 / 0 => +inf, -6 / 0 => -inf
598 ? ($x->binf($x->{sign}),$self->bnan()) : $x->binf($x->{sign})
599 if ($x->{sign} =~ /^[+-]$/ && $y->is_zero());
601 $y = $class->new($y) if ref($y) ne $class; # promote bigints
603 # print "mbf bdiv $x ",ref($x)," ",$y," ",ref($y),"\n";
604 # we need to limit the accuracy to protect against overflow
605 my ($scale) = $x->_scale_a($accuracy,$rnd_mode,$a,$r); # ignore $p
609 # simulate old behaviour
610 $scale = $div_scale+1; # one more for proper riund
611 $a = $div_scale; # and round to it
612 $fallback = 1; # to clear a/p afterwards
614 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
615 $scale = $lx if $lx > $scale;
616 $scale = $ly if $ly > $scale;
617 #print "scale $scale $lx $ly\n";
618 my $diff = $ly - $lx;
619 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
621 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
623 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
625 # check for / +-1 ( +/- 1E0)
628 return wantarray ? ($x,$self->bzero()) : $x;
631 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
632 #print "self: $self x: $x ref(x) ", ref($x)," m: $x->{_m}\n";
633 # my $scale_10 = 10 ** $scale; $x->{_m}->bmul($scale_10);
634 $x->{_m}->blsft($scale,10);
635 #print "m: $x->{_m} $y->{_m}\n";
636 $x->{_m}->bdiv( $y->{_m} ); # a/c
637 #print "m: $x->{_m}\n";
638 #print "e: $x->{_e} $y->{_e}",$scale,"\n";
639 $x->{_e}->bsub($y->{_e}); # b-d
640 #print "e: $x->{_e}\n";
641 $x->{_e}->bsub($scale); # correct for 10**scale
642 #print "after div: m: $x->{_m} e: $x->{_e}\n";
643 $x->bnorm(); # remove trailing 0's
644 #print "after div: m: $x->{_m} e: $x->{_e}\n";
645 $x->round($a,$p,$r); # then round accordingly
648 # clear a/p after round, since user did not request it
655 my $rem = $x->copy();
656 $rem->bmod($y,$a,$p,$r);
659 # clear a/p after round, since user did not request it
670 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
671 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
673 return $x->bnan() if ($x->{sign} eq $nan || $y->is_nan() || $y->is_zero());
674 return $x->bzero() if $y->is_one();
676 # XXX tels: not done yet
677 return $x->round($a,$p,$r,$y);
682 # calculate square root; this should probably
683 # use a different test to see whether the accuracy we want is...
684 my ($self,$x,$a,$p,$r) = objectify(1,@_);
686 return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
687 return $x if $x->{sign} eq '+inf'; # +inf
688 return $x if $x->is_zero() || $x == 1;
690 # we need to limit the accuracy to protect against overflow
691 my ($scale) = $x->_scale_a($accuracy,$rnd_mode,$a,$r); # ignore $p
695 # simulate old behaviour
696 $scale = $div_scale+1; # one more for proper riund
697 $a = $div_scale; # and round to it
698 $fallback = 1; # to clear a/p afterwards
700 my $lx = $x->{_m}->length();
701 $scale = $lx if $scale < $lx;
702 my $e = Math::BigFloat->new("1E-$scale"); # make test variable
703 return $x->bnan() if $e->sign() eq 'NaN';
705 # start with some reasonable guess
706 #$x *= 10 ** ($len - $org->{_e}); $x /= 2; # !?!?
709 my $gs = Math::BigFloat->new('1'. ('0' x $lx));
711 # print "first guess: $gs (x $x) scale $scale\n";
715 my $two = Math::BigFloat->new(2);
716 $x = Math::BigFloat->new($x) if ref($x) ne $class; # promote BigInts
720 return $x->bnan() if $gs->is_zero();
721 $r = $y->copy(); $r->bdiv($gs,$scale);
723 $x->bdiv($two,$scale);
724 $diff = $x->copy()->bsub($gs)->babs();
730 # clear a/p after round, since user did not request it
739 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
740 # compute power of two numbers, second arg is used as integer
741 # modifies first argument
743 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
745 return $x if $x->{sign} =~ /^[+-]inf$/;
746 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
747 return $x->bone() if $y->is_zero();
748 return $x if $x->is_one() || $y->is_one();
749 my $y1 = $y->as_number(); # make bigint
752 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
753 return $y1->is_odd() ? $x : $x->babs(1);
755 return $x if $x->is_zero() && $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
756 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
757 return $x->binf() if $x->is_zero() && $y->{sign} eq '-';
759 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
763 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
765 if ($y->{sign} eq '-')
767 # modify $x in place!
768 my $z = $x->copy(); $x->bzero()->binc();
769 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
771 return $x->round($a,$p,$r,$y);
774 ###############################################################################
779 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
780 # $n == 0 means round to integer
781 # expects and returns normalized numbers!
782 my $x = shift; $x = $class->new($x) unless ref $x;
784 return $x if $x->modify('bfround');
786 my ($scale,$mode) = $x->_scale_p($precision,$rnd_mode,@_);
787 return $x if !defined $scale; # no-op
789 # never round a 0, +-inf, NaN
790 return $x if $x->{sign} !~ /^[+-]$/ || $x->is_zero();
791 # print "MBF bfround $x to scale $scale mode $mode\n";
795 # print "bfround scale $scale e $x->{_e}\n";
796 # round right from the '.'
797 return $x if $x->{_e} >= 0; # nothing to round
798 $scale = -$scale; # positive for simplicity
799 my $len = $x->{_m}->length(); # length of mantissa
800 my $dad = -$x->{_e}; # digits after dot
801 my $zad = 0; # zeros after dot
802 $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
803 # print "scale $scale dad $dad zad $zad len $len\n";
805 # number bsstr len zad dad
807 # 0.0123 123e-4 3 1 4
810 # 1.2345 12345e-4 5 0 4
812 # do not round after/right of the $dad
813 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
815 # round to zero if rounding inside the $zad, but not for last zero like:
816 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
821 if ($scale == $zad) # for 0.006, scale -2 and trunc
827 # adjust round-point to be inside mantissa
830 $scale = $scale-$zad;
834 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
835 $scale = $dbd+$scale;
838 # print "round to $x->{_m} to $scale\n";
842 # 123 => 100 means length(123) = 3 - $scale (2) => 1
844 # calculate digits before dot
845 my $dbt = $x->{_m}->length(); $dbt += $x->{_e} if $x->{_e}->sign() eq '-';
846 if (($scale > $dbt) && ($dbt < 0))
848 # if not enough digits before dot, round to zero
851 if (($scale >= 0) && ($dbt == 0))
853 # 0.49->bfround(1): scale == 1, dbt == 0: => 0.0
854 # 0.51->bfround(0): scale == 0, dbt == 0: => 1.0
855 # 0.5->bfround(0): scale == 0, dbt == 0: => 0
856 # 0.05->bfround(0): scale == 0, dbt == 0: => 0
857 # print "$scale $dbt $x->{_m}\n";
858 $scale = -$x->{_m}->length();
862 # correct by subtracting scale
863 $scale = $dbt - $scale;
867 $scale = $x->{_m}->length() - $scale;
870 # print "using $scale for $x->{_m} with '$mode'\n";
871 # pass sign to bround for rounding modes '+inf' and '-inf'
872 $x->{_m}->{sign} = $x->{sign};
873 $x->{_m}->bround($scale,$mode);
874 $x->{_m}->{sign} = '+'; # fix sign back
880 # accuracy: preserve $N digits, and overwrite the rest with 0's
881 my $x = shift; $x = $class->new($x) unless ref $x;
882 my ($scale,$mode) = $x->_scale_a($accuracy,$rnd_mode,@_);
883 return $x if !defined $scale; # no-op
885 return $x if $x->modify('bround');
887 # print "bround $scale $mode\n";
888 # 0 => return all digits, scale < 0 makes no sense
889 return $x if ($scale <= 0);
890 # never round a 0, +-inf, NaN
891 return $x if $x->{sign} !~ /^[+-]$/ || $x->is_zero();
893 # if $e longer than $m, we have 0.0000xxxyyy style number, and must
894 # subtract the delta from scale, to simulate keeping the zeros
895 # -5 +5 => 1; -10 +5 => -4
896 my $delta = $x->{_e} + $x->{_m}->length() + 1;
897 # removed by tlr, since causes problems with fraction tests:
898 # $scale += $delta if $delta < 0;
900 # if we should keep more digits than the mantissa has, do nothing
901 return $x if $x->{_m}->length() <= $scale;
903 # pass sign to bround for '+inf' and '-inf' rounding modes
904 $x->{_m}->{sign} = $x->{sign};
905 $x->{_m}->bround($scale,$mode); # round mantissa
906 $x->{_m}->{sign} = '+'; # fix sign back
907 $x->bnorm(); # del trailing zeros gen. by bround()
912 # return integer less or equal then $x
913 my ($self,$x,$a,$p,$r) = objectify(1,@_);
915 return $x if $x->modify('bfloor');
917 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
919 # if $x has digits after dot
920 if ($x->{_e}->{sign} eq '-')
922 $x->{_m}->brsft(-$x->{_e},10);
924 $x-- if $x->{sign} eq '-';
926 return $x->round($a,$p,$r);
931 # return integer greater or equal then $x
932 my ($self,$x,$a,$p,$r) = objectify(1,@_);
934 return $x if $x->modify('bceil');
935 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
937 # if $x has digits after dot
938 if ($x->{_e}->{sign} eq '-')
940 $x->{_m}->brsft(-$x->{_e},10);
942 $x++ if $x->{sign} eq '+';
944 return $x->round($a,$p,$r);
947 ###############################################################################
951 # going trough AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
956 # make fxxx and bxxx work
958 my $name = $AUTOLOAD;
960 $name =~ s/.*:://; # split package
962 if (!method_valid($name))
966 #&{$class."::SUPER->$name"}(@_);
967 # delayed load of Carp and avoid recursion
969 Carp::croak ("Can't call $class\-\>$name, not a valid method");
972 my $bname = $name; $bname =~ s/^f/b/;
973 *{$class."\:\:$name"} = \&$bname;
979 # return a copy of the exponent
981 $self = $class->new($self) unless ref $self;
983 return bnan() if $self->is_nan();
984 return $self->{_e}->copy();
989 # return a copy of the mantissa
991 $self = $class->new($self) unless ref $self;
993 return bnan() if $self->is_nan();
994 my $m = $self->{_m}->copy(); # faster than going via bstr()
995 $m->bneg() if $self->{sign} eq '-';
1002 # return a copy of both the exponent and the mantissa
1004 $self = $class->new($self) unless ref $self;
1006 return (bnan(),bnan()) if $self->is_nan();
1007 my $m = $self->{_m}->copy(); # faster than going via bstr()
1008 $m->bneg() if $self->{sign} eq '-';
1009 return ($m,$self->{_e}->copy());
1012 ##############################################################################
1013 # private stuff (internal use only)
1017 # internal speedup, set argument to 1, or create a +/- 1
1018 my $self = shift; $self = ref($self) if ref($self);
1019 my $x = {}; bless $x, $self;
1020 $x->{_m} = Math::BigInt->new(1);
1021 $x->{_e} = Math::BigInt->new(0);
1022 $x->{sign} = shift || '+';
1029 #print "import $self\n";
1030 for ( my $i = 0; $i < @_ ; $i++ )
1032 if ( $_[$i] eq ':constant' )
1034 # this rest causes overlord er load to step in
1035 # print "overload @_\n";
1036 overload::constant float => sub { $self->new(shift); };
1037 splice @_, $i, 1; last;
1040 # any non :constant stuff is handled by our parent, Exporter
1041 # even if @_ is empty, to give it a chance
1042 #$self->SUPER::import(@_); # does not work (would call MBI)
1043 $self->export_to_level(1,$self,@_); # need this instead
1048 # adjust m and e so that m is smallest possible
1049 # round number according to accuracy and precision settings
1052 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
1054 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
1057 $x->{_m}->brsft($zeros,10); $x->{_e} += $zeros;
1059 # for something like 0Ey, set y to 1
1060 $x->{sign} = '+', $x->{_e}->bzero()->binc() if $x->{_m}->is_zero();
1061 $x->{_m}->{_f} = MB_NEVER_ROUND;
1062 $x->{_e}->{_f} = MB_NEVER_ROUND;
1063 return $x; # MBI bnorm is no-op
1066 ##############################################################################
1067 # internal calculation routines
1071 # return a bigint representation of this BigFloat number
1072 my ($self,$x) = objectify(1,@_);
1075 if ($x->{_e}->is_zero())
1077 $z = $x->{_m}->copy();
1078 $z->{sign} = $x->{sign};
1081 $z = $x->{_m}->copy();
1084 $z->brsft(-$x->{_e},10);
1088 $z->blsft($x->{_e},10);
1090 $z->{sign} = $x->{sign};
1096 my $x = shift; $x = $class->new($x) unless ref $x;
1098 my $len = $x->{_m}->length();
1099 $len += $x->{_e} if $x->{_e}->sign() eq '+';
1102 my $t = Math::BigInt::bzero();
1103 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1114 Math::BigFloat - Arbitrary size floating point math package
1121 $x = Math::BigInt->new($str); # defaults to 0
1122 $nan = Math::BigInt->bnan(); # create a NotANumber
1123 $zero = Math::BigInt->bzero();# create a "+0"
1126 $x->is_zero(); # return whether arg is zero or not
1127 $x->is_nan(); # return whether arg is NaN or not
1128 $x->is_one(); # true if arg is +1
1129 $x->is_one('-'); # true if arg is -1
1130 $x->is_odd(); # true if odd, false for even
1131 $x->is_even(); # true if even, false for odd
1132 $x->is_positive(); # true if >= 0
1133 $x->is_negative(); # true if < 0
1134 $x->is_inf(sign) # true if +inf or -inf (sign default '+')
1135 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
1136 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
1137 $x->sign(); # return the sign, either +,- or NaN
1139 # The following all modify their first argument:
1142 $x->bzero(); # set $i to 0
1143 $x->bnan(); # set $i to NaN
1145 $x->bneg(); # negation
1146 $x->babs(); # absolute value
1147 $x->bnorm(); # normalize (no-op)
1148 $x->bnot(); # two's complement (bit wise not)
1149 $x->binc(); # increment x by 1
1150 $x->bdec(); # decrement x by 1
1152 $x->badd($y); # addition (add $y to $x)
1153 $x->bsub($y); # subtraction (subtract $y from $x)
1154 $x->bmul($y); # multiplication (multiply $x by $y)
1155 $x->bdiv($y); # divide, set $i to quotient
1156 # return (quo,rem) or quo if scalar
1158 $x->bmod($y); # modulus
1159 $x->bpow($y); # power of arguments (a**b)
1160 $x->blsft($y); # left shift
1161 $x->brsft($y); # right shift
1162 # return (quo,rem) or quo if scalar
1164 $x->band($y); # bit-wise and
1165 $x->bior($y); # bit-wise inclusive or
1166 $x->bxor($y); # bit-wise exclusive or
1167 $x->bnot(); # bit-wise not (two's complement)
1169 $x->bround($N); # accuracy: preserver $N digits
1170 $x->bfround($N); # precision: round to the $Nth digit
1172 # The following do not modify their arguments:
1174 bgcd(@values); # greatest common divisor
1175 blcm(@values); # lowest common multiplicator
1177 $x->bstr(); # return string
1178 $x->bsstr(); # return string in scientific notation
1180 $x->exponent(); # return exponent as BigInt
1181 $x->mantissa(); # return mantissa as BigInt
1182 $x->parts(); # return (mantissa,exponent) as BigInt
1184 $x->length(); # number of digits (w/o sign and '.')
1185 ($l,$f) = $x->length(); # number of digits, and length of fraction
1189 All operators (inlcuding basic math operations) are overloaded if you
1190 declare your big floating point numbers as
1192 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1194 Operations with overloaded operators preserve the arguments, which is
1195 exactly what you expect.
1197 =head2 Canonical notation
1199 Input to these routines are either BigFloat objects, or strings of the
1200 following four forms:
1214 C</^[+-]\d+E[+-]?\d+$/>
1218 C</^[+-]\d*\.\d+E[+-]?\d+$/>
1222 all with optional leading and trailing zeros and/or spaces. Additonally,
1223 numbers are allowed to have an underscore between any two digits.
1225 Empty strings as well as other illegal numbers results in 'NaN'.
1227 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1228 are always stored in normalized form. On a string, it creates a BigFloat
1233 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1235 The string output will always have leading and trailing zeros stripped and drop
1236 a plus sign. C<bstr()> will give you always the form with a decimal point,
1237 while C<bsstr()> (for scientific) gives you the scientific notation.
1239 Input bstr() bsstr()
1241 ' -123 123 123' '-123123123' '-123123123E0'
1242 '00.0123' '0.0123' '123E-4'
1243 '123.45E-2' '1.2345' '12345E-4'
1244 '10E+3' '10000' '1E4'
1246 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1247 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1248 return either undef, <0, 0 or >0 and are suited for sort.
1250 Actual math is done by using BigInts to represent the mantissa and exponent.
1251 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1252 represent the result when input arguments are not numbers, as well as
1253 the result of dividing by zero.
1255 =head2 C<mantissa()>, C<exponent()> and C<parts()>
1257 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1258 as BigInts such that:
1260 $m = $x->mantissa();
1261 $e = $x->exponent();
1262 $y = $m * ( 10 ** $e );
1263 print "ok\n" if $x == $y;
1265 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
1267 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
1269 Currently the mantissa is reduced as much as possible, favouring higher
1270 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
1271 This might change in the future, so do not depend on it.
1273 =head2 Accuracy vs. Precision
1275 See also: L<Rounding|Rounding>.
1277 Math::BigFloat supports both precision and accuracy. (here should follow
1278 a short description of both).
1280 Precision: digits after the '.', laber, schwad
1281 Accuracy: Significant digits blah blah
1283 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
1284 a operation consumes all resources, each operation produces no more than
1285 C<Math::BigFloat::precision()> digits.
1287 In case the result of one operation has more precision than specified,
1288 it is rounded. The rounding mode taken is either the default mode, or the one
1289 supplied to the operation after the I<scale>:
1291 $x = Math::BigFloat->new(2);
1292 Math::BigFloat::precision(5); # 5 digits max
1293 $y = $x->copy()->bdiv(3); # will give 0.66666
1294 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1295 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
1296 Math::BigFloat::round_mode('zero');
1297 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1303 =item ffround ( +$scale )
1305 Rounds to the $scale'th place left from the '.', counting from the dot.
1306 The first digit is numbered 1.
1308 =item ffround ( -$scale )
1310 Rounds to the $scale'th place right from the '.', counting from the dot.
1314 Rounds to an integer.
1316 =item fround ( +$scale )
1318 Preserves accuracy to $scale digits from the left (aka significant digits)
1319 and pads the rest with zeros. If the number is between 1 and -1, the
1320 significant digits count from the first non-zero after the '.'
1322 =item fround ( -$scale ) and fround ( 0 )
1324 These are effetively no-ops.
1328 All rounding functions take as a second parameter a rounding mode from one of
1329 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
1331 The default rounding mode is 'even'. By using
1332 C<< Math::BigFloat::round_mode($rnd_mode); >> you can get and set the default
1333 mode for subsequent rounding. The usage of C<$Math::BigFloat::$rnd_mode> is
1334 no longer supported.
1335 The second parameter to the round functions then overrides the default
1338 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
1339 'trunc' as rounding mode to make it equivalent to:
1344 You can override this by passing the desired rounding mode as parameter to
1347 $x = Math::BigFloat->new(2.5);
1348 $y = $x->as_number('odd'); # $y = 3
1354 =head1 Autocreating constants
1356 After C<use Math::BigFloat ':constant'> all the floating point constants
1357 in the given scope are converted to C<Math::BigFloat>. This conversion
1358 happens at compile time.
1362 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
1364 prints the value of C<2E-100>. Note that without conversion of
1365 constants the expression 2E-100 will be calculated as normal floating point
1374 The following does not work yet:
1376 $m = $x->mantissa();
1377 $e = $x->exponent();
1378 $y = $m * ( 10 ** $e );
1379 print "ok\n" if $x == $y;
1383 There is no fmod() function yet.
1391 =item stringify, bstr()
1393 Both stringify and bstr() now drop the leading '+'. The old code would return
1394 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
1395 reasoning and details.
1399 The following will probably not do what you expect:
1401 print $c->bdiv(123.456),"\n";
1403 It prints both quotient and reminder since print works in list context. Also,
1404 bdiv() will modify $c, so be carefull. You probably want to use
1406 print $c / 123.456,"\n";
1407 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
1411 =item Modifying and =
1415 $x = Math::BigFloat->new(5);
1418 It will not do what you think, e.g. making a copy of $x. Instead it just makes
1419 a second reference to the B<same> object and stores it in $y. Thus anything
1420 that modifies $x will modify $y, and vice versa.
1423 print "$x, $y\n"; # prints '10, 10'
1425 If you want a true copy of $x, use:
1429 See also the documentation in L<overload> regarding C<=>.
1433 C<bpow()> now modifies the first argument, unlike the old code which left
1434 it alone and only returned the result. This is to be consistent with
1435 C<badd()> etc. The first will modify $x, the second one won't:
1437 print bpow($x,$i),"\n"; # modify $x
1438 print $x->bpow($i),"\n"; # ditto
1439 print $x ** $i,"\n"; # leave $x alone
1445 This program is free software; you may redistribute it and/or modify it under
1446 the same terms as Perl itself.
1450 Mark Biggar, overloaded interface by Ilya Zakharevich.
1451 Completely rewritten by Tels http://bloodgate.com in 2001.