3 # mark.biggar@TrustedSysLabs.com
5 # The following hash values are internally used:
6 # _e: exponent (BigInt)
7 # _m: mantissa (absolute BigInt)
8 # sign: +,-,"NaN" if not a number
11 # _cow: Copy-On-Write (NRY)
13 package Math::BigFloat;
18 use Math::BigInt qw/trace objectify/;
19 @ISA = qw( Exporter Math::BigInt);
20 # can not export bneg/babs since the are only in MBI
23 badd bmul bdiv bmod bnorm bsub
24 bgcd blcm bround bfround
25 bpow bnan bzero bfloor bceil
26 bacmp bstr binc bdec bint binf
27 is_odd is_even is_nan is_inf
33 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $rnd_mode/;
34 my $class = "Math::BigFloat";
39 $class->bcmp($_[1],$_[0]) :
40 $class->bcmp($_[0],$_[1])},
41 'int' => sub { $_[0]->copy()->bround(0,'trunc'); },
46 # set to 1 for tracing
48 # constant for easier life
50 my $ten = Math::BigInt->new(10); # shortcut for speed
52 # Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
60 my %methods = map { $_ => 1 }
61 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
62 fabs fneg fint fcmp fzero fnan finc fdec
65 sub method_valid { return exists $methods{$_[0]||''}; }
68 ##############################################################################
73 # create a new BigFloat object from a string or another bigfloat object.
76 # sign => sign (+/-), or "NaN"
81 my $wanted = shift; # avoid numify call by not using || here
82 return $class->bzero() if !defined $wanted; # default to 0
83 return $wanted->copy() if ref($wanted) eq $class;
85 my $round = shift; $round = 0 if !defined $round; # no rounding as default
86 my $self = {}; bless $self, $class;
88 if (ref($wanted) eq 'Math::BigInt')
90 $self->{_m} = $wanted;
91 $self->{_e} = Math::BigInt->new(0);
93 $self->{sign} = $wanted->sign();
97 # handle '+inf', '-inf' first
98 if ($wanted =~ /^[+-]inf$/)
100 $self->{_e} = Math::BigInt->new(0);
101 $self->{_m} = Math::BigInt->new(0);
102 $self->{sign} = $wanted;
105 #print "new string '$wanted'\n";
106 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
109 die "$wanted is not a number initialized to $class" if !$NaNOK;
110 $self->{_e} = Math::BigInt->new(0);
111 $self->{_m} = Math::BigInt->new(0);
112 $self->{sign} = $nan;
116 # make integer from mantissa by adjusting exp, then convert to bigint
117 $self->{_e} = Math::BigInt->new("$$es$$ev"); # exponent
118 $self->{_m} = Math::BigInt->new("$$mis$$miv$$mfv"); # create mantissa
119 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
120 $self->{_e} -= CORE::length($$mfv);
121 $self->{sign} = $self->{_m}->sign(); $self->{_m}->babs();
123 #print "$wanted => $self->{sign} $self->{value}->[0]\n";
124 $self->bnorm(); # first normalize
125 # if any of the globals is set, round to them and thus store them insid $self
126 $self->round($accuracy,$precision,$rnd_mode)
127 if defined $accuracy || defined $precision;
131 # some shortcuts for easier life
134 # exportable version of new
136 return $class->new(@_);
141 # exportable version of new
143 return $class->new(@_,0)->bround(0,'trunc');
148 # create a bigfloat 'NaN', if given a BigFloat, set it to 'NaN'
150 $self = $class if !defined $self;
153 my $c = $self; $self = {}; bless $self, $c;
155 $self->{_e} = new Math::BigInt 0;
156 $self->{_m} = new Math::BigInt 0;
157 $self->{sign} = $nan;
164 # create a bigfloat '+-inf', if given a BigFloat, set it to '+-inf'
166 my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
168 $self = $class if !defined $self;
171 my $c = $self; $self = {}; bless $self, $c;
173 $self->{_e} = new Math::BigInt 0;
174 $self->{_m} = new Math::BigInt 0;
175 $self->{sign} = $sign.'inf';
182 # create a bigfloat '+0', if given a BigFloat, set it to 0
184 $self = $class if !defined $self;
187 my $c = $self; $self = {}; bless $self, $c;
189 $self->{_m} = new Math::BigInt 0;
190 $self->{_e} = new Math::BigInt 1;
196 ##############################################################################
197 # string conversation
201 # (ref to BFLOAT or num_str ) return num_str
202 # Convert number from internal format to (non-scientific) string format.
203 # internal format is always normalized (no leading zeros, "-0" => "+0")
205 my ($self,$x) = objectify(1,@_);
207 #return "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
208 #return "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
209 return $x->{sign} if $x->{sign} !~ /^[+-]$/;
210 return '0' if $x->is_zero();
212 my $es = $x->{_m}->bstr();
213 if ($x->{_e}->is_zero())
215 $es = $x->{sign}.$es if $x->{sign} eq '-';
219 if ($x->{_e}->sign() eq '-')
221 if ($x->{_e} <= -CORE::length($es))
223 # print "style: 0.xxxx\n";
224 my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
225 $es = '0.'. ('0' x $r) . $es;
229 # print "insert '.' at $x->{_e} in '$es'\n";
230 substr($es,$x->{_e},0) = '.';
236 $es .= '0' x $x->{_e};
238 $es = $x->{sign}.$es if $x->{sign} eq '-';
244 # (ref to BFLOAT or num_str ) return num_str
245 # Convert number from internal format to scientific string format.
246 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
248 my ($self,$x) = objectify(1,@_);
250 return "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
251 return "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
252 return $x->{sign} if $x->{sign} !~ /^[+-]$/;
253 my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
255 return $x->{_m}->bstr().$sep.$x->{_e}->bstr();
260 # Make a number from a BigFloat object
261 # simple return string and let Perl's atoi() handle the rest
263 my ($self,$x) = objectify(1,@_);
267 ##############################################################################
268 # public stuff (usually prefixed with "b")
270 # really? Just for exporting them is not what I had in mind
273 # $class->SUPER::babs($class,@_);
277 # $class->SUPER::bneg($class,@_);
281 # $class->SUPER::bnot($class,@_);
286 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
287 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
288 my ($self,$x,$y) = objectify(2,@_);
289 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
292 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-';
293 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
295 return 0 if $x->is_zero() && $y->is_zero(); # 0 <=> 0
296 return -1 if $x->is_zero() && $y->{sign} eq '+'; # 0 <=> +y
297 return 1 if $y->is_zero() && $x->{sign} eq '+'; # +x <=> 0
299 # adjust so that exponents are equal
300 my $lx = $x->{_m}->length() + $x->{_e};
301 my $ly = $y->{_m}->length() + $y->{_e};
302 # print "x $x y $y lx $lx ly $ly\n";
303 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
304 # print "$l $x->{sign}\n";
305 return $l if $l != 0;
307 # lens are equal, so compare mantissa, if equal, compare exponents
308 # this assumes normaized numbers (no trailing zeros etc)
309 my $rc = $x->{_m} <=> $y->{_m} || $x->{_e} <=> $y->{_e};
310 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
316 # Compares 2 values, ignoring their signs.
317 # Returns one of undef, <0, =0, >0. (suitable for sort)
318 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
319 my ($self,$x,$y) = objectify(2,@_);
320 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
322 # signs are ignored, so check length
323 # length(x) is length(m)+e aka length of non-fraction part
324 # the longer one is bigger
325 my $l = $x->length() - $y->length();
327 return $l if $l != 0;
328 #print "equal lengths\n";
330 # if both are equal long, make full compare
331 # first compare only the mantissa
332 # if mantissa are equal, compare fractions
334 return $x->{_m} <=> $y->{_m} || $x->{_e} <=> $y->{_e};
339 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
340 # return result as BFLOAT
342 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
344 #print "add $x ",ref($x)," $y ",ref($y),"\n";
345 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
347 # speed: no add for 0+y or x+0
348 return $x if $y->is_zero(); # x+0
349 if ($x->is_zero()) # 0+y
351 # make copy, clobbering up x (modify in place!)
352 $x->{_e} = $y->{_e}->copy();
353 $x->{_m} = $y->{_m}->copy();
354 $x->{sign} = $y->{sign} || $nan;
355 return $x->round($a,$p,$r,$y);
358 # take lower of the two e's and adapt m1 to it to match m2
359 my $e = $y->{_e}; $e = Math::BigInt::bzero() if !defined $e; # if no BFLOAT
361 my $add = $y->{_m}->copy();
365 #print "\$x->{_m}: $x->{_m} ";
366 #print "\$x->{_e}: $x->{_e}\n";
367 my $e1 = $e->copy()->babs();
368 $x->{_m} *= (10 ** $e1);
369 $x->{_e} += $e; # need the sign of e
370 #$x->{_m} += $y->{_m};
371 #print "\$x->{_m}: $x->{_m} ";
372 #print "\$x->{_e}: $x->{_e}\n";
377 #print "\$x->{_m}: $x->{_m} \$y->{_m}: $y->{_m} \$e: $e ",ref($e),"\n";
379 #$x->{_m} += $y->{_m} * (10 ** $e);
380 #print "\$x->{_m}: $x->{_m}\n";
382 # else: both e are same, so leave them
383 #print "badd $x->{sign}$x->{_m} + $y->{sign}$add\n";
385 $x->{_m}->{sign} = $x->{sign};
386 $add->{sign} = $y->{sign};
390 $x->{sign} = $x->{_m}->{sign};
391 $x->{_m}->{sign} = '+';
392 return $x->round($a,$p,$r,$y);
397 # (BINT or num_str, BINT or num_str) return num_str
398 # subtract second arg from first, modify first
399 my ($self,$x,$y) = objectify(2,@_);
402 $x->badd($y->bneg()); # badd does not leave internal zeros
403 $y->bneg(); # refix y, assumes no one reads $y in between
409 # increment arg by one
410 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,@_);
420 $x->badd($self->_one('-'))->round($a,$p,$r);
425 # (BINT or num_str, BINT or num_str) return BINT
426 # does not modify arguments, but returns new object
427 # Lowest Common Multiplicator
430 my ($self,@arg) = objectify(0,@_);
431 my $x = $self->new(shift @arg);
432 while (@arg) { $x = _lcm($x,shift @arg); }
438 # (BINT or num_str, BINT or num_str) return BINT
439 # does not modify arguments, but returns new object
440 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
443 my ($self,@arg) = objectify(0,@_);
444 my $x = $self->new(shift @arg);
445 while (@arg) { $x = _gcd($x,shift @arg); }
451 # return true if arg (BINT or num_str) is zero (array '+', '0')
452 my $x = shift; $x = $class->new($x) unless ref $x;
453 #my ($self,$x) = objectify(1,@_);
455 return ($x->{sign} ne $nan && $x->{_m}->is_zero());
460 # return true if arg (BINT or num_str) is +1 (array '+', '1')
461 # or -1 if signis given
462 my $x = shift; $x = $class->new($x) unless ref $x;
463 #my ($self,$x) = objectify(1,@_);
464 my $sign = $_[2] || '+';
465 return ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
470 # return true if arg (BINT or num_str) is odd or -1 if even
471 my $x = shift; $x = $class->new($x) unless ref $x;
472 #my ($self,$x) = objectify(1,@_);
473 return ($x->{sign} ne $nan && $x->{_e}->is_zero() && $x->{_m}->is_odd());
478 # return true if arg (BINT or num_str) is even or -1 if odd
479 my $x = shift; $x = $class->new($x) unless ref $x;
480 #my ($self,$x) = objectify(1,@_);
481 return 0 if $x->{sign} eq $nan; # NaN isn't
482 return 1 if $x->{_m}->is_zero(); # 0 is
483 return ($x->{_e}->is_zero() && $x->{_m}->is_even());
488 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
489 # (BINT or num_str, BINT or num_str) return BINT
490 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
493 #print "mul $x->{_m}e$x->{_e} $y->{_m}e$y->{_e}\n";
494 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
497 # aEb * cEd = (a*c)E(b+d)
498 $x->{_m} = $x->{_m} * $y->{_m};
499 #print "m: $x->{_m}\n";
500 $x->{_e} = $x->{_e} + $y->{_e};
501 #print "e: $x->{_m}\n";
503 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
504 #print "s: $x->{sign}\n";
505 return $x->round($a,$p,$r,$y);
510 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
511 # (BFLOAT,BFLOAT) (quo,rem) or BINT (only rem)
512 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
514 return wantarray ? ($x->bnan(),bnan()) : $x->bnan()
515 if ($x->{sign} eq $nan || $y->is_nan() || $y->is_zero());
517 # we need to limit the accuracy to protect against overflow
518 my ($scale,$mode) = $x->_scale_a($accuracy,$rnd_mode,$a,$r); # ignore $p
519 my $add = 1; # for proper rounding
523 $fallback = 1; $scale = $div_scale; # simulate old behaviour
525 #print "div_scale $div_scale\n";
526 my $lx = $x->{_m}->length();
527 $scale = $lx if $lx > $scale;
528 my $ly = $y->{_m}->length();
529 $scale = $ly if $ly > $scale;
530 #print "scale $scale $lx $ly\n";
531 #$scale = $scale - $lx + $ly;
532 #print "scale $scale\n";
533 $scale += $add; # calculate some more digits for proper rounding
535 # print "bdiv $x $y scale $scale xl $lx yl $ly\n";
537 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
539 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
541 # check for / +-1 ( +/- 1E0)
544 return wantarray ? ($x,$self->bzero()) : $x;
547 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
548 #print "self: $self x: $x ref(x) ", ref($x)," m: $x->{_m}\n";
549 # my $scale_10 = 10 ** $scale; $x->{_m}->bmul($scale_10);
550 $x->{_m}->blsft($scale,10);
551 #print "m: $x->{_m}\n";
552 $x->{_m}->bdiv( $y->{_m} ); # a/c
553 #print "m: $x->{_m}\n";
554 #print "e: $x->{_e} $y->{_e}",$scale,"\n";
555 $x->{_e}->bsub($y->{_e}); # b-d
556 #print "e: $x->{_e}\n";
557 $x->{_e}->bsub($scale); # correct for 10**scale
558 #print "e: $x->{_e}\n";
559 $x->bnorm(); # remove trailing zeros
561 # print "round $x to -$scale (-$add) mode $mode\n";
562 #print "$x ",scalar ref($x), "=> $t",scalar ref($t),"\n";
565 $scale -= $add; $x->round($scale,undef,$r); # round to less
569 return $x->round($a,$p,$r,$y);
573 my $rem = $x->copy();
574 $rem->bmod($y,$a,$p,$r);
575 return ($x,$rem->round($scale,undef,$r)) if $fallback;
576 return ($x,$rem->round($a,$p,$r,$y));
583 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
584 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
586 return $x->bnan() if ($x->{sign} eq $nan || $y->is_nan() || $y->is_zero());
587 return $x->bzero() if $y->is_one();
589 # XXX tels: not done yet
590 return $x->round($a,$p,$r,$y);
595 # calculate square root
596 # this should use a different test to see wether the accuracy we want is...
597 my ($self,$x,$a,$p,$r) = objectify(1,@_);
599 # we need to limit the accuracy to protect against overflow
600 my ($scale,$mode) = $x->_scale_a($accuracy,$rnd_mode,$a,$r); # ignore $p
601 $scale = $div_scale if (!defined $scale); # simulate old behaviour
602 # print "scale $scale\n";
604 return $x->bnan() if ($x->sign() eq '-') || ($x->sign() eq $nan);
605 return $x if $x->is_zero() || $x == 1;
607 my $len = $x->{_m}->length();
608 $scale = $len if $scale < $len;
609 print "scale $scale\n";
610 $scale += 1; # because we need more than $scale to later round
611 my $e = Math::BigFloat->new("1E-$scale"); # make test variable
612 return $x->bnan() if $e->sign() eq 'NaN';
614 # print "$scale $e\n";
616 my $gs = Math::BigFloat->new(100); # first guess
617 my $org = $x->copy();
619 # start with some reasonable guess
620 #$x *= 10 ** ($len - $org->{_e});
622 #my $gs = Math::BigFloat->new(1);
623 # print "first guess: $gs (x $x)\n";
627 my $two = Math::BigFloat->new(2);
628 $x = Math::BigFloat->new($x) if ref($x) ne $class; # promote BigInts
633 return $x->bnan() if $gs->is_zero();
635 #print "$y / $gs = ",$r," ref(\$r) ",ref($r),"\n";
636 my $r = $y->copy(); $r->bdiv($gs,$scale); # $scale);
638 $x->bdiv($two,$scale); # $scale *= 2;
639 $diff = $x->copy()->bsub($gs)->babs();
640 #print "gs: $gs x: $x \n";
642 # print "$x $org $scale $gs\n";
645 #$x += $y->bdiv($x, $scale); # need only $gs scale
648 print "x $x diff $diff $e\n";
650 $x->bnorm($scale-1,undef,$mode);
655 # set to a specific 'small' value, internal usage
659 $x->{sign} = $nan, return if $v !~ /^[-+]?[0-9]+$/;
660 $x->{_m}->{value} = [abs($v)];
661 $x->{_e}->{value} = [0];
662 $x->{sign} = '+'; $x->{sign} = '-' if $v < 0;
668 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
669 # compute power of two numbers, second arg is used as integer
670 # modifies first argument
672 my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
674 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
675 return $x->_one() if $y->is_zero();
676 return $x if $x->is_one() || $y->is_one();
677 my $y1 = $y->as_number(); # make bigint
680 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
681 return $y1->is_odd() ? $x : $x->_set(1); # $x->babs() would work to
683 return $x if $x->is_zero() && $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
684 # 0 ** -y => 1 / (0 ** y) => / 0!
685 return $x->bnan() if $x->is_zero() && $y->{sign} eq '-';
687 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
691 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
693 if ($y->{sign} eq '-')
695 # modify $x in place!
696 my $z = $x->copy(); $x->_set(1);
697 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
699 return $x->round($a,$p,$r,$y);
702 ###############################################################################
707 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
708 # $n == 0 means round to integer
709 # expects and returns normalized numbers!
710 my $x = shift; $x = $class->new($x) unless ref $x;
712 return $x if $x->modify('bfround');
714 my ($scale,$mode) = $x->_scale_p($precision,$rnd_mode,@_);
715 return $x if !defined $scale; # no-op
717 # print "MBF bfround $x to scale $scale mode $mode\n";
718 return $x if $x->is_nan() or $x->is_zero();
722 # print "bfround scale $scale e $x->{_e}\n";
723 # round right from the '.'
724 return $x if $x->{_e} >= 0; # nothing to round
725 $scale = -$scale; # positive for simplicity
726 my $len = $x->{_m}->length(); # length of mantissa
727 my $dad = -$x->{_e}; # digits after dot
728 my $zad = 0; # zeros after dot
729 $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
730 # print "scale $scale dad $dad zad $zad len $len\n";
732 # number bsstr len zad dad
734 # 0.0123 123e-4 3 1 4
737 # 1.2345 12345e-4 5 0 4
739 # do not round after/right of the $dad
740 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
742 # round to zero if rounding inside the $zad, but not for last zero like:
743 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
746 $x->{_m} = Math::BigInt->new(0);
747 $x->{_e} = Math::BigInt->new(1);
751 if ($scale == $zad) # for 0.006, scale -2 and trunc
757 # adjust round-point to be inside mantissa
760 $scale = $scale-$zad;
764 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
765 $scale = $dbd+$scale;
768 # print "round to $x->{_m} to $scale\n";
772 # 123 => 100 means length(123) = 3 - $scale (2) => 1
774 # calculate digits before dot
775 my $dbt = $x->{_m}->length(); $dbt += $x->{_e} if $x->{_e}->sign() eq '-';
776 if (($scale > $dbt) && ($dbt < 0))
778 # if not enough digits before dot, round to zero
779 $x->{_m} = Math::BigInt->new(0);
780 $x->{_e} = Math::BigInt->new(1);
784 if (($scale >= 0) && ($dbt == 0))
786 # 0.49->bfround(1): scale == 1, dbt == 0: => 0.0
787 # 0.51->bfround(0): scale == 0, dbt == 0: => 1.0
788 # 0.5->bfround(0): scale == 0, dbt == 0: => 0
789 # 0.05->bfround(0): scale == 0, dbt == 0: => 0
790 # print "$scale $dbt $x->{_m}\n";
791 $scale = -$x->{_m}->length();
795 # correct by subtracting scale
796 $scale = $dbt - $scale;
800 $scale = $x->{_m}->length() - $scale;
803 #print "using $scale for $x->{_m} with '$mode'\n";
804 # pass sign to bround for '+inf' and '-inf' rounding modes
805 $x->{_m}->{sign} = $x->{sign};
806 $x->{_m}->bround($scale,$mode);
807 $x->{_m}->{sign} = '+'; # fix sign back
813 # accuracy: preserve $N digits, and overwrite the rest with 0's
814 my $x = shift; $x = $class->new($x) unless ref $x;
815 my ($scale,$mode) = $x->_scale_a($accuracy,$rnd_mode,@_);
816 return $x if !defined $scale; # no-op
818 return $x if $x->modify('bround');
820 # print "bround $scale $mode\n";
821 # 0 => return all digits, scale < 0 makes no sense
822 return $x if ($scale <= 0);
823 return $x if $x->is_nan() or $x->is_zero(); # never round a 0
825 # if $e longer than $m, we have 0.0000xxxyyy style number, and must
826 # subtract the delta from scale, to simulate keeping the zeros
827 # -5 +5 => 1; -10 +5 => -4
828 my $delta = $x->{_e} + $x->{_m}->length() + 1;
829 # removed by tlr, since causes problems with fraction tests:
830 # $scale += $delta if $delta < 0;
832 # if we should keep more digits than the mantissa has, do nothing
833 return $x if $x->{_m}->length() <= $scale;
835 # pass sign to bround for '+inf' and '-inf' rounding modes
836 $x->{_m}->{sign} = $x->{sign};
837 $x->{_m}->bround($scale,$mode); # round mantissa
838 $x->{_m}->{sign} = '+'; # fix sign back
839 return $x->bnorm(); # del trailing zeros gen. by bround()
844 # return integer less or equal then $x
845 my ($self,$x,$a,$p,$r) = objectify(1,@_);
847 return $x if $x->modify('bfloor');
849 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
851 # if $x has digits after dot
852 if ($x->{_e}->{sign} eq '-')
854 $x->{_m}->brsft(-$x->{_e},10);
856 $x-- if $x->{sign} eq '-';
858 return $x->round($a,$p,$r);
863 # return integer greater or equal then $x
864 my ($self,$x,$a,$p,$r) = objectify(1,@_);
866 return $x if $x->modify('bceil');
867 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
869 # if $x has digits after dot
870 if ($x->{_e}->{sign} eq '-')
872 $x->{_m}->brsft(-$x->{_e},10);
874 $x++ if $x->{sign} eq '+';
876 return $x->round($a,$p,$r);
879 ###############################################################################
883 # going trough AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
888 # make fxxx and bxxx work
890 my $name = $AUTOLOAD;
892 $name =~ s/.*:://; # split package
894 if (!method_valid($name))
898 #&{$class."::SUPER->$name"}(@_);
899 # delayed load of Carp and avoid recursion
901 Carp::croak ("Can't call $class\-\>$name, not a valid method");
904 my $bname = $name; $bname =~ s/^f/b/;
905 *{$class."\:\:$name"} = \&$bname;
911 # return a copy of the exponent
913 $self = $class->new($self) unless ref $self;
915 return bnan() if $self->is_nan();
916 return $self->{_e}->copy();
921 # return a copy of the mantissa
923 $self = $class->new($self) unless ref $self;
925 return bnan() if $self->is_nan();
926 my $m = $self->{_m}->copy(); # faster than going via bstr()
927 $m->bneg() if $self->{sign} eq '-';
934 # return a copy of both the exponent and the mantissa
936 $self = $class->new($self) unless ref $self;
938 return (bnan(),bnan()) if $self->is_nan();
939 my $m = $self->{_m}->copy(); # faster than going via bstr()
940 $m->bneg() if $self->{sign} eq '-';
941 return ($m,$self->{_e}->copy());
944 ##############################################################################
945 # private stuff (internal use only)
949 # internal speedup, set argument to 1, or create a +/- 1
950 # uses internal knowledge about MBI, thus (bad)
952 my $x = $self->bzero();
953 $x->{_m}->{value} = [ 1 ]; $x->{_m}->{sign} = '+';
954 $x->{_e}->{value} = [ 0 ]; $x->{_e}->{sign} = '+';
955 $x->{sign} = shift || '+';
962 #print "import $self\n";
963 for ( my $i = 0; $i < @_ ; $i++ )
965 if ( $_[$i] eq ':constant' )
967 # this rest causes overlord er load to step in
968 # print "overload @_\n";
969 overload::constant float => sub { $self->new(shift); };
970 splice @_, $i, 1; last;
973 # any non :constant stuff is handled by our parent, Exporter
974 # even if @_ is empty, to give it a chance
975 #$self->SUPER::import(@_); # does not work (would call MBI)
976 $self->export_to_level(1,$self,@_); # need this instead
981 # adjust m and e so that m is smallest possible
982 # round number according to accuracy and precision settings
985 return $x if $x->is_nan();
987 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
990 $x->{_m}->brsft($zeros,10); $x->{_e} += $zeros;
992 # for something like 0Ey, set y to 1
993 $x->{_e}->bzero()->binc() if $x->{_m}->is_zero();
994 return $x->SUPER::bnorm(@_); # call MBI bnorm for round
997 ##############################################################################
998 # internal calculation routines
1002 # return a bigint representation of this BigFloat number
1003 my ($self,$x) = objectify(1,@_);
1006 if ($x->{_e}->is_zero())
1008 $z = $x->{_m}->copy();
1009 $z->{sign} = $x->{sign};
1015 my $y = $x->{_m} / ($ten ** $x->{_e});
1017 $y->{sign} = $x->{sign};
1020 $z = $x->{_m} * ($ten ** $x->{_e});
1021 $z->{sign} = $x->{sign};
1027 my $x = shift; $x = $class->new($x) unless ref $x;
1029 my $len = $x->{_m}->length();
1030 $len += $x->{_e} if $x->{_e}->sign() eq '+';
1033 my $t = Math::BigInt::bzero();
1034 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1045 Math::BigFloat - Arbitrary size floating point math package
1052 $x = Math::BigInt->new($str); # defaults to 0
1053 $nan = Math::BigInt->bnan(); # create a NotANumber
1054 $zero = Math::BigInt->bzero();# create a "+0"
1057 $x->is_zero(); # return whether arg is zero or not
1058 $x->is_one(); # return true if arg is +1
1059 $x->is_one('-'); # return true if arg is -1
1060 $x->is_odd(); # return true if odd, false for even
1061 $x->is_even(); # return true if even, false for odd
1062 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
1063 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
1064 $x->sign(); # return the sign, either +,- or NaN
1066 # The following all modify their first argument:
1069 $x->bzero(); # set $i to 0
1070 $x->bnan(); # set $i to NaN
1072 $x->bneg(); # negation
1073 $x->babs(); # absolute value
1074 $x->bnorm(); # normalize (no-op)
1075 $x->bnot(); # two's complement (bit wise not)
1076 $x->binc(); # increment x by 1
1077 $x->bdec(); # decrement x by 1
1079 $x->badd($y); # addition (add $y to $x)
1080 $x->bsub($y); # subtraction (subtract $y from $x)
1081 $x->bmul($y); # multiplication (multiply $x by $y)
1082 $x->bdiv($y); # divide, set $i to quotient
1083 # return (quo,rem) or quo if scalar
1085 $x->bmod($y); # modulus
1086 $x->bpow($y); # power of arguments (a**b)
1087 $x->blsft($y); # left shift
1088 $x->brsft($y); # right shift
1089 # return (quo,rem) or quo if scalar
1091 $x->band($y); # bit-wise and
1092 $x->bior($y); # bit-wise inclusive or
1093 $x->bxor($y); # bit-wise exclusive or
1094 $x->bnot(); # bit-wise not (two's complement)
1096 $x->bround($N); # accuracy: preserver $N digits
1097 $x->bfround($N); # precision: round to the $Nth digit
1099 # The following do not modify their arguments:
1101 bgcd(@values); # greatest common divisor
1102 blcm(@values); # lowest common multiplicator
1104 $x->bstr(); # return string
1105 $x->bsstr(); # return string in scientific notation
1107 $x->exponent(); # return exponent as BigInt
1108 $x->mantissa(); # return mantissa as BigInt
1109 $x->parts(); # return (mantissa,exponent) as BigInt
1111 $x->length(); # number of digits (w/o sign and '.')
1112 ($l,$f) = $x->length(); # number of digits, and length of fraction
1116 All operators (inlcuding basic math operations) are overloaded if you
1117 declare your big floating point numbers as
1119 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1121 Operations with overloaded operators preserve the arguments, which is
1122 exactly what you expect.
1124 =head2 Canonical notation
1126 Input to these routines are either BigFloat objects, or strings of the
1127 following four forms:
1141 C</^[+-]\d+E[+-]?\d+$/>
1145 C</^[+-]\d*\.\d+E[+-]?\d+$/>
1149 all with optional leading and trailing zeros and/or spaces. Additonally,
1150 numbers are allowed to have an underscore between any two digits.
1152 Empty strings as well as other illegal numbers results in 'NaN'.
1154 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1155 are always stored in normalized form. On a string, it creates a BigFloat
1160 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1162 The string output will always have leading and trailing zeros stripped and drop
1163 a plus sign. C<bstr()> will give you always the form with a decimal point,
1164 while C<bsstr()> (for scientific) gives you the scientific notation.
1166 Input bstr() bsstr()
1168 ' -123 123 123' '-123123123' '-123123123E0'
1169 '00.0123' '0.0123' '123E-4'
1170 '123.45E-2' '1.2345' '12345E-4'
1171 '10E+3' '10000' '1E4'
1173 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1174 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1175 return either undef, <0, 0 or >0 and are suited for sort.
1177 Actual math is done by using BigInts to represent the mantissa and exponent.
1178 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1179 represent the result when input arguments are not numbers, as well as
1180 the result of dividing by zero.
1182 =head2 C<mantissa()>, C<exponent()> and C<parts()>
1184 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1185 as BigInts such that:
1187 $m = $x->mantissa();
1188 $e = $x->exponent();
1189 $y = $m * ( 10 ** $e );
1190 print "ok\n" if $x == $y;
1192 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
1194 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
1196 Currently the mantissa is reduced as much as possible, favouring higher
1197 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
1198 This might change in the future, so do not depend on it.
1200 =head2 Accuracy vs. Precision
1202 See also: L<Rounding|Rounding>.
1204 Math::BigFloat supports both precision and accuracy. (here should follow
1205 a short description of both).
1207 Precision: digits after the '.', laber, schwad
1208 Accuracy: Significant digits blah blah
1210 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
1211 a operation consumes all resources, each operation produces no more than
1212 C<Math::BigFloat::precision()> digits.
1214 In case the result of one operation has more precision than specified,
1215 it is rounded. The rounding mode taken is either the default mode, or the one
1216 supplied to the operation after the I<scale>:
1218 $x = Math::BigFloat->new(2);
1219 Math::BigFloat::precision(5); # 5 digits max
1220 $y = $x->copy()->bdiv(3); # will give 0.66666
1221 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1222 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
1223 Math::BigFloat::round_mode('zero');
1224 $y = $x->copy()->bdiv(3,6); # will give 0.666666
1230 =item ffround ( +$scale ) rounds to the $scale'th place left from the '.', counting from the dot. The first digit is numbered 1.
1232 =item ffround ( -$scale ) rounds to the $scale'th place right from the '.', counting from the dot
1234 =item ffround ( 0 ) rounds to an integer
1236 =item fround ( +$scale ) preserves accuracy to $scale digits from the left (aka significant digits) and paddes the rest with zeros. If the number is between 1 and -1, the significant digits count from the first non-zero after the '.'
1238 =item fround ( -$scale ) and fround ( 0 ) are a no-ops
1242 All rounding functions take as a second parameter a rounding mode from one of
1243 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
1245 The default rounding mode is 'even'. By using
1246 C<< Math::BigFloat::round_mode($rnd_mode); >> you can get and set the default
1247 mode for subsequent rounding. The usage of C<$Math::BigFloat::$rnd_mode> is
1248 no longer supported.
1249 The second parameter to the round functions then overrides the default
1252 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
1253 'trunc' as rounding mode to make it equivalent to:
1258 You can override this by passing the desired rounding mode as parameter to
1261 $x = Math::BigFloat->new(2.5);
1262 $y = $x->as_number('odd'); # $y = 3
1266 use Math::BigFloat qw(bstr bint);
1268 $x = bstr("1234") # string "1234"
1269 $x = "$x"; # same as bstr()
1270 $x = bneg("1234") # BigFloat "-1234"
1271 $x = Math::BigFloat->bneg("1234"); # BigFloat "1234"
1272 $x = Math::BigFloat->babs("-12345"); # BigFloat "12345"
1273 $x = Math::BigFloat->bnorm("-0 00"); # BigFloat "0"
1274 $x = bint(1) + bint(2); # BigFloat "3"
1275 $x = bint(1) + "2"; # ditto (auto-BigFloatify of "2")
1276 $x = bint(1); # BigFloat "1"
1277 $x = $x + 5 / 2; # BigFloat "3"
1278 $x = $x ** 3; # BigFloat "27"
1279 $x *= 2; # BigFloat "54"
1280 $x = new Math::BigFloat; # BigFloat "0"
1281 $x--; # BigFloat "-1"
1283 =head1 Autocreating constants
1285 After C<use Math::BigFloat ':constant'> all the floating point constants
1286 in the given scope are converted to C<Math::BigFloat>. This conversion
1287 happens at compile time.
1291 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
1293 prints the value of C<2E-100>. Note that without conversion of
1294 constants the expression 2E-100 will be calculated as normal floating point
1299 Greatly enhanced ;o)
1308 The following does not work yet:
1310 $m = $x->mantissa();
1311 $e = $x->exponent();
1312 $y = $m * ( 10 ** $e );
1313 print "ok\n" if $x == $y;
1317 There is no fmod() function yet.
1325 =item stringify, bstr()
1327 Both stringify and bstr() now drop the leading '+'. The old code would return
1328 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
1329 reasoning and details.
1333 The following will probably not do what you expect:
1335 print $c->bdiv(123.456),"\n";
1337 It prints both quotient and reminder since print works in list context. Also,
1338 bdiv() will modify $c, so be carefull. You probably want to use
1340 print $c / 123.456,"\n";
1341 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
1345 =item Modifying and =
1349 $x = Math::BigFloat->new(5);
1352 It will not do what you think, e.g. making a copy of $x. Instead it just makes
1353 a second reference to the B<same> object and stores it in $y. Thus anything
1354 that modifies $x will modify $y, and vice versa.
1357 print "$x, $y\n"; # prints '10, 10'
1359 If you want a true copy of $x, use:
1363 See also the documentation in L<overload> regarding C<=>.
1367 C<bpow()> now modifies the first argument, unlike the old code which left
1368 it alone and only returned the result. This is to be consistent with
1369 C<badd()> etc. The first will modify $x, the second one won't:
1371 print bpow($x,$i),"\n"; # modify $x
1372 print $x->bpow($i),"\n"; # ditto
1373 print $x ** $i,"\n"; # leave $x alone
1379 This program is free software; you may redistribute it and/or modify it under
1380 the same terms as Perl itself.
1384 Mark Biggar, overloaded interface by Ilya Zakharevich.
1385 Completely rewritten by Tels http://bloodgate.com in 2001.