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 bint 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
52 my $ten = Math::BigInt->new(10); # shortcut for speed
54 # Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
62 my %methods = map { $_ => 1 }
63 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
64 fabs fneg fint fcmp fzero fnan finc fdec
67 sub method_valid { return exists $methods{$_[0]||''}; }
70 ##############################################################################
75 # create a new BigFloat object from a string or another bigfloat object.
78 # sign => sign (+/-), or "NaN"
82 my $wanted = shift; # avoid numify call by not using || here
83 return $class->bzero() if !defined $wanted; # default to 0
84 return $wanted->copy() if ref($wanted) eq $class;
86 my $round = shift; $round = 0 if !defined $round; # no rounding as default
87 my $self = {}; bless $self, $class;
88 # shortcut for bigints and its subclasses
89 if ((ref($wanted)) && (ref($wanted) ne $class))
91 $self->{_m} = $wanted->as_number(); # get us a bigint copy
92 $self->{_e} = Math::BigInt->new(0);
94 $self->{sign} = $wanted->sign();
95 return $self->bnorm();
98 # handle '+inf', '-inf' first
99 if ($wanted =~ /^[+-]inf$/)
101 $self->{_e} = Math::BigInt->new(0);
102 $self->{_m} = Math::BigInt->new(0);
103 $self->{sign} = $wanted;
104 return $self->bnorm();
106 #print "new string '$wanted'\n";
107 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
110 die "$wanted is not a number initialized to $class" if !$NaNOK;
111 $self->{_e} = Math::BigInt->new(0);
112 $self->{_m} = Math::BigInt->new(0);
113 $self->{sign} = $nan;
117 # make integer from mantissa by adjusting exp, then convert to bigint
118 $self->{_e} = Math::BigInt->new("$$es$$ev"); # exponent
119 $self->{_m} = Math::BigInt->new("$$mis$$miv$$mfv"); # create mantissa
120 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
121 $self->{_e} -= CORE::length($$mfv);
122 $self->{sign} = $self->{_m}->sign(); $self->{_m}->babs();
124 #print "$wanted => $self->{sign} $self->{value}\n";
125 $self->bnorm(); # first normalize
126 # if any of the globals is set, round to them and thus store them insid $self
127 $self->round($accuracy,$precision,$rnd_mode)
128 if defined $accuracy || defined $precision;
132 # some shortcuts for easier life
135 # exportable version of new
136 return $class->new(@_);
141 # exportable version of new
142 return $class->new(@_,0)->bround(0,'trunc');
147 # create a bigfloat 'NaN', if given a BigFloat, set it to 'NaN'
149 $self = $class if !defined $self;
152 my $c = $self; $self = {}; bless $self, $c;
154 $self->{_e} = new Math::BigInt 0;
155 $self->{_m} = new Math::BigInt 0;
156 $self->{sign} = $nan;
162 # create a bigfloat '+-inf', if given a BigFloat, set it to '+-inf'
164 my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
166 $self = $class if !defined $self;
169 my $c = $self; $self = {}; bless $self, $c;
171 $self->{_e} = new Math::BigInt 0;
172 $self->{_m} = new Math::BigInt 0;
173 $self->{sign} = $sign.'inf';
179 # create a bigfloat '+0', if given a BigFloat, set it to 0
181 $self = $class if !defined $self;
184 my $c = $self; $self = {}; bless $self, $c;
186 $self->{_m} = new Math::BigInt 0;
187 $self->{_e} = new Math::BigInt 1;
192 ##############################################################################
193 # string conversation
197 # (ref to BFLOAT or num_str ) return num_str
198 # Convert number from internal format to (non-scientific) string format.
199 # internal format is always normalized (no leading zeros, "-0" => "+0")
200 my ($self,$x) = objectify(1,@_);
202 #return "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
203 #return "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
204 return $x->{sign} if $x->{sign} !~ /^[+-]$/;
205 return '0' if $x->is_zero();
207 my $es = $x->{_m}->bstr();
208 if ($x->{_e}->is_zero())
210 $es = $x->{sign}.$es if $x->{sign} eq '-';
214 if ($x->{_e}->sign() eq '-')
216 if ($x->{_e} <= -CORE::length($es))
218 # print "style: 0.xxxx\n";
219 my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
220 $es = '0.'. ('0' x $r) . $es;
224 # print "insert '.' at $x->{_e} in '$es'\n";
225 substr($es,$x->{_e},0) = '.';
231 $es .= '0' x $x->{_e};
233 $es = $x->{sign}.$es if $x->{sign} eq '-';
239 # (ref to BFLOAT or num_str ) return num_str
240 # Convert number from internal format to scientific string format.
241 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
242 my ($self,$x) = objectify(1,@_);
244 return "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
245 return "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
246 return $x->{sign} if $x->{sign} !~ /^[+-]$/;
247 my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
249 return $x->{_m}->bstr().$sep.$x->{_e}->bstr();
254 # Make a number from a BigFloat object
255 # simple return string and let Perl's atoi() handle the rest
256 my ($self,$x) = objectify(1,@_);
260 ##############################################################################
261 # public stuff (usually prefixed with "b")
263 # really? Just for exporting them is not what I had in mind
266 # $class->SUPER::babs($class,@_);
270 # $class->SUPER::bneg($class,@_);
274 # $class->SUPER::bnot($class,@_);
279 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
280 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
281 my ($self,$x,$y) = objectify(2,@_);
283 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
285 # handle +-inf and NaN
286 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
287 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
288 return +1 if $x->{sign} eq '+inf';
289 return -1 if $x->{sign} eq '-inf';
290 return -1 if $y->{sign} eq '+inf';
291 return +1 if $y->{sign} eq '-inf';
294 # check sign for speed first
295 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-';
296 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
298 return 0 if $x->is_zero() && $y->is_zero(); # 0 <=> 0
299 return -1 if $x->is_zero() && $y->{sign} eq '+'; # 0 <=> +y
300 return 1 if $y->is_zero() && $x->{sign} eq '+'; # +x <=> 0
302 # adjust so that exponents are equal
303 my $lx = $x->{_m}->length() + $x->{_e};
304 my $ly = $y->{_m}->length() + $y->{_e};
305 # print "x $x y $y lx $lx ly $ly\n";
306 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
307 # print "$l $x->{sign}\n";
308 return $l if $l != 0;
310 # lengths are equal, so compare mantissa, if equal, compare exponents
311 # this assumes normalized numbers (no trailing zeros etc!)
312 my $rc = $x->{_m} <=> $y->{_m} || $x->{_e} <=> $y->{_e};
313 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
319 # Compares 2 values, ignoring their signs.
320 # 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,@_);
323 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
325 # signs are ignored, so check length
326 # length(x) is length(m)+e aka length of non-fraction part
327 # the longer one is bigger
328 my $l = $x->length() - $y->length();
330 return $l if $l != 0;
331 #print "equal lengths\n";
333 # if both are equal long, make full compare
334 # first compare only the mantissa
335 # if mantissa are equal, compare fractions
337 return $x->{_m} <=> $y->{_m} || $x->{_e} <=> $y->{_e};
342 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
343 # return result as BFLOAT
344 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
346 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
348 # speed: no add for 0+y or x+0
349 return $x if $y->is_zero(); # x+0
350 if ($x->is_zero()) # 0+y
352 # make copy, clobbering up x (modify in place!)
353 $x->{_e} = $y->{_e}->copy();
354 $x->{_m} = $y->{_m}->copy();
355 $x->{sign} = $y->{sign} || $nan;
356 return $x->round($a,$p,$r,$y);
359 # take lower of the two e's and adapt m1 to it to match m2
360 my $e = $y->{_e}; $e = Math::BigInt::bzero() if !defined $e; # if no BFLOAT
362 my $add = $y->{_m}->copy();
366 #print "\$x->{_m}: $x->{_m} ";
367 #print "\$x->{_e}: $x->{_e}\n";
368 my $e1 = $e->copy()->babs();
369 $x->{_m} *= (10 ** $e1);
370 $x->{_e} += $e; # need the sign of e
371 #$x->{_m} += $y->{_m};
372 #print "\$x->{_m}: $x->{_m} ";
373 #print "\$x->{_e}: $x->{_e}\n";
378 #print "\$x->{_m}: $x->{_m} \$y->{_m}: $y->{_m} \$e: $e ",ref($e),"\n";
380 #$x->{_m} += $y->{_m} * (10 ** $e);
381 #print "\$x->{_m}: $x->{_m}\n";
383 # else: both e are same, so leave them
384 #print "badd $x->{sign}$x->{_m} + $y->{sign}$add\n";
386 $x->{_m}->{sign} = $x->{sign};
387 $add->{sign} = $y->{sign};
391 $x->{sign} = $x->{_m}->{sign};
392 $x->{_m}->{sign} = '+';
393 #$x->bnorm(); # delete trailing zeros
394 return $x->round($a,$p,$r,$y);
399 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
400 # subtract second arg from first, modify first
401 my ($self,$x,$y) = objectify(2,@_);
403 $x->badd($y->bneg()); # badd does not leave internal zeros
404 $y->bneg(); # refix y, assumes no one reads $y in between
405 return $x; # badd() already normalized and rounded
410 # increment arg by one
411 my ($self,$x,$a,$p,$r) = objectify(1,@_);
412 $x->badd($self->_one())->round($a,$p,$r);
417 # decrement arg by one
418 my ($self,$x,$a,$p,$r) = objectify(1,@_);
419 $x->badd($self->_one('-'))->round($a,$p,$r);
424 # (BINT or num_str, BINT or num_str) return BINT
425 # does not modify arguments, but returns new object
426 # Lowest Common Multiplicator
428 my ($self,@arg) = objectify(0,@_);
429 my $x = $self->new(shift @arg);
430 while (@arg) { $x = _lcm($x,shift @arg); }
436 # (BINT or num_str, BINT or num_str) return BINT
437 # does not modify arguments, but returns new object
438 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
440 my ($self,@arg) = objectify(0,@_);
441 my $x = $self->new(shift @arg);
442 while (@arg) { $x = _gcd($x,shift @arg); }
448 # return true if arg (BINT or num_str) is zero (array '+', '0')
449 my $x = shift; $x = $class->new($x) unless ref $x;
450 #my ($self,$x) = objectify(1,@_);
451 return ($x->{sign} ne $nan && $x->{_m}->is_zero());
456 # return true if arg (BINT or num_str) is +1 (array '+', '1')
457 # or -1 if signis given
458 my $x = shift; $x = $class->new($x) unless ref $x;
459 #my ($self,$x) = objectify(1,@_);
460 my $sign = $_[2] || '+';
461 return ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
466 # return true if arg (BINT or num_str) is odd or false if even
467 my $x = shift; $x = $class->new($x) unless ref $x;
468 #my ($self,$x) = objectify(1,@_);
470 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
471 return ($x->{_e}->is_zero() && $x->{_m}->is_odd());
476 # return true if arg (BINT or num_str) is even or false if odd
477 my $x = shift; $x = $class->new($x) unless ref $x;
478 #my ($self,$x) = objectify(1,@_);
480 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
481 return 1 if $x->{_m}->is_zero(); # 0e1 is even
482 return ($x->{_e}->is_zero() && $x->{_m}->is_even()); # 123.45 is never
487 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
488 # (BINT or num_str, BINT or num_str) return BINT
489 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
491 # print "mbf bmul $x->{_m}e$x->{_e} $y->{_m}e$y->{_e}\n";
492 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
494 # aEb * cEd = (a*c)E(b+d)
495 $x->{_m} = $x->{_m} * $y->{_m};
496 #print "m: $x->{_m}\n";
497 $x->{_e} = $x->{_e} + $y->{_e};
498 #print "e: $x->{_m}\n";
500 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
501 #print "s: $x->{sign}\n";
503 return $x->round($a,$p,$r,$y);
508 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
509 # (BFLOAT,BFLOAT) (quo,rem) or BINT (only rem)
510 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
512 return wantarray ? ($x->bnan(),bnan()) : $x->bnan()
513 if ($x->{sign} eq $nan || $y->is_nan() || $y->is_zero());
515 $y = $class->new($y) if ref($y) ne $class; # promote bigints
517 # print "mbf bdiv $x ",ref($x)," ",$y," ",ref($y),"\n";
518 # we need to limit the accuracy to protect against overflow
519 my ($scale) = $x->_scale_a($accuracy,$rnd_mode,$a,$r); # ignore $p
522 # simulate old behaviour
523 $scale = $div_scale+1; # one more for proper riund
524 $a = $div_scale; # and round to it
526 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
527 $scale = $lx if $lx > $scale;
528 $scale = $ly if $ly > $scale;
529 #print "scale $scale $lx $ly\n";
530 my $diff = $ly - $lx;
531 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
533 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
535 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
537 # check for / +-1 ( +/- 1E0)
540 return wantarray ? ($x,$self->bzero()) : $x;
543 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
544 #print "self: $self x: $x ref(x) ", ref($x)," m: $x->{_m}\n";
545 # my $scale_10 = 10 ** $scale; $x->{_m}->bmul($scale_10);
546 $x->{_m}->blsft($scale,10);
547 #print "m: $x->{_m} $y->{_m}\n";
548 $x->{_m}->bdiv( $y->{_m} ); # a/c
549 #print "m: $x->{_m}\n";
550 #print "e: $x->{_e} $y->{_e}",$scale,"\n";
551 $x->{_e}->bsub($y->{_e}); # b-d
552 #print "e: $x->{_e}\n";
553 $x->{_e}->bsub($scale); # correct for 10**scale
554 #print "after div: m: $x->{_m} e: $x->{_e}\n";
555 $x->bnorm(); # remove trailing 0's
556 #print "after div: m: $x->{_m} e: $x->{_e}\n";
557 $x->round($a,$p,$r); # then round accordingly
561 my $rem = $x->copy();
562 $rem->bmod($y,$a,$p,$r);
570 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
571 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
573 return $x->bnan() if ($x->{sign} eq $nan || $y->is_nan() || $y->is_zero());
574 return $x->bzero() if $y->is_one();
576 # XXX tels: not done yet
577 return $x->round($a,$p,$r,$y);
582 # calculate square root; this should probably
583 # use a different test to see whether the accuracy we want is...
584 my ($self,$x,$a,$p,$r) = objectify(1,@_);
586 return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
587 return $x if $x->{sign} eq '+inf'; # +inf
588 return $x if $x->is_zero() || $x == 1;
590 # we need to limit the accuracy to protect against overflow
591 my ($scale) = $x->_scale_a($accuracy,$rnd_mode,$a,$r); # ignore $p
594 # simulate old behaviour
595 $scale = $div_scale+1; # one more for proper riund
596 $a = $div_scale; # and round to it
598 my $lx = $x->{_m}->length();
599 $scale = $lx if $scale < $lx;
600 my $e = Math::BigFloat->new("1E-$scale"); # make test variable
601 return $x->bnan() if $e->sign() eq 'NaN';
603 # start with some reasonable guess
604 #$x *= 10 ** ($len - $org->{_e}); $x /= 2; # !?!?
606 my $gs = Math::BigFloat->new('1'. ('0' x $lx));
608 # print "first guess: $gs (x $x) scale $scale\n";
612 my $two = Math::BigFloat->new(2);
613 $x = Math::BigFloat->new($x) if ref($x) ne $class; # promote BigInts
617 return $x->bnan() if $gs->is_zero();
618 $r = $y->copy(); $r->bdiv($gs,$scale);
620 $x->bdiv($two,$scale);
621 $diff = $x->copy()->bsub($gs)->babs();
629 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
630 # compute power of two numbers, second arg is used as integer
631 # modifies first argument
633 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
635 return $x if $x->{sign} =~ /^[+-]inf$/;
636 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
637 return $x->bzero()->binc() if $y->is_zero();
638 return $x if $x->is_one() || $y->is_one();
639 my $y1 = $y->as_number(); # make bigint
642 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
643 return $y1->is_odd() ? $x : $x->babs(1);
645 return $x if $x->is_zero() && $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
646 # 0 ** -y => 1 / (0 ** y) => / 0!
647 return $x->bnan() if $x->is_zero() && $y->{sign} eq '-';
649 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
653 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
655 if ($y->{sign} eq '-')
657 # modify $x in place!
658 my $z = $x->copy(); $x->bzero()->binc();
659 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
661 return $x->round($a,$p,$r,$y);
664 ###############################################################################
669 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
670 # $n == 0 means round to integer
671 # expects and returns normalized numbers!
672 my $x = shift; $x = $class->new($x) unless ref $x;
674 return $x if $x->modify('bfround');
676 my ($scale,$mode) = $x->_scale_p($precision,$rnd_mode,@_);
677 return $x if !defined $scale; # no-op
679 # print "MBF bfround $x to scale $scale mode $mode\n";
680 return $x if $x->is_nan() or $x->is_zero();
684 # print "bfround scale $scale e $x->{_e}\n";
685 # round right from the '.'
686 return $x if $x->{_e} >= 0; # nothing to round
687 $scale = -$scale; # positive for simplicity
688 my $len = $x->{_m}->length(); # length of mantissa
689 my $dad = -$x->{_e}; # digits after dot
690 my $zad = 0; # zeros after dot
691 $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
692 # print "scale $scale dad $dad zad $zad len $len\n";
694 # number bsstr len zad dad
696 # 0.0123 123e-4 3 1 4
699 # 1.2345 12345e-4 5 0 4
701 # do not round after/right of the $dad
702 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
704 # round to zero if rounding inside the $zad, but not for last zero like:
705 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
708 $x->{_m} = Math::BigInt->new(0);
709 $x->{_e} = Math::BigInt->new(1);
713 if ($scale == $zad) # for 0.006, scale -2 and trunc
719 # adjust round-point to be inside mantissa
722 $scale = $scale-$zad;
726 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
727 $scale = $dbd+$scale;
730 # print "round to $x->{_m} to $scale\n";
734 # 123 => 100 means length(123) = 3 - $scale (2) => 1
736 # calculate digits before dot
737 my $dbt = $x->{_m}->length(); $dbt += $x->{_e} if $x->{_e}->sign() eq '-';
738 if (($scale > $dbt) && ($dbt < 0))
740 # if not enough digits before dot, round to zero
741 $x->{_m} = Math::BigInt->new(0);
742 $x->{_e} = Math::BigInt->new(1);
746 if (($scale >= 0) && ($dbt == 0))
748 # 0.49->bfround(1): scale == 1, dbt == 0: => 0.0
749 # 0.51->bfround(0): scale == 0, dbt == 0: => 1.0
750 # 0.5->bfround(0): scale == 0, dbt == 0: => 0
751 # 0.05->bfround(0): scale == 0, dbt == 0: => 0
752 # print "$scale $dbt $x->{_m}\n";
753 $scale = -$x->{_m}->length();
757 # correct by subtracting scale
758 $scale = $dbt - $scale;
762 $scale = $x->{_m}->length() - $scale;
765 #print "using $scale for $x->{_m} with '$mode'\n";
766 # pass sign to bround for '+inf' and '-inf' rounding modes
767 $x->{_m}->{sign} = $x->{sign};
768 $x->{_m}->bround($scale,$mode);
769 $x->{_m}->{sign} = '+'; # fix sign back
775 # accuracy: preserve $N digits, and overwrite the rest with 0's
776 my $x = shift; $x = $class->new($x) unless ref $x;
777 my ($scale,$mode) = $x->_scale_a($accuracy,$rnd_mode,@_);
778 return $x if !defined $scale; # no-op
780 return $x if $x->modify('bround');
782 # print "bround $scale $mode\n";
783 # 0 => return all digits, scale < 0 makes no sense
784 return $x if ($scale <= 0);
785 return $x if $x->is_nan() or $x->is_zero(); # never round a 0
787 # if $e longer than $m, we have 0.0000xxxyyy style number, and must
788 # subtract the delta from scale, to simulate keeping the zeros
789 # -5 +5 => 1; -10 +5 => -4
790 my $delta = $x->{_e} + $x->{_m}->length() + 1;
791 # removed by tlr, since causes problems with fraction tests:
792 # $scale += $delta if $delta < 0;
794 # if we should keep more digits than the mantissa has, do nothing
795 return $x if $x->{_m}->length() <= $scale;
797 # pass sign to bround for '+inf' and '-inf' rounding modes
798 $x->{_m}->{sign} = $x->{sign};
799 $x->{_m}->bround($scale,$mode); # round mantissa
800 $x->{_m}->{sign} = '+'; # fix sign back
801 return $x->bnorm(); # del trailing zeros gen. by bround()
806 # return integer less or equal then $x
807 my ($self,$x,$a,$p,$r) = objectify(1,@_);
809 return $x if $x->modify('bfloor');
811 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
813 # if $x has digits after dot
814 if ($x->{_e}->{sign} eq '-')
816 $x->{_m}->brsft(-$x->{_e},10);
818 $x-- if $x->{sign} eq '-';
820 return $x->round($a,$p,$r);
825 # return integer greater or equal then $x
826 my ($self,$x,$a,$p,$r) = objectify(1,@_);
828 return $x if $x->modify('bceil');
829 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
831 # if $x has digits after dot
832 if ($x->{_e}->{sign} eq '-')
834 $x->{_m}->brsft(-$x->{_e},10);
836 $x++ if $x->{sign} eq '+';
838 return $x->round($a,$p,$r);
841 ###############################################################################
845 # going trough AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
850 # make fxxx and bxxx work
852 my $name = $AUTOLOAD;
854 $name =~ s/.*:://; # split package
856 if (!method_valid($name))
860 #&{$class."::SUPER->$name"}(@_);
861 # delayed load of Carp and avoid recursion
863 Carp::croak ("Can't call $class\-\>$name, not a valid method");
866 my $bname = $name; $bname =~ s/^f/b/;
867 *{$class."\:\:$name"} = \&$bname;
873 # return a copy of the exponent
875 $self = $class->new($self) unless ref $self;
877 return bnan() if $self->is_nan();
878 return $self->{_e}->copy();
883 # return a copy of the mantissa
885 $self = $class->new($self) unless ref $self;
887 return bnan() if $self->is_nan();
888 my $m = $self->{_m}->copy(); # faster than going via bstr()
889 $m->bneg() if $self->{sign} eq '-';
896 # return a copy of both the exponent and the mantissa
898 $self = $class->new($self) unless ref $self;
900 return (bnan(),bnan()) if $self->is_nan();
901 my $m = $self->{_m}->copy(); # faster than going via bstr()
902 $m->bneg() if $self->{sign} eq '-';
903 return ($m,$self->{_e}->copy());
906 ##############################################################################
907 # private stuff (internal use only)
911 # internal speedup, set argument to 1, or create a +/- 1
912 my $self = shift; $self = ref($self) if ref($self);
913 my $x = {}; bless $x, $self;
914 $x->{_m} = Math::BigInt->new(1);
915 $x->{_e} = Math::BigInt->new(0);
916 $x->{sign} = shift || '+';
923 #print "import $self\n";
924 for ( my $i = 0; $i < @_ ; $i++ )
926 if ( $_[$i] eq ':constant' )
928 # this rest causes overlord er load to step in
929 # print "overload @_\n";
930 overload::constant float => sub { $self->new(shift); };
931 splice @_, $i, 1; last;
934 # any non :constant stuff is handled by our parent, Exporter
935 # even if @_ is empty, to give it a chance
936 #$self->SUPER::import(@_); # does not work (would call MBI)
937 $self->export_to_level(1,$self,@_); # need this instead
942 # adjust m and e so that m is smallest possible
943 # round number according to accuracy and precision settings
946 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
948 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
951 $x->{_m}->brsft($zeros,10); $x->{_e} += $zeros;
953 # for something like 0Ey, set y to 1
954 $x->{_e}->bzero()->binc() if $x->{_m}->is_zero();
955 $x->{_m}->{_f} = MB_NEVER_ROUND;
956 $x->{_e}->{_f} = MB_NEVER_ROUND;
957 return $x; # MBI bnorm is no-op
960 ##############################################################################
961 # internal calculation routines
965 # return a bigint representation of this BigFloat number
966 my ($self,$x) = objectify(1,@_);
969 if ($x->{_e}->is_zero())
971 $z = $x->{_m}->copy();
972 $z->{sign} = $x->{sign};
975 $z = $x->{_m}->copy();
978 $z->brsft(-$x->{_e},10);
982 $z->blsft($x->{_e},10);
984 $z->{sign} = $x->{sign};
990 my $x = shift; $x = $class->new($x) unless ref $x;
992 my $len = $x->{_m}->length();
993 $len += $x->{_e} if $x->{_e}->sign() eq '+';
996 my $t = Math::BigInt::bzero();
997 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1008 Math::BigFloat - Arbitrary size floating point math package
1015 $x = Math::BigInt->new($str); # defaults to 0
1016 $nan = Math::BigInt->bnan(); # create a NotANumber
1017 $zero = Math::BigInt->bzero();# create a "+0"
1020 $x->is_zero(); # return whether arg is zero or not
1021 $x->is_nan(); # return whether arg is NaN or not
1022 $x->is_one(); # true if arg is +1
1023 $x->is_one('-'); # true if arg is -1
1024 $x->is_odd(); # true if odd, false for even
1025 $x->is_even(); # true if even, false for odd
1026 $x->is_positive(); # true if >= 0
1027 $x->is_negative(); # true if < 0
1028 $x->is_inf(sign) # true if +inf or -inf (sign default '+')
1029 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
1030 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
1031 $x->sign(); # return the sign, either +,- or NaN
1033 # The following all modify their first argument:
1036 $x->bzero(); # set $i to 0
1037 $x->bnan(); # set $i to NaN
1039 $x->bneg(); # negation
1040 $x->babs(); # absolute value
1041 $x->bnorm(); # normalize (no-op)
1042 $x->bnot(); # two's complement (bit wise not)
1043 $x->binc(); # increment x by 1
1044 $x->bdec(); # decrement x by 1
1046 $x->badd($y); # addition (add $y to $x)
1047 $x->bsub($y); # subtraction (subtract $y from $x)
1048 $x->bmul($y); # multiplication (multiply $x by $y)
1049 $x->bdiv($y); # divide, set $i to quotient
1050 # return (quo,rem) or quo if scalar
1052 $x->bmod($y); # modulus
1053 $x->bpow($y); # power of arguments (a**b)
1054 $x->blsft($y); # left shift
1055 $x->brsft($y); # right shift
1056 # return (quo,rem) or quo if scalar
1058 $x->band($y); # bit-wise and
1059 $x->bior($y); # bit-wise inclusive or
1060 $x->bxor($y); # bit-wise exclusive or
1061 $x->bnot(); # bit-wise not (two's complement)
1063 $x->bround($N); # accuracy: preserver $N digits
1064 $x->bfround($N); # precision: round to the $Nth digit
1066 # The following do not modify their arguments:
1068 bgcd(@values); # greatest common divisor
1069 blcm(@values); # lowest common multiplicator
1071 $x->bstr(); # return string
1072 $x->bsstr(); # return string in scientific notation
1074 $x->exponent(); # return exponent as BigInt
1075 $x->mantissa(); # return mantissa as BigInt
1076 $x->parts(); # return (mantissa,exponent) as BigInt
1078 $x->length(); # number of digits (w/o sign and '.')
1079 ($l,$f) = $x->length(); # number of digits, and length of fraction
1083 All operators (inlcuding basic math operations) are overloaded if you
1084 declare your big floating point numbers as
1086 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1088 Operations with overloaded operators preserve the arguments, which is
1089 exactly what you expect.
1091 =head2 Canonical notation
1093 Input to these routines are either BigFloat objects, or strings of the
1094 following four forms:
1108 C</^[+-]\d+E[+-]?\d+$/>
1112 C</^[+-]\d*\.\d+E[+-]?\d+$/>
1116 all with optional leading and trailing zeros and/or spaces. Additonally,
1117 numbers are allowed to have an underscore between any two digits.
1119 Empty strings as well as other illegal numbers results in 'NaN'.
1121 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1122 are always stored in normalized form. On a string, it creates a BigFloat
1127 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1129 The string output will always have leading and trailing zeros stripped and drop
1130 a plus sign. C<bstr()> will give you always the form with a decimal point,
1131 while C<bsstr()> (for scientific) gives you the scientific notation.
1133 Input bstr() bsstr()
1135 ' -123 123 123' '-123123123' '-123123123E0'
1136 '00.0123' '0.0123' '123E-4'
1137 '123.45E-2' '1.2345' '12345E-4'
1138 '10E+3' '10000' '1E4'
1140 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1141 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1142 return either undef, <0, 0 or >0 and are suited for sort.
1144 Actual math is done by using BigInts to represent the mantissa and exponent.
1145 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1146 represent the result when input arguments are not numbers, as well as
1147 the result of dividing by zero.
1149 =head2 C<mantissa()>, C<exponent()> and C<parts()>
1151 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1152 as BigInts such that:
1154 $m = $x->mantissa();
1155 $e = $x->exponent();
1156 $y = $m * ( 10 ** $e );
1157 print "ok\n" if $x == $y;
1159 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
1161 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
1163 Currently the mantissa is reduced as much as possible, favouring higher
1164 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
1165 This might change in the future, so do not depend on it.
1167 =head2 Accuracy vs. Precision
1169 See also: L<Rounding|Rounding>.
1171 Math::BigFloat supports both precision and accuracy. (here should follow
1172 a short description of both).
1174 Precision: digits after the '.', laber, schwad
1175 Accuracy: Significant digits blah blah
1177 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
1178 a operation consumes all resources, each operation produces no more than
1179 C<Math::BigFloat::precision()> digits.
1181 In case the result of one operation has more precision than specified,
1182 it is rounded. The rounding mode taken is either the default mode, or the one
1183 supplied to the operation after the I<scale>:
1185 $x = Math::BigFloat->new(2);
1186 Math::BigFloat::precision(5); # 5 digits max
1187 $y = $x->copy()->bdiv(3); # will give 0.66666
1188 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1189 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
1190 Math::BigFloat::round_mode('zero');
1191 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1197 =item ffround ( +$scale )
1199 Rounds to the $scale'th place left from the '.', counting from the dot.
1200 The first digit is numbered 1.
1202 =item ffround ( -$scale )
1204 Rounds to the $scale'th place right from the '.', counting from the dot.
1208 Rounds to an integer.
1210 =item fround ( +$scale )
1212 Preserves accuracy to $scale digits from the left (aka significant digits)
1213 and pads the rest with zeros. If the number is between 1 and -1, the
1214 significant digits count from the first non-zero after the '.'
1216 =item fround ( -$scale ) and fround ( 0 )
1218 These are effetively no-ops.
1222 All rounding functions take as a second parameter a rounding mode from one of
1223 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
1225 The default rounding mode is 'even'. By using
1226 C<< Math::BigFloat::round_mode($rnd_mode); >> you can get and set the default
1227 mode for subsequent rounding. The usage of C<$Math::BigFloat::$rnd_mode> is
1228 no longer supported.
1229 The second parameter to the round functions then overrides the default
1232 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
1233 'trunc' as rounding mode to make it equivalent to:
1238 You can override this by passing the desired rounding mode as parameter to
1241 $x = Math::BigFloat->new(2.5);
1242 $y = $x->as_number('odd'); # $y = 3
1246 use Math::BigFloat qw(bstr bint);
1248 $x = bstr("1234") # string "1234"
1249 $x = "$x"; # same as bstr()
1250 $x = bneg("1234") # BigFloat "-1234"
1251 $x = Math::BigFloat->bneg("1234"); # BigFloat "1234"
1252 $x = Math::BigFloat->babs("-12345"); # BigFloat "12345"
1253 $x = Math::BigFloat->bnorm("-0 00"); # BigFloat "0"
1254 $x = bint(1) + bint(2); # BigFloat "3"
1255 $x = bint(1) + "2"; # ditto (auto-BigFloatify of "2")
1256 $x = bint(1); # BigFloat "1"
1257 $x = $x + 5 / 2; # BigFloat "3"
1258 $x = $x ** 3; # BigFloat "27"
1259 $x *= 2; # BigFloat "54"
1260 $x = new Math::BigFloat; # BigFloat "0"
1261 $x--; # BigFloat "-1"
1263 =head1 Autocreating constants
1265 After C<use Math::BigFloat ':constant'> all the floating point constants
1266 in the given scope are converted to C<Math::BigFloat>. This conversion
1267 happens at compile time.
1271 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
1273 prints the value of C<2E-100>. Note that without conversion of
1274 constants the expression 2E-100 will be calculated as normal floating point
1279 Greatly enhanced ;o)
1288 The following does not work yet:
1290 $m = $x->mantissa();
1291 $e = $x->exponent();
1292 $y = $m * ( 10 ** $e );
1293 print "ok\n" if $x == $y;
1297 There is no fmod() function yet.
1305 =item stringify, bstr()
1307 Both stringify and bstr() now drop the leading '+'. The old code would return
1308 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
1309 reasoning and details.
1313 The following will probably not do what you expect:
1315 print $c->bdiv(123.456),"\n";
1317 It prints both quotient and reminder since print works in list context. Also,
1318 bdiv() will modify $c, so be carefull. You probably want to use
1320 print $c / 123.456,"\n";
1321 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
1325 =item Modifying and =
1329 $x = Math::BigFloat->new(5);
1332 It will not do what you think, e.g. making a copy of $x. Instead it just makes
1333 a second reference to the B<same> object and stores it in $y. Thus anything
1334 that modifies $x will modify $y, and vice versa.
1337 print "$x, $y\n"; # prints '10, 10'
1339 If you want a true copy of $x, use:
1343 See also the documentation in L<overload> regarding C<=>.
1347 C<bpow()> now modifies the first argument, unlike the old code which left
1348 it alone and only returned the result. This is to be consistent with
1349 C<badd()> etc. The first will modify $x, the second one won't:
1351 print bpow($x,$i),"\n"; # modify $x
1352 print $x->bpow($i),"\n"; # ditto
1353 print $x ** $i,"\n"; # leave $x alone
1359 This program is free software; you may redistribute it and/or modify it under
1360 the same terms as Perl itself.
1364 Mark Biggar, overloaded interface by Ilya Zakharevich.
1365 Completely rewritten by Tels http://bloodgate.com in 2001.