1 package Math::BigFloat;
4 # Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
7 # The following hash values are internally used:
8 # _e : exponent (ref to $CALC object)
9 # _m : mantissa (ref to $CALC object)
11 # sign : +,-,+inf,-inf, or "NaN" if not a number
19 @ISA = qw(Exporter Math::BigInt);
22 # $_trap_inf/$_trap_nan are internal and should never be accessed from outside
23 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode
24 $upgrade $downgrade $_trap_nan $_trap_inf/;
25 my $class = "Math::BigFloat";
28 '<=>' => sub { my $rc = $_[2] ?
29 ref($_[0])->bcmp($_[1],$_[0]) :
30 ref($_[0])->bcmp($_[0],$_[1]);
31 $rc = 1 unless defined $rc;
34 # we need '>=' to get things like "1 >= NaN" right:
35 '>=' => sub { my $rc = $_[2] ?
36 ref($_[0])->bcmp($_[1],$_[0]) :
37 ref($_[0])->bcmp($_[0],$_[1]);
38 # if there was a NaN involved, return false
39 return '' unless defined $rc;
42 'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
45 ##############################################################################
46 # global constants, flags and assorted stuff
48 # the following are public, but their usage is not recommended. Use the
49 # accessor methods instead.
51 # class constants, use Class->constant_name() to access
52 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
59 # the package we are using for our private parts, defaults to:
60 # Math::BigInt->config()->{lib}
61 my $MBI = 'Math::BigInt::FastCalc';
63 # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
65 # the same for infinity
68 # constant for easier life
71 my $IMPORT = 0; # was import() called yet? used to make require work
73 # some digits of accuracy for blog(undef,10); which we use in blog() for speed
75 '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
76 my $LOG_10_A = length($LOG_10)-1;
79 '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
80 my $LOG_2_A = length($LOG_2)-1;
81 my $HALF = '0.5'; # made into an object if nec.
83 ##############################################################################
84 # the old code had $rnd_mode, so we need to support it, too
86 sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
87 sub FETCH { return $round_mode; }
88 sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
92 # when someone sets $rnd_mode, we catch this and check the value to see
93 # whether it is valid or not.
94 $rnd_mode = 'even'; tie $rnd_mode, 'Math::BigFloat';
96 # we need both of them in this package:
97 *as_int = \&as_number;
100 ##############################################################################
103 # valid method aliases for AUTOLOAD
104 my %methods = map { $_ => 1 }
105 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
106 fint facmp fcmp fzero fnan finf finc fdec ffac fneg
107 fceil ffloor frsft flsft fone flog froot fexp
109 # valid methods that can be handed up (for AUTOLOAD)
110 my %hand_ups = map { $_ => 1 }
111 qw / is_nan is_inf is_negative is_positive is_pos is_neg
112 accuracy precision div_scale round_mode fabs fnot
113 objectify upgrade downgrade
118 sub _method_alias { exists $methods{$_[0]||''}; }
119 sub _method_hand_up { exists $hand_ups{$_[0]||''}; }
122 ##############################################################################
127 # create a new BigFloat object from a string or another bigfloat object.
130 # sign => sign (+/-), or "NaN"
132 my ($class,$wanted,@r) = @_;
134 # avoid numify-calls by not using || on $wanted!
135 return $class->bzero() if !defined $wanted; # default to 0
136 return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
138 $class->import() if $IMPORT == 0; # make require work
140 my $self = {}; bless $self, $class;
141 # shortcut for bigints and its subclasses
142 if ((ref($wanted)) && UNIVERSAL::can( $wanted, "as_number"))
144 $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy
145 $self->{_e} = $MBI->_zero();
147 $self->{sign} = $wanted->sign();
148 return $self->bnorm();
150 # else: got a string or something maskerading as number (with overload)
152 # handle '+inf', '-inf' first
153 if ($wanted =~ /^[+-]?inf\z/)
155 return $downgrade->new($wanted) if $downgrade;
157 $self->{sign} = $wanted; # set a default sign for bstr()
158 return $self->binf($wanted);
161 # shortcut for simple forms like '12' that neither have trailing nor leading
163 if ($wanted =~ /^([+-]?)([1-9][0-9]*[1-9])$/)
165 $self->{_e} = $MBI->_zero();
167 $self->{sign} = $1 || '+';
168 $self->{_m} = $MBI->_new($2);
169 return $self->round(@r) if !$downgrade;
172 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
178 Carp::croak ("$wanted is not a number initialized to $class");
181 return $downgrade->bnan() if $downgrade;
183 $self->{_e} = $MBI->_zero();
185 $self->{_m} = $MBI->_zero();
186 $self->{sign} = $nan;
190 # make integer from mantissa by adjusting exp, then convert to int
191 $self->{_e} = $MBI->_new($$ev); # exponent
192 $self->{_es} = $$es || '+';
193 my $mantissa = "$$miv$$mfv"; # create mant.
194 $mantissa =~ s/^0+(\d)/$1/; # strip leading zeros
195 $self->{_m} = $MBI->_new($mantissa); # create mant.
197 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
198 if (CORE::length($$mfv) != 0)
200 my $len = $MBI->_new( CORE::length($$mfv));
201 ($self->{_e}, $self->{_es}) =
202 _e_sub ($self->{_e}, $len, $self->{_es}, '+');
204 # we can only have trailing zeros on the mantissa if $$mfv eq ''
207 # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
208 # because that is faster, especially when _m is not stored in base 10.
209 my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/;
212 my $z = $MBI->_new($zeros);
213 # turn '120e2' into '12e3'
214 $MBI->_rsft ( $self->{_m}, $z, 10);
215 ($self->{_e}, $self->{_es}) =
216 _e_add ( $self->{_e}, $z, $self->{_es}, '+');
219 $self->{sign} = $$mis;
221 # for something like 0Ey, set y to 1, and -0 => +0
222 # Check $$miv for being '0' and $$mfv eq '', because otherwise _m could not
223 # have become 0. That's faster than to call $MBI->_is_zero().
224 $self->{sign} = '+', $self->{_e} = $MBI->_one()
225 if $$miv eq '0' and $$mfv eq '';
227 return $self->round(@r) if !$downgrade;
229 # if downgrade, inf, NaN or integers go down
231 if ($downgrade && $self->{_es} eq '+')
233 if ($MBI->_is_zero( $self->{_e} ))
235 return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
237 return $downgrade->new($self->bsstr());
239 $self->bnorm()->round(@r); # first normalize, then round
247 # if two arguments, the first one is the class to "swallow" subclasses
255 return unless ref($x); # only for objects
257 my $self = {}; bless $self,$c;
259 $self->{sign} = $x->{sign};
260 $self->{_es} = $x->{_es};
261 $self->{_m} = $MBI->_copy($x->{_m});
262 $self->{_e} = $MBI->_copy($x->{_e});
263 $self->{_a} = $x->{_a} if defined $x->{_a};
264 $self->{_p} = $x->{_p} if defined $x->{_p};
270 # used by parent class bone() to initialize number to NaN
276 my $class = ref($self);
277 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
280 $IMPORT=1; # call our import only once
281 $self->{_m} = $MBI->_zero();
282 $self->{_e} = $MBI->_zero();
288 # used by parent class bone() to initialize number to +-inf
294 my $class = ref($self);
295 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
298 $IMPORT=1; # call our import only once
299 $self->{_m} = $MBI->_zero();
300 $self->{_e} = $MBI->_zero();
306 # used by parent class bone() to initialize number to 1
308 $IMPORT=1; # call our import only once
309 $self->{_m} = $MBI->_one();
310 $self->{_e} = $MBI->_zero();
316 # used by parent class bone() to initialize number to 0
318 $IMPORT=1; # call our import only once
319 $self->{_m} = $MBI->_zero();
320 $self->{_e} = $MBI->_one();
326 my ($self,$class) = @_;
327 return if $class =~ /^Math::BigInt/; # we aren't one of these
328 UNIVERSAL::isa($self,$class);
333 # return (later set?) configuration data as hash ref
334 my $class = shift || 'Math::BigFloat';
336 my $cfg = $class->SUPER::config(@_);
338 # now we need only to override the ones that are different from our parent
339 $cfg->{class} = $class;
344 ##############################################################################
345 # string conversation
349 # (ref to BFLOAT or num_str ) return num_str
350 # Convert number from internal format to (non-scientific) string format.
351 # internal format is always normalized (no leading zeros, "-0" => "+0")
352 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
354 if ($x->{sign} !~ /^[+-]$/)
356 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
360 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
363 my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
366 $es = $MBI->_str($x->{_m});
367 $len = CORE::length($es);
368 my $e = $MBI->_num($x->{_e});
369 $e = -$e if $x->{_es} eq '-';
373 # if _e is bigger than a scalar, the following will blow your memory
376 my $r = abs($e) - $len;
377 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
381 substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
382 $cad = -$cad if $x->{_es} eq '-';
388 $es .= '0' x $e; $len += $e; $cad = 0;
392 $es = '-'.$es if $x->{sign} eq '-';
393 # if set accuracy or precision, pad with zeros on the right side
394 if ((defined $x->{_a}) && ($not_zero))
396 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
397 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
398 $zeros = $x->{_a} - $len if $cad != $len;
399 $es .= $dot.'0' x $zeros if $zeros > 0;
401 elsif ((($x->{_p} || 0) < 0))
403 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
404 my $zeros = -$x->{_p} + $cad;
405 $es .= $dot.'0' x $zeros if $zeros > 0;
412 # (ref to BFLOAT or num_str ) return num_str
413 # Convert number from internal format to scientific string format.
414 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
415 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
417 if ($x->{sign} !~ /^[+-]$/)
419 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
422 my $sep = 'e'.$x->{_es};
423 my $sign = $x->{sign}; $sign = '' if $sign eq '+';
424 $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
429 # Make a number from a BigFloat object
430 # simple return a string and let Perl's atoi()/atof() handle the rest
431 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
435 ##############################################################################
436 # public stuff (usually prefixed with "b")
440 # (BINT or num_str) return BINT
441 # negate number or make a negated number from string
442 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
444 return $x if $x->modify('bneg');
446 # for +0 dont negate (to have always normalized +0). Does nothing for 'NaN'
447 $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
452 # XXX TODO this must be overwritten and return NaN for non-integer values
453 # band(), bior(), bxor(), too
456 # $class->SUPER::bnot($class,@_);
461 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
464 my ($self,$x,$y) = (ref($_[0]),@_);
465 # objectify is costly, so avoid it
466 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
468 ($self,$x,$y) = objectify(2,@_);
471 return $upgrade->bcmp($x,$y) if defined $upgrade &&
472 ((!$x->isa($self)) || (!$y->isa($self)));
474 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
476 # handle +-inf and NaN
477 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
478 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
479 return +1 if $x->{sign} eq '+inf';
480 return -1 if $x->{sign} eq '-inf';
481 return -1 if $y->{sign} eq '+inf';
485 # check sign for speed first
486 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
487 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
490 my $xz = $x->is_zero();
491 my $yz = $y->is_zero();
492 return 0 if $xz && $yz; # 0 <=> 0
493 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
494 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
496 # adjust so that exponents are equal
497 my $lxm = $MBI->_len($x->{_m});
498 my $lym = $MBI->_len($y->{_m});
499 # the numify somewhat limits our length, but makes it much faster
500 my ($xes,$yes) = (1,1);
501 $xes = -1 if $x->{_es} ne '+';
502 $yes = -1 if $y->{_es} ne '+';
503 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
504 my $ly = $lym + $yes * $MBI->_num($y->{_e});
505 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
506 return $l <=> 0 if $l != 0;
508 # lengths (corrected by exponent) are equal
509 # so make mantissa equal length by padding with zero (shift left)
510 my $diff = $lxm - $lym;
511 my $xm = $x->{_m}; # not yet copy it
515 $ym = $MBI->_copy($y->{_m});
516 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
520 $xm = $MBI->_copy($x->{_m});
521 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
523 my $rc = $MBI->_acmp($xm,$ym);
524 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
530 # Compares 2 values, ignoring their signs.
531 # Returns one of undef, <0, =0, >0. (suitable for sort)
534 my ($self,$x,$y) = (ref($_[0]),@_);
535 # objectify is costly, so avoid it
536 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
538 ($self,$x,$y) = objectify(2,@_);
541 return $upgrade->bacmp($x,$y) if defined $upgrade &&
542 ((!$x->isa($self)) || (!$y->isa($self)));
544 # handle +-inf and NaN's
545 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
547 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
548 return 0 if ($x->is_inf() && $y->is_inf());
549 return 1 if ($x->is_inf() && !$y->is_inf());
554 my $xz = $x->is_zero();
555 my $yz = $y->is_zero();
556 return 0 if $xz && $yz; # 0 <=> 0
557 return -1 if $xz && !$yz; # 0 <=> +y
558 return 1 if $yz && !$xz; # +x <=> 0
560 # adjust so that exponents are equal
561 my $lxm = $MBI->_len($x->{_m});
562 my $lym = $MBI->_len($y->{_m});
563 my ($xes,$yes) = (1,1);
564 $xes = -1 if $x->{_es} ne '+';
565 $yes = -1 if $y->{_es} ne '+';
566 # the numify somewhat limits our length, but makes it much faster
567 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
568 my $ly = $lym + $yes * $MBI->_num($y->{_e});
570 return $l <=> 0 if $l != 0;
572 # lengths (corrected by exponent) are equal
573 # so make mantissa equal-length by padding with zero (shift left)
574 my $diff = $lxm - $lym;
575 my $xm = $x->{_m}; # not yet copy it
579 $ym = $MBI->_copy($y->{_m});
580 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
584 $xm = $MBI->_copy($x->{_m});
585 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
587 $MBI->_acmp($xm,$ym);
592 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
593 # return result as BFLOAT
596 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
597 # objectify is costly, so avoid it
598 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
600 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
603 return $x if $x->modify('badd');
605 # inf and NaN handling
606 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
609 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
611 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
613 # +inf++inf or -inf+-inf => same, rest is NaN
614 return $x if $x->{sign} eq $y->{sign};
617 # +-inf + something => +inf; something +-inf => +-inf
618 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
622 return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
623 ((!$x->isa($self)) || (!$y->isa($self)));
625 # speed: no add for 0+y or x+0
626 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
627 if ($x->is_zero()) # 0+y
629 # make copy, clobbering up x (modify in place!)
630 $x->{_e} = $MBI->_copy($y->{_e});
631 $x->{_es} = $y->{_es};
632 $x->{_m} = $MBI->_copy($y->{_m});
633 $x->{sign} = $y->{sign} || $nan;
634 return $x->round($a,$p,$r,$y);
637 # take lower of the two e's and adapt m1 to it to match m2
639 $e = $MBI->_zero() if !defined $e; # if no BFLOAT?
640 $e = $MBI->_copy($e); # make copy (didn't do it yet)
644 ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
646 my $add = $MBI->_copy($y->{_m});
648 if ($es eq '-') # < 0
650 $MBI->_lsft( $x->{_m}, $e, 10);
651 ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
653 elsif (!$MBI->_is_zero($e)) # > 0
655 $MBI->_lsft($add, $e, 10);
657 # else: both e are the same, so just leave them
659 if ($x->{sign} eq $y->{sign})
662 $x->{_m} = $MBI->_add($x->{_m}, $add);
666 ($x->{_m}, $x->{sign}) =
667 _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
670 # delete trailing zeros, then round
671 $x->bnorm()->round($a,$p,$r,$y);
674 # sub bsub is inherited from Math::BigInt!
678 # increment arg by one
679 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
681 return $x if $x->modify('binc');
683 if ($x->{_es} eq '-')
685 return $x->badd($self->bone(),@r); # digits after dot
688 if (!$MBI->_is_zero($x->{_e})) # _e == 0 for NaN, inf, -inf
690 # 1e2 => 100, so after the shift below _m has a '0' as last digit
691 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
692 $x->{_e} = $MBI->_zero(); # normalize
694 # we know that the last digit of $x will be '1' or '9', depending on the
698 if ($x->{sign} eq '+')
700 $MBI->_inc($x->{_m});
701 return $x->bnorm()->bround(@r);
703 elsif ($x->{sign} eq '-')
705 $MBI->_dec($x->{_m});
706 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
707 return $x->bnorm()->bround(@r);
709 # inf, nan handling etc
710 $x->badd($self->bone(),@r); # badd() does round
715 # decrement arg by one
716 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
718 return $x if $x->modify('bdec');
720 if ($x->{_es} eq '-')
722 return $x->badd($self->bone('-'),@r); # digits after dot
725 if (!$MBI->_is_zero($x->{_e}))
727 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
728 $x->{_e} = $MBI->_zero(); # normalize
732 my $zero = $x->is_zero();
734 if (($x->{sign} eq '-') || $zero)
736 $MBI->_inc($x->{_m});
737 $x->{sign} = '-' if $zero; # 0 => 1 => -1
738 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
739 return $x->bnorm()->round(@r);
742 elsif ($x->{sign} eq '+')
744 $MBI->_dec($x->{_m});
745 return $x->bnorm()->round(@r);
747 # inf, nan handling etc
748 $x->badd($self->bone('-'),@r); # does round
755 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
757 return $x if $x->modify('blog');
759 # $base > 0, $base != 1; if $base == undef default to $base == e
762 # we need to limit the accuracy to protect against overflow
765 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
767 # also takes care of the "error in _find_round_parameters?" case
768 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
770 # no rounding at all, so must use fallback
771 if (scalar @params == 0)
773 # simulate old behaviour
774 $params[0] = $self->div_scale(); # and round to it as accuracy
775 $params[1] = undef; # P = undef
776 $scale = $params[0]+4; # at least four more for proper round
777 $params[2] = $r; # round mode by caller or undef
778 $fallback = 1; # to clear a/p afterwards
782 # the 4 below is empirical, and there might be cases where it is not
784 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
787 return $x->bzero(@params) if $x->is_one();
788 # base not defined => base == Euler's number e
791 # make object, since we don't feed it through objectify() to still get the
792 # case of $base == undef
793 $base = $self->new($base) unless ref($base);
794 # $base > 0; $base != 1
795 return $x->bnan() if $base->is_zero() || $base->is_one() ||
796 $base->{sign} ne '+';
797 # if $x == $base, we know the result must be 1.0
798 if ($x->bcmp($base) == 0)
800 $x->bone('+',@params);
803 # clear a/p after round, since user did not request it
804 delete $x->{_a}; delete $x->{_p};
810 # when user set globals, they would interfere with our calculation, so
811 # disable them and later re-enable them
813 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
814 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
815 # we also need to disable any set A or P on $x (_find_round_parameters took
816 # them already into account), since these would interfere, too
817 delete $x->{_a}; delete $x->{_p};
818 # need to disable $upgrade in BigInt, to avoid deep recursion
819 local $Math::BigInt::upgrade = undef;
820 local $Math::BigFloat::downgrade = undef;
822 # upgrade $x if $x is not a BigFloat (handle BigInt input)
824 if (!$x->isa('Math::BigFloat'))
826 $x = Math::BigFloat->new($x);
832 # If the base is defined and an integer, try to calculate integer result
833 # first. This is very fast, and in case the real result was found, we can
835 if (defined $base && $base->is_int() && $x->is_int())
837 my $i = $MBI->_copy( $x->{_m} );
838 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
839 my $int = Math::BigInt->bzero();
841 $int->blog($base->as_number());
843 if ($base->as_number()->bpow($int) == $x)
845 # found result, return it
846 $x->{_m} = $int->{value};
847 $x->{_e} = $MBI->_zero();
856 # base is undef, so base should be e (Euler's number), so first calculate the
857 # log to base e (using reduction by 10 (and probably 2)):
858 $self->_log_10($x,$scale);
860 # and if a different base was requested, convert it
863 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
864 # not ln, but some other base (don't modify $base)
865 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
869 # shortcut to not run through _find_round_parameters again
870 if (defined $params[0])
872 $x->bround($params[0],$params[2]); # then round accordingly
876 $x->bfround($params[1],$params[2]); # then round accordingly
880 # clear a/p after round, since user did not request it
881 delete $x->{_a}; delete $x->{_p};
884 $$abr = $ab; $$pbr = $pb;
891 # Given D (digits in decimal), compute N so that N! (N factorial) is
892 # at least D digits long. D should be at least 50.
895 # two constants for the Ramanujan estimate of ln(N!)
896 my $lg2 = log(2 * 3.14159265) / 2;
899 # D = 50 => N => 42, so L = 40 and R = 50
900 my $l = 40; my $r = $d;
902 # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
903 $l = $l->numify if ref($l);
904 $r = $r->numify if ref($r);
905 $lg2 = $lg2->numify if ref($lg2);
906 $lg10 = $lg10->numify if ref($lg10);
908 # binary search for the right value (could this be written as the reverse of lg(n!)?)
911 my $n = int(($r - $l) / 2) + $l;
913 int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
914 $ramanujan > $d ? $r = $n : $l = $n;
921 # Calculate n over k (binomial coefficient or "choose" function) as integer.
923 my ($self,$x,$y,@r) = (ref($_[0]),@_);
925 # objectify is costly, so avoid it
926 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
928 ($self,$x,$y,@r) = objectify(2,@_);
931 return $x if $x->modify('bnok');
933 return $x->bnan() if $x->is_nan() || $y->is_nan();
934 return $x->binf() if $x->is_inf();
936 my $u = $x->as_int();
937 $u->bnok($y->as_int());
939 $x->{_m} = $u->{value};
940 $x->{_e} = $MBI->_zero();
948 # Calculate e ** X (Euler's number to the power of X)
949 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
951 return $x if $x->modify('bexp');
953 return $x->binf() if $x->{sign} eq '+inf';
954 return $x->bzero() if $x->{sign} eq '-inf';
956 # we need to limit the accuracy to protect against overflow
959 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
961 # also takes care of the "error in _find_round_parameters?" case
962 return $x if $x->{sign} eq 'NaN';
964 # no rounding at all, so must use fallback
965 if (scalar @params == 0)
967 # simulate old behaviour
968 $params[0] = $self->div_scale(); # and round to it as accuracy
969 $params[1] = undef; # P = undef
970 $scale = $params[0]+4; # at least four more for proper round
971 $params[2] = $r; # round mode by caller or undef
972 $fallback = 1; # to clear a/p afterwards
976 # the 4 below is empirical, and there might be cases where it's not enough...
977 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
980 return $x->bone(@params) if $x->is_zero();
982 if (!$x->isa('Math::BigFloat'))
984 $x = Math::BigFloat->new($x);
988 # when user set globals, they would interfere with our calculation, so
989 # disable them and later re-enable them
991 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
992 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
993 # we also need to disable any set A or P on $x (_find_round_parameters took
994 # them already into account), since these would interfere, too
995 delete $x->{_a}; delete $x->{_p};
996 # need to disable $upgrade in BigInt, to avoid deep recursion
997 local $Math::BigInt::upgrade = undef;
998 local $Math::BigFloat::downgrade = undef;
1000 my $x_org = $x->copy();
1002 # We use the following Taylor series:
1005 # e = 1 + --- + --- + --- + --- ...
1008 # The difference for each term is X and N, which would result in:
1009 # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term
1011 # But it is faster to compute exp(1) and then raising it to the
1012 # given power, esp. if $x is really big and an integer because:
1014 # * The numerator is always 1, making the computation faster
1015 # * the series converges faster in the case of x == 1
1016 # * We can also easily check when we have reached our limit: when the
1017 # term to be added is smaller than "1E$scale", we can stop - f.i.
1018 # scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
1019 # * we can compute the *exact* result by simulating bigrat math:
1021 # 1 1 gcd(3,4) = 1 1*24 + 1*6 5
1022 # - + - = ---------- = --
1025 # We do not compute the gcd() here, but simple do:
1027 # - + - = --------- = --
1031 # a c a*d + c*b and note that c is always 1 and d = (b*f)
1035 # This leads to: which can be reduced by b to:
1036 # a 1 a*b*f + b a*f + 1
1037 # - + - = --------- = -------
1040 # The first terms in the series are:
1042 # 1 1 1 1 1 1 1 1 13700
1043 # -- + -- + -- + -- + -- + --- + --- + ---- = -----
1044 # 1 1 2 6 24 120 720 5040 5040
1046 # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!
1050 # set $x directly from a cached string form
1051 $x->{_m} = $MBI->_new(
1052 "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
1055 $x->{_e} = $MBI->_new(79);
1059 # compute A and B so that e = A / B.
1061 # After some terms we end up with this, so we use it as a starting point:
1062 my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
1063 my $F = $MBI->_new(42); my $step = 42;
1065 # Compute how many steps we need to take to get $A and $B sufficiently big
1066 my $steps = _len_to_steps($scale - 4);
1067 # print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
1068 while ($step++ <= $steps)
1070 # calculate $a * $f + 1
1071 $A = $MBI->_mul($A, $F);
1072 $A = $MBI->_inc($A);
1074 $F = $MBI->_inc($F);
1076 # compute $B as factorial of $steps (this is faster than doing it manually)
1077 my $B = $MBI->_fac($MBI->_new($steps));
1079 # print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
1081 # compute A/B with $scale digits in the result (truncate, not round)
1082 $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
1083 $A = $MBI->_div( $A, $B );
1088 $x->{_e} = $MBI->_new($scale);
1091 # $x contains now an estimate of e, with some surplus digits, so we can round
1092 if (!$x_org->is_one())
1094 # raise $x to the wanted power and round it in one step:
1095 $x->bpow($x_org, @params);
1099 # else just round the already computed result
1100 delete $x->{_a}; delete $x->{_p};
1101 # shortcut to not run through _find_round_parameters again
1102 if (defined $params[0])
1104 $x->bround($params[0],$params[2]); # then round accordingly
1108 $x->bfround($params[1],$params[2]); # then round accordingly
1113 # clear a/p after round, since user did not request it
1114 delete $x->{_a}; delete $x->{_p};
1117 $$abr = $ab; $$pbr = $pb;
1119 $x; # return modified $x
1124 # internal log function to calculate ln() based on Taylor series.
1125 # Modifies $x in place.
1126 my ($self,$x,$scale) = @_;
1128 # in case of $x == 1, result is 0
1129 return $x->bzero() if $x->is_one();
1131 # XXX TODO: rewrite this in a similiar manner to bexp()
1133 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
1137 # Taylor: | u 1 u^3 1 u^5 |
1138 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
1139 # |_ v 3 v^3 5 v^5 _|
1141 # This takes much more steps to calculate the result and is thus not used
1144 # Taylor: | u 1 u^2 1 u^3 |
1145 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
1146 # |_ x 2 x^2 3 x^3 _|
1148 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
1150 $v = $x->copy(); $v->binc(); # v = x+1
1151 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
1152 $x->bdiv($v,$scale); # first term: u/v
1153 $below = $v->copy();
1155 $u *= $u; $v *= $v; # u^2, v^2
1156 $below->bmul($v); # u^3, v^3
1158 $factor = $self->new(3); $f = $self->new(2);
1160 my $steps = 0 if DEBUG;
1161 $limit = $self->new("1E-". ($scale-1));
1164 # we calculate the next term, and add it to the last
1165 # when the next term is below our limit, it won't affect the outcome
1166 # anymore, so we stop
1168 # calculating the next term simple from over/below will result in quite
1169 # a time hog if the input has many digits, since over and below will
1170 # accumulate more and more digits, and the result will also have many
1171 # digits, but in the end it is rounded to $scale digits anyway. So if we
1172 # round $over and $below first, we save a lot of time for the division
1173 # (not with log(1.2345), but try log (123**123) to see what I mean. This
1174 # can introduce a rounding error if the division result would be f.i.
1175 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
1176 # if we truncated $over and $below we might get 0.12345. Does this matter
1177 # for the end result? So we give $over and $below 4 more digits to be
1178 # on the safe side (unscientific error handling as usual... :+D
1180 $next = $over->copy->bround($scale+4)->bdiv(
1181 $below->copy->bmul($factor)->bround($scale+4),
1185 ## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
1187 last if $next->bacmp($limit) <= 0;
1189 delete $next->{_a}; delete $next->{_p};
1191 # calculate things for the next term
1192 $over *= $u; $below *= $v; $factor->badd($f);
1195 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
1198 print "took $steps steps\n" if DEBUG;
1199 $x->bmul($f); # $x *= 2
1204 # Internal log function based on reducing input to the range of 0.1 .. 9.99
1205 # and then "correcting" the result to the proper one. Modifies $x in place.
1206 my ($self,$x,$scale) = @_;
1208 # Taking blog() from numbers greater than 10 takes a *very long* time, so we
1209 # break the computation down into parts based on the observation that:
1210 # blog(X*Y) = blog(X) + blog(Y)
1211 # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
1212 # $x is the faster it gets. Since 2*$x takes about 10 times as
1213 # long, we make it faster by about a factor of 100 by dividing $x by 10.
1215 # The same observation is valid for numbers smaller than 0.1, e.g. computing
1216 # log(1) is fastest, and the further away we get from 1, the longer it takes.
1217 # So we also 'break' this down by multiplying $x with 10 and subtract the
1218 # log(10) afterwards to get the correct result.
1220 # To get $x even closer to 1, we also divide by 2 and then use log(2) to
1221 # correct for this. For instance if $x is 2.4, we use the formula:
1222 # blog(2.4 * 2) == blog (1.2) + blog(2)
1223 # and thus calculate only blog(1.2) and blog(2), which is faster in total
1224 # than calculating blog(2.4).
1226 # In addition, the values for blog(2) and blog(10) are cached.
1228 # Calculate nr of digits before dot:
1229 my $dbd = $MBI->_num($x->{_e});
1230 $dbd = -$dbd if $x->{_es} eq '-';
1231 $dbd += $MBI->_len($x->{_m});
1233 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
1234 # infinite recursion
1236 my $calc = 1; # do some calculation?
1238 # disable the shortcut for 10, since we need log(10) and this would recurse
1240 if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
1242 $dbd = 0; # disable shortcut
1243 # we can use the cached value in these cases
1244 if ($scale <= $LOG_10_A)
1246 $x->bzero(); $x->badd($LOG_10); # modify $x in place
1247 $calc = 0; # no need to calc, but round
1249 # if we can't use the shortcut, we continue normally
1253 # disable the shortcut for 2, since we maybe have it cached
1254 if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
1256 $dbd = 0; # disable shortcut
1257 # we can use the cached value in these cases
1258 if ($scale <= $LOG_2_A)
1260 $x->bzero(); $x->badd($LOG_2); # modify $x in place
1261 $calc = 0; # no need to calc, but round
1263 # if we can't use the shortcut, we continue normally
1267 # if $x = 0.1, we know the result must be 0-log(10)
1268 if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
1269 $MBI->_is_one($x->{_m}))
1271 $dbd = 0; # disable shortcut
1272 # we can use the cached value in these cases
1273 if ($scale <= $LOG_10_A)
1275 $x->bzero(); $x->bsub($LOG_10);
1276 $calc = 0; # no need to calc, but round
1280 return if $calc == 0; # already have the result
1282 # default: these correction factors are undef and thus not used
1283 my $l_10; # value of ln(10) to A of $scale
1284 my $l_2; # value of ln(2) to A of $scale
1286 my $two = $self->new(2);
1288 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1289 # so don't do this shortcut for 1 or 0
1290 if (($dbd > 1) || ($dbd < 0))
1292 # convert our cached value to an object if not already (avoid doing this
1293 # at import() time, since not everybody needs this)
1294 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1296 #print "x = $x, dbd = $dbd, calc = $calc\n";
1297 # got more than one digit before the dot, or more than one zero after the
1299 # log(123) == log(1.23) + log(10) * 2
1300 # log(0.0123) == log(1.23) - log(10) * 2
1302 if ($scale <= $LOG_10_A)
1305 $l_10 = $LOG_10->copy(); # copy for mul
1309 # else: slower, compute and cache result
1310 # also disable downgrade for this code path
1311 local $Math::BigFloat::downgrade = undef;
1313 # shorten the time to calculate log(10) based on the following:
1314 # log(1.25 * 8) = log(1.25) + log(8)
1315 # = log(1.25) + log(2) + log(2) + log(2)
1317 # first get $l_2 (and possible compute and cache log(2))
1318 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1319 if ($scale <= $LOG_2_A)
1322 $l_2 = $LOG_2->copy(); # copy() for the mul below
1326 # else: slower, compute and cache result
1327 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1328 $LOG_2 = $l_2->copy(); # cache the result for later
1329 # the copy() is for mul below
1333 # now calculate log(1.25):
1334 $l_10 = $self->new('1.25'); $self->_log($l_10, $scale); # scale+4, actually
1336 # log(1.25) + log(2) + log(2) + log(2):
1340 $LOG_10 = $l_10->copy(); # cache the result for later
1341 # the copy() is for mul below
1344 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
1345 $l_10->bmul( $self->new($dbd)); # log(10) * (digits_before_dot-1)
1352 ($x->{_e}, $x->{_es}) =
1353 _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
1357 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1359 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1360 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1362 $HALF = $self->new($HALF) unless ref($HALF);
1364 my $twos = 0; # default: none (0 times)
1365 while ($x->bacmp($HALF) <= 0) # X <= 0.5
1367 $twos--; $x->bmul($two);
1369 while ($x->bacmp($two) >= 0) # X >= 2
1371 $twos++; $x->bdiv($two,$scale+4); # keep all digits
1373 # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
1374 # So calculate correction factor based on ln(2):
1377 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1378 if ($scale <= $LOG_2_A)
1381 $l_2 = $LOG_2->copy(); # copy() for the mul below
1385 # else: slower, compute and cache result
1386 # also disable downgrade for this code path
1387 local $Math::BigFloat::downgrade = undef;
1388 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1389 $LOG_2 = $l_2->copy(); # cache the result for later
1390 # the copy() is for mul below
1393 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
1396 $self->_log($x,$scale); # need to do the "normal" way
1397 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
1398 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
1400 # all done, $x contains now the result
1406 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1407 # does not modify arguments, but returns new object
1408 # Lowest Common Multiplicator
1410 my ($self,@arg) = objectify(0,@_);
1411 my $x = $self->new(shift @arg);
1412 while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); }
1418 # (BINT or num_str, BINT or num_str) return BINT
1419 # does not modify arguments, but returns new object
1422 $y = __PACKAGE__->new($y) if !ref($y);
1424 my $x = $y->copy()->babs(); # keep arguments
1426 return $x->bnan() if $x->{sign} !~ /^[+-]$/ # x NaN?
1427 || !$x->is_int(); # only for integers now
1431 my $t = shift; $t = $self->new($t) if !ref($t);
1432 $y = $t->copy()->babs();
1434 return $x->bnan() if $y->{sign} !~ /^[+-]$/ # y NaN?
1435 || !$y->is_int(); # only for integers now
1437 # greatest common divisor
1438 while (! $y->is_zero())
1440 ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1443 last if $x->is_one();
1448 ##############################################################################
1452 # Internal helper sub to take two positive integers and their signs and
1453 # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1454 # output ($CALC,('+'|'-'))
1455 my ($x,$y,$xs,$ys) = @_;
1457 # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1460 $x = $MBI->_add ($x, $y ); # a+b
1461 # the sign follows $xs
1465 my $a = $MBI->_acmp($x,$y);
1468 $x = $MBI->_sub ($x , $y); # abs sub
1472 $x = $MBI->_zero(); # result is 0
1477 $x = $MBI->_sub ( $y, $x, 1 ); # abs sub
1485 # Internal helper sub to take two positive integers and their signs and
1486 # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1487 # output ($CALC,('+'|'-'))
1488 my ($x,$y,$xs,$ys) = @_;
1492 _e_add($x,$y,$xs,$ys); # call add (does subtract now)
1495 ###############################################################################
1496 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1500 # return true if arg (BFLOAT or num_str) is an integer
1501 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1503 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1504 $x->{_es} eq '+'; # 1e-1 => no integer
1510 # return true if arg (BFLOAT or num_str) is zero
1511 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1513 return 1 if $x->{sign} eq '+' && $MBI->_is_zero($x->{_m});
1519 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1520 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1522 $sign = '+' if !defined $sign || $sign ne '-';
1524 if ($x->{sign} eq $sign &&
1525 $MBI->_is_zero($x->{_e}) && $MBI->_is_one($x->{_m}));
1531 # return true if arg (BFLOAT or num_str) is odd or false if even
1532 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1534 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1535 ($MBI->_is_zero($x->{_e}) && $MBI->_is_odd($x->{_m}));
1541 # return true if arg (BINT or num_str) is even or false if odd
1542 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1544 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
1545 return 1 if ($x->{_es} eq '+' # 123.45 is never
1546 && $MBI->_is_even($x->{_m})); # but 1200 is
1552 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
1553 # (BINT or num_str, BINT or num_str) return BINT
1556 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1557 # objectify is costly, so avoid it
1558 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1560 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1563 return $x if $x->modify('bmul');
1565 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1568 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1570 return $x->bnan() if $x->is_zero() || $y->is_zero();
1571 # result will always be +-inf:
1572 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1573 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1574 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1575 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1576 return $x->binf('-');
1579 return $x->bzero() if $x->is_zero() || $y->is_zero();
1581 return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
1582 ((!$x->isa($self)) || (!$y->isa($self)));
1584 # aEb * cEd = (a*c)E(b+d)
1585 $MBI->_mul($x->{_m},$y->{_m});
1586 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1589 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1590 return $x->bnorm()->round($a,$p,$r,$y);
1595 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1596 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1599 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1600 # objectify is costly, so avoid it
1601 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1603 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1606 return $x if $x->modify('bdiv');
1608 return $self->_div_inf($x,$y)
1609 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1611 # x== 0 # also: or y == 1 or y == -1
1612 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1615 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1617 # we need to limit the accuracy to protect against overflow
1619 my (@params,$scale);
1620 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1622 return $x if $x->is_nan(); # error in _find_round_parameters?
1624 # no rounding at all, so must use fallback
1625 if (scalar @params == 0)
1627 # simulate old behaviour
1628 $params[0] = $self->div_scale(); # and round to it as accuracy
1629 $scale = $params[0]+4; # at least four more for proper round
1630 $params[2] = $r; # round mode by caller or undef
1631 $fallback = 1; # to clear a/p afterwards
1635 # the 4 below is empirical, and there might be cases where it is not
1637 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1640 my $rem; $rem = $self->bzero() if wantarray;
1642 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1644 my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
1645 $scale = $lx if $lx > $scale;
1646 $scale = $ly if $ly > $scale;
1647 my $diff = $ly - $lx;
1648 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
1650 # already handled inf/NaN/-inf above:
1652 # check that $y is not 1 nor -1 and cache the result:
1653 my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));
1655 # flipping the sign of $y will also flip the sign of $x for the special
1656 # case of $x->bsub($x); so we can catch it below:
1657 my $xsign = $x->{sign};
1658 $y->{sign} =~ tr/+-/-+/;
1660 if ($xsign ne $x->{sign})
1662 # special case of $x /= $x results in 1
1663 $x->bone(); # "fixes" also sign of $y, since $x is $y
1667 # correct $y's sign again
1668 $y->{sign} =~ tr/+-/-+/;
1669 # continue with normal div code:
1671 # make copy of $x in case of list context for later reminder calculation
1672 if (wantarray && $y_not_one)
1677 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1679 # check for / +-1 ( +/- 1E0)
1682 # promote BigInts and it's subclasses (except when already a BigFloat)
1683 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1685 # calculate the result to $scale digits and then round it
1686 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1687 $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1688 $MBI->_div ($x->{_m},$y->{_m}); # a/c
1690 # correct exponent of $x
1691 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1692 # correct for 10**scale
1693 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1694 $x->bnorm(); # remove trailing 0's
1696 } # ende else $x != $y
1698 # shortcut to not run through _find_round_parameters again
1699 if (defined $params[0])
1701 delete $x->{_a}; # clear before round
1702 $x->bround($params[0],$params[2]); # then round accordingly
1706 delete $x->{_p}; # clear before round
1707 $x->bfround($params[1],$params[2]); # then round accordingly
1711 # clear a/p after round, since user did not request it
1712 delete $x->{_a}; delete $x->{_p};
1719 $rem->bmod($y,@params); # copy already done
1723 # clear a/p after round, since user did not request it
1724 delete $rem->{_a}; delete $rem->{_p};
1733 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1736 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1737 # objectify is costly, so avoid it
1738 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1740 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1743 return $x if $x->modify('bmod');
1745 # handle NaN, inf, -inf
1746 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1748 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1749 $x->{sign} = $re->{sign};
1750 $x->{_e} = $re->{_e};
1751 $x->{_m} = $re->{_m};
1752 return $x->round($a,$p,$r,$y);
1756 return $x->bnan() if $x->is_zero();
1760 return $x->bzero() if $x->is_zero()
1762 # check that $y == +1 or $y == -1:
1763 ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));
1765 my $cmp = $x->bacmp($y); # equal or $x < $y?
1766 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1768 # only $y of the operands negative?
1769 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1771 $x->{sign} = $y->{sign}; # calc sign first
1772 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1774 my $ym = $MBI->_copy($y->{_m});
1777 $MBI->_lsft( $ym, $y->{_e}, 10)
1778 if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1780 # if $y has digits after dot
1781 my $shifty = 0; # correct _e of $x by this
1782 if ($y->{_es} eq '-') # has digits after dot
1784 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1785 $shifty = $MBI->_num($y->{_e}); # no more digits after dot
1786 $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1788 # $ym is now mantissa of $y based on exponent 0
1790 my $shiftx = 0; # correct _e of $x by this
1791 if ($x->{_es} eq '-') # has digits after dot
1793 # 123.4 % 20 => 1234 % 200
1794 $shiftx = $MBI->_num($x->{_e}); # no more digits after dot
1795 $MBI->_lsft($ym, $x->{_e}, 10); # 123 => 1230
1797 # 123e1 % 20 => 1230 % 20
1798 if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1800 $MBI->_lsft( $x->{_m}, $x->{_e},10); # es => '+' here
1803 $x->{_e} = $MBI->_new($shiftx);
1805 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1806 $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1808 # now mantissas are equalized, exponent of $x is adjusted, so calc result
1810 $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1812 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1815 if ($neg != 0) # one of them negative => correct in place
1818 $x->{_m} = $r->{_m};
1819 $x->{_e} = $r->{_e};
1820 $x->{_es} = $r->{_es};
1821 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1825 $x->round($a,$p,$r,$y); # round and return
1830 # calculate $y'th root of $x
1833 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1834 # objectify is costly, so avoid it
1835 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1837 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1840 return $x if $x->modify('broot');
1842 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1843 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1844 $y->{sign} !~ /^\+$/;
1846 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1848 # we need to limit the accuracy to protect against overflow
1850 my (@params,$scale);
1851 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1853 return $x if $x->is_nan(); # error in _find_round_parameters?
1855 # no rounding at all, so must use fallback
1856 if (scalar @params == 0)
1858 # simulate old behaviour
1859 $params[0] = $self->div_scale(); # and round to it as accuracy
1860 $scale = $params[0]+4; # at least four more for proper round
1861 $params[2] = $r; # iound mode by caller or undef
1862 $fallback = 1; # to clear a/p afterwards
1866 # the 4 below is empirical, and there might be cases where it is not
1868 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1871 # when user set globals, they would interfere with our calculation, so
1872 # disable them and later re-enable them
1874 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1875 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1876 # we also need to disable any set A or P on $x (_find_round_parameters took
1877 # them already into account), since these would interfere, too
1878 delete $x->{_a}; delete $x->{_p};
1879 # need to disable $upgrade in BigInt, to avoid deep recursion
1880 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1882 # remember sign and make $x positive, since -4 ** (1/2) => -2
1883 my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
1886 if ($y->isa('Math::BigFloat'))
1888 $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
1892 $is_two = ($y == 2);
1895 # normal square root if $y == 2:
1898 $x->bsqrt($scale+4);
1900 elsif ($y->is_one('-'))
1903 my $u = $self->bone()->bdiv($x,$scale);
1904 # copy private parts over
1905 $x->{_m} = $u->{_m};
1906 $x->{_e} = $u->{_e};
1907 $x->{_es} = $u->{_es};
1911 # calculate the broot() as integer result first, and if it fits, return
1912 # it rightaway (but only if $x and $y are integer):
1914 my $done = 0; # not yet
1915 if ($y->is_int() && $x->is_int())
1917 my $i = $MBI->_copy( $x->{_m} );
1918 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1919 my $int = Math::BigInt->bzero();
1921 $int->broot($y->as_number());
1923 if ($int->copy()->bpow($y) == $x)
1925 # found result, return it
1926 $x->{_m} = $int->{value};
1927 $x->{_e} = $MBI->_zero();
1935 my $u = $self->bone()->bdiv($y,$scale+4);
1936 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
1937 $x->bpow($u,$scale+4); # el cheapo
1940 $x->bneg() if $sign == 1;
1942 # shortcut to not run through _find_round_parameters again
1943 if (defined $params[0])
1945 $x->bround($params[0],$params[2]); # then round accordingly
1949 $x->bfround($params[1],$params[2]); # then round accordingly
1953 # clear a/p after round, since user did not request it
1954 delete $x->{_a}; delete $x->{_p};
1957 $$abr = $ab; $$pbr = $pb;
1963 # calculate square root
1964 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1966 return $x if $x->modify('bsqrt');
1968 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
1969 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
1970 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
1972 # we need to limit the accuracy to protect against overflow
1974 my (@params,$scale);
1975 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1977 return $x if $x->is_nan(); # error in _find_round_parameters?
1979 # no rounding at all, so must use fallback
1980 if (scalar @params == 0)
1982 # simulate old behaviour
1983 $params[0] = $self->div_scale(); # and round to it as accuracy
1984 $scale = $params[0]+4; # at least four more for proper round
1985 $params[2] = $r; # round mode by caller or undef
1986 $fallback = 1; # to clear a/p afterwards
1990 # the 4 below is empirical, and there might be cases where it is not
1992 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1995 # when user set globals, they would interfere with our calculation, so
1996 # disable them and later re-enable them
1998 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1999 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2000 # we also need to disable any set A or P on $x (_find_round_parameters took
2001 # them already into account), since these would interfere, too
2002 delete $x->{_a}; delete $x->{_p};
2003 # need to disable $upgrade in BigInt, to avoid deep recursion
2004 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2006 my $i = $MBI->_copy( $x->{_m} );
2007 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2008 my $xas = Math::BigInt->bzero();
2011 my $gs = $xas->copy()->bsqrt(); # some guess
2013 if (($x->{_es} ne '-') # guess can't be accurate if there are
2014 # digits after the dot
2015 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
2017 # exact result, copy result over to keep $x
2018 $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
2020 # shortcut to not run through _find_round_parameters again
2021 if (defined $params[0])
2023 $x->bround($params[0],$params[2]); # then round accordingly
2027 $x->bfround($params[1],$params[2]); # then round accordingly
2031 # clear a/p after round, since user did not request it
2032 delete $x->{_a}; delete $x->{_p};
2034 # re-enable A and P, upgrade is taken care of by "local"
2035 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
2039 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
2040 # of the result by multipyling the input by 100 and then divide the integer
2041 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
2043 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
2044 my $y1 = $MBI->_copy($x->{_m});
2046 my $length = $MBI->_len($y1);
2048 # Now calculate how many digits the result of sqrt(y1) would have
2049 my $digits = int($length / 2);
2051 # But we need at least $scale digits, so calculate how many are missing
2052 my $shift = $scale - $digits;
2054 # That should never happen (we take care of integer guesses above)
2055 # $shift = 0 if $shift < 0;
2057 # Multiply in steps of 100, by shifting left two times the "missing" digits
2058 my $s2 = $shift * 2;
2060 # We now make sure that $y1 has the same odd or even number of digits than
2061 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
2062 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
2063 # steps of 10. The length of $x does not count, since an even or odd number
2064 # of digits before the dot is not changed by adding an even number of digits
2065 # after the dot (the result is still odd or even digits long).
2066 $s2++ if $MBI->_is_odd($x->{_e});
2068 $MBI->_lsft( $y1, $MBI->_new($s2), 10);
2070 # now take the square root and truncate to integer
2071 $y1 = $MBI->_sqrt($y1);
2073 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
2074 # result, which is than later rounded to the desired scale.
2076 # calculate how many zeros $x had after the '.' (or before it, depending
2077 # on sign of $dat, the result should have half as many:
2078 my $dat = $MBI->_num($x->{_e});
2079 $dat = -$dat if $x->{_es} eq '-';
2084 # no zeros after the dot (e.g. 1.23, 0.49 etc)
2085 # preserve half as many digits before the dot than the input had
2086 # (but round this "up")
2087 $dat = int(($dat+1)/2);
2091 $dat = int(($dat)/2);
2093 $dat -= $MBI->_len($y1);
2097 $x->{_e} = $MBI->_new( $dat );
2102 $x->{_e} = $MBI->_new( $dat );
2108 # shortcut to not run through _find_round_parameters again
2109 if (defined $params[0])
2111 $x->bround($params[0],$params[2]); # then round accordingly
2115 $x->bfround($params[1],$params[2]); # then round accordingly
2119 # clear a/p after round, since user did not request it
2120 delete $x->{_a}; delete $x->{_p};
2123 $$abr = $ab; $$pbr = $pb;
2129 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2130 # compute factorial number, modifies first argument
2133 my ($self,$x,@r) = (ref($_[0]),@_);
2134 # objectify is costly, so avoid it
2135 ($self,$x,@r) = objectify(1,@_) if !ref($x);
2138 return $x if $x->modify('bfac') || $x->{sign} eq '+inf';
2141 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
2142 ($x->{_es} ne '+')); # digits after dot?
2144 # use BigInt's bfac() for faster calc
2145 if (! $MBI->_is_zero($x->{_e}))
2147 $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
2148 $x->{_e} = $MBI->_zero(); # normalize
2151 $MBI->_fac($x->{_m}); # calculate factorial
2152 $x->bnorm()->round(@r); # norm again and round result
2157 # Calculate a power where $y is a non-integer, like 2 ** 0.5
2158 my ($x,$y,$a,$p,$r) = @_;
2161 # if $y == 0.5, it is sqrt($x)
2162 $HALF = $self->new($HALF) unless ref($HALF);
2163 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp($HALF) == 0;
2166 # a ** x == e ** (x * ln a)
2170 # Taylor: | u u^2 u^3 |
2171 # x ** y = 1 + | --- + --- + ----- + ... |
2174 # we need to limit the accuracy to protect against overflow
2176 my ($scale,@params);
2177 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2179 return $x if $x->is_nan(); # error in _find_round_parameters?
2181 # no rounding at all, so must use fallback
2182 if (scalar @params == 0)
2184 # simulate old behaviour
2185 $params[0] = $self->div_scale(); # and round to it as accuracy
2186 $params[1] = undef; # disable P
2187 $scale = $params[0]+4; # at least four more for proper round
2188 $params[2] = $r; # round mode by caller or undef
2189 $fallback = 1; # to clear a/p afterwards
2193 # the 4 below is empirical, and there might be cases where it is not
2195 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2198 # when user set globals, they would interfere with our calculation, so
2199 # disable them and later re-enable them
2201 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2202 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2203 # we also need to disable any set A or P on $x (_find_round_parameters took
2204 # them already into account), since these would interfere, too
2205 delete $x->{_a}; delete $x->{_p};
2206 # need to disable $upgrade in BigInt, to avoid deep recursion
2207 local $Math::BigInt::upgrade = undef;
2209 my ($limit,$v,$u,$below,$factor,$next,$over);
2211 $u = $x->copy()->blog(undef,$scale)->bmul($y);
2212 $v = $self->bone(); # 1
2213 $factor = $self->new(2); # 2
2214 $x->bone(); # first term: 1
2216 $below = $v->copy();
2219 $limit = $self->new("1E-". ($scale-1));
2223 # we calculate the next term, and add it to the last
2224 # when the next term is below our limit, it won't affect the outcome
2225 # anymore, so we stop:
2226 $next = $over->copy()->bdiv($below,$scale);
2227 last if $next->bacmp($limit) <= 0;
2229 # calculate things for the next term
2230 $over *= $u; $below *= $factor; $factor->binc();
2232 last if $x->{sign} !~ /^[-+]$/;
2237 # shortcut to not run through _find_round_parameters again
2238 if (defined $params[0])
2240 $x->bround($params[0],$params[2]); # then round accordingly
2244 $x->bfround($params[1],$params[2]); # then round accordingly
2248 # clear a/p after round, since user did not request it
2249 delete $x->{_a}; delete $x->{_p};
2252 $$abr = $ab; $$pbr = $pb;
2258 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2259 # compute power of two numbers, second arg is used as integer
2260 # modifies first argument
2263 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2264 # objectify is costly, so avoid it
2265 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2267 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2270 return $x if $x->modify('bpow');
2272 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
2273 return $x if $x->{sign} =~ /^[+-]inf$/;
2275 # cache the result of is_zero
2276 my $y_is_zero = $y->is_zero();
2277 return $x->bone() if $y_is_zero;
2278 return $x if $x->is_one() || $y->is_one();
2280 my $x_is_zero = $x->is_zero();
2281 return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int(); # non-integer power
2283 my $y1 = $y->as_number()->{value}; # make MBI part
2286 if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
2288 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
2289 return $MBI->_is_odd($y1) ? $x : $x->babs(1);
2293 return $x->bone() if $y_is_zero;
2294 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
2295 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
2300 $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
2302 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
2303 $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
2304 $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
2306 $x->{sign} = $new_sign;
2308 if ($y->{sign} eq '-')
2310 # modify $x in place!
2311 my $z = $x->copy(); $x->bone();
2312 return scalar $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
2314 $x->round($a,$p,$r,$y);
2317 ###############################################################################
2318 # rounding functions
2322 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
2323 # $n == 0 means round to integer
2324 # expects and returns normalized numbers!
2325 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2327 my ($scale,$mode) = $x->_scale_p(@_);
2328 return $x if !defined $scale || $x->modify('bfround'); # no-op
2330 # never round a 0, +-inf, NaN
2333 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
2336 return $x if $x->{sign} !~ /^[+-]$/;
2338 # don't round if x already has lower precision
2339 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
2341 $x->{_p} = $scale; # remember round in any case
2342 delete $x->{_a}; # and clear A
2345 # round right from the '.'
2347 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round
2349 $scale = -$scale; # positive for simplicity
2350 my $len = $MBI->_len($x->{_m}); # length of mantissa
2352 # the following poses a restriction on _e, but if _e is bigger than a
2353 # scalar, you got other problems (memory etc) anyway
2354 my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e}))); # digits after dot
2355 my $zad = 0; # zeros after dot
2356 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
2358 # p rint "scale $scale dad $dad zad $zad len $len\n";
2359 # number bsstr len zad dad
2360 # 0.123 123e-3 3 0 3
2361 # 0.0123 123e-4 3 1 4
2364 # 1.2345 12345e-4 5 0 4
2366 # do not round after/right of the $dad
2367 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
2369 # round to zero if rounding inside the $zad, but not for last zero like:
2370 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
2371 return $x->bzero() if $scale < $zad;
2372 if ($scale == $zad) # for 0.006, scale -3 and trunc
2378 # adjust round-point to be inside mantissa
2381 $scale = $scale-$zad;
2385 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
2386 $scale = $dbd+$scale;
2392 # round left from the '.'
2394 # 123 => 100 means length(123) = 3 - $scale (2) => 1
2396 my $dbt = $MBI->_len($x->{_m});
2398 my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
2399 # should be the same, so treat it as this
2400 $scale = 1 if $scale == 0;
2401 # shortcut if already integer
2402 return $x if $scale == 1 && $dbt <= $dbd;
2403 # maximum digits before dot
2408 # not enough digits before dot, so round to zero
2411 elsif ( $scale == $dbd )
2418 $scale = $dbd - $scale;
2421 # pass sign to bround for rounding modes '+inf' and '-inf'
2422 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2423 $m->bround($scale,$mode);
2424 $x->{_m} = $m->{value}; # get our mantissa back
2430 # accuracy: preserve $N digits, and overwrite the rest with 0's
2431 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2433 if (($_[0] || 0) < 0)
2435 require Carp; Carp::croak ('bround() needs positive accuracy');
2438 my ($scale,$mode) = $x->_scale_a(@_);
2439 return $x if !defined $scale || $x->modify('bround'); # no-op
2441 # scale is now either $x->{_a}, $accuracy, or the user parameter
2442 # test whether $x already has lower accuracy, do nothing in this case
2443 # but do round if the accuracy is the same, since a math operation might
2444 # want to round a number with A=5 to 5 digits afterwards again
2445 return $x if defined $x->{_a} && $x->{_a} < $scale;
2447 # scale < 0 makes no sense
2448 # scale == 0 => keep all digits
2449 # never round a +-inf, NaN
2450 return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
2452 # 1: never round a 0
2453 # 2: if we should keep more digits than the mantissa has, do nothing
2454 if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
2456 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
2460 # pass sign to bround for '+inf' and '-inf' rounding modes
2461 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2463 $m->bround($scale,$mode); # round mantissa
2464 $x->{_m} = $m->{value}; # get our mantissa back
2465 $x->{_a} = $scale; # remember rounding
2466 delete $x->{_p}; # and clear P
2467 $x->bnorm(); # del trailing zeros gen. by bround()
2472 # return integer less or equal then $x
2473 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2475 return $x if $x->modify('bfloor');
2477 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2479 # if $x has digits after dot
2480 if ($x->{_es} eq '-')
2482 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2483 $x->{_e} = $MBI->_zero(); # trunc/norm
2484 $x->{_es} = '+'; # abs e
2485 $MBI->_inc($x->{_m}) if $x->{sign} eq '-'; # increment if negative
2487 $x->round($a,$p,$r);
2492 # return integer greater or equal then $x
2493 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2495 return $x if $x->modify('bceil');
2496 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2498 # if $x has digits after dot
2499 if ($x->{_es} eq '-')
2501 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2502 $x->{_e} = $MBI->_zero(); # trunc/norm
2503 $x->{_es} = '+'; # abs e
2504 $MBI->_inc($x->{_m}) if $x->{sign} eq '+'; # increment if positive
2506 $x->round($a,$p,$r);
2511 # shift right by $y (divide by power of $n)
2514 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2515 # objectify is costly, so avoid it
2516 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2518 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2521 return $x if $x->modify('brsft');
2522 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2524 $n = 2 if !defined $n; $n = $self->new($n);
2527 return $x->blsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
2529 # the following call to bdiv() will return either quo or (quo,reminder):
2530 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
2535 # shift left by $y (multiply by power of $n)
2538 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2539 # objectify is costly, so avoid it
2540 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2542 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2545 return $x if $x->modify('blsft');
2546 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2548 $n = 2 if !defined $n; $n = $self->new($n);
2551 return $x->brsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
2553 $x->bmul($n->bpow($y),$a,$p,$r,$y);
2556 ###############################################################################
2560 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
2565 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
2566 # or falling back to MBI::bxxx()
2567 my $name = $AUTOLOAD;
2569 $name =~ s/(.*):://; # split package
2570 my $c = $1 || $class;
2572 $c->import() if $IMPORT == 0;
2573 if (!_method_alias($name))
2577 # delayed load of Carp and avoid recursion
2579 Carp::croak ("$c: Can't call a method without name");
2581 if (!_method_hand_up($name))
2583 # delayed load of Carp and avoid recursion
2585 Carp::croak ("Can't call $c\-\>$name, not a valid method");
2587 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
2589 return &{"Math::BigInt"."::$name"}(@_);
2591 my $bname = $name; $bname =~ s/^f/b/;
2599 # return a copy of the exponent
2600 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2602 if ($x->{sign} !~ /^[+-]$/)
2604 my $s = $x->{sign}; $s =~ s/^[+-]//;
2605 return Math::BigInt->new($s); # -inf, +inf => +inf
2607 Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
2612 # return a copy of the mantissa
2613 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2615 if ($x->{sign} !~ /^[+-]$/)
2617 my $s = $x->{sign}; $s =~ s/^[+]//;
2618 return Math::BigInt->new($s); # -inf, +inf => +inf
2620 my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
2621 $m->bneg() if $x->{sign} eq '-';
2628 # return a copy of both the exponent and the mantissa
2629 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2631 if ($x->{sign} !~ /^[+-]$/)
2633 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2634 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2636 my $m = Math::BigInt->bzero();
2637 $m->{value} = $MBI->_copy($x->{_m});
2638 $m->bneg() if $x->{sign} eq '-';
2639 ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
2642 ##############################################################################
2643 # private stuff (internal use only)
2649 my $lib = ''; my @a;
2650 my $lib_kind = 'try';
2652 for ( my $i = 0; $i < $l ; $i++)
2654 if ( $_[$i] eq ':constant' )
2656 # This causes overlord er load to step in. 'binary' and 'integer'
2657 # are handled by BigInt.
2658 overload::constant float => sub { $self->new(shift); };
2660 elsif ($_[$i] eq 'upgrade')
2662 # this causes upgrading
2663 $upgrade = $_[$i+1]; # or undef to disable
2666 elsif ($_[$i] eq 'downgrade')
2668 # this causes downgrading
2669 $downgrade = $_[$i+1]; # or undef to disable
2672 elsif ($_[$i] =~ /^(lib|try|only)\z/)
2674 # alternative library
2675 $lib = $_[$i+1] || ''; # default Calc
2676 $lib_kind = $1; # lib, try or only
2679 elsif ($_[$i] eq 'with')
2681 # alternative class for our private parts()
2682 # XXX: no longer supported
2683 # $MBI = $_[$i+1] || 'Math::BigInt';
2692 $lib =~ tr/a-zA-Z0-9,://cd; # restrict to sane characters
2693 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2694 my $mbilib = eval { Math::BigInt->config()->{lib} };
2695 if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
2697 # MBI already loaded
2698 Math::BigInt->import( $lib_kind, "$lib,$mbilib", 'objectify');
2702 # MBI not loaded, or with ne "Math::BigInt::Calc"
2703 $lib .= ",$mbilib" if defined $mbilib;
2704 $lib =~ s/^,//; # don't leave empty
2706 # replacement library can handle lib statement, but also could ignore it
2708 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2709 # used in the same script, or eval inside import(). So we require MBI:
2710 require Math::BigInt;
2711 Math::BigInt->import( $lib_kind => $lib, 'objectify' );
2715 require Carp; Carp::croak ("Couldn't load $lib: $! $@");
2717 # find out which one was actually loaded
2718 $MBI = Math::BigInt->config()->{lib};
2720 # register us with MBI to get notified of future lib changes
2721 Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
2723 # any non :constant stuff is handled by our parent, Exporter
2724 # even if @_ is empty, to give it a chance
2725 $self->SUPER::import(@a); # for subclasses
2726 $self->export_to_level(1,$self,@a); # need this, too
2731 # adjust m and e so that m is smallest possible
2732 # round number according to accuracy and precision settings
2733 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
2735 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2737 my $zeros = $MBI->_zeros($x->{_m}); # correct for trailing zeros
2740 my $z = $MBI->_new($zeros);
2741 $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
2742 if ($x->{_es} eq '-')
2744 if ($MBI->_acmp($x->{_e},$z) >= 0)
2746 $x->{_e} = $MBI->_sub ($x->{_e}, $z);
2747 $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
2751 $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
2757 $x->{_e} = $MBI->_add ($x->{_e}, $z);
2762 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
2763 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
2764 $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
2765 if $MBI->_is_zero($x->{_m});
2768 $x; # MBI bnorm is no-op, so dont call it
2771 ##############################################################################
2775 # return number as hexadecimal string (only for integers defined)
2776 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2778 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2779 return '0x0' if $x->is_zero();
2781 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2783 my $z = $MBI->_copy($x->{_m});
2784 if (! $MBI->_is_zero($x->{_e})) # > 0
2786 $MBI->_lsft( $z, $x->{_e},10);
2788 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2794 # return number as binary digit string (only for integers defined)
2795 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2797 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2798 return '0b0' if $x->is_zero();
2800 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2802 my $z = $MBI->_copy($x->{_m});
2803 if (! $MBI->_is_zero($x->{_e})) # > 0
2805 $MBI->_lsft( $z, $x->{_e},10);
2807 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2813 # return number as octal digit string (only for integers defined)
2814 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2816 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2817 return '0' if $x->is_zero();
2819 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2821 my $z = $MBI->_copy($x->{_m});
2822 if (! $MBI->_is_zero($x->{_e})) # > 0
2824 $MBI->_lsft( $z, $x->{_e},10);
2826 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2832 # return copy as a bigint representation of this BigFloat number
2833 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2835 return $x if $x->modify('as_number');
2837 my $z = $MBI->_copy($x->{_m});
2838 if ($x->{_es} eq '-') # < 0
2840 $MBI->_rsft( $z, $x->{_e},10);
2842 elsif (! $MBI->_is_zero($x->{_e})) # > 0
2844 $MBI->_lsft( $z, $x->{_e},10);
2846 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2853 my $class = ref($x) || $x;
2854 $x = $class->new(shift) unless ref($x);
2856 return 1 if $MBI->_is_zero($x->{_m});
2858 my $len = $MBI->_len($x->{_m});
2859 $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
2863 $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
2874 Math::BigFloat - Arbitrary size floating point math package
2881 $x = Math::BigFloat->new($str); # defaults to 0
2882 $nan = Math::BigFloat->bnan(); # create a NotANumber
2883 $zero = Math::BigFloat->bzero(); # create a +0
2884 $inf = Math::BigFloat->binf(); # create a +inf
2885 $inf = Math::BigFloat->binf('-'); # create a -inf
2886 $one = Math::BigFloat->bone(); # create a +1
2887 $one = Math::BigFloat->bone('-'); # create a -1
2890 $x->is_zero(); # true if arg is +0
2891 $x->is_nan(); # true if arg is NaN
2892 $x->is_one(); # true if arg is +1
2893 $x->is_one('-'); # true if arg is -1
2894 $x->is_odd(); # true if odd, false for even
2895 $x->is_even(); # true if even, false for odd
2896 $x->is_pos(); # true if >= 0
2897 $x->is_neg(); # true if < 0
2898 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2900 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2901 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2902 $x->sign(); # return the sign, either +,- or NaN
2903 $x->digit($n); # return the nth digit, counting from right
2904 $x->digit(-$n); # return the nth digit, counting from left
2906 # The following all modify their first argument. If you want to preserve
2907 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2908 # necessary when mixing $a = $b assignments with non-overloaded math.
2911 $x->bzero(); # set $i to 0
2912 $x->bnan(); # set $i to NaN
2913 $x->bone(); # set $x to +1
2914 $x->bone('-'); # set $x to -1
2915 $x->binf(); # set $x to inf
2916 $x->binf('-'); # set $x to -inf
2918 $x->bneg(); # negation
2919 $x->babs(); # absolute value
2920 $x->bnorm(); # normalize (no-op)
2921 $x->bnot(); # two's complement (bit wise not)
2922 $x->binc(); # increment x by 1
2923 $x->bdec(); # decrement x by 1
2925 $x->badd($y); # addition (add $y to $x)
2926 $x->bsub($y); # subtraction (subtract $y from $x)
2927 $x->bmul($y); # multiplication (multiply $x by $y)
2928 $x->bdiv($y); # divide, set $x to quotient
2929 # return (quo,rem) or quo if scalar
2931 $x->bmod($y); # modulus ($x % $y)
2932 $x->bpow($y); # power of arguments ($x ** $y)
2933 $x->blsft($y, $n); # left shift by $y places in base $n
2934 $x->brsft($y, $n); # right shift by $y places in base $n
2935 # returns (quo,rem) or quo if in scalar context
2937 $x->blog(); # logarithm of $x to base e (Euler's number)
2938 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
2939 $x->bexp(); # calculate e ** $x where e is Euler's number
2941 $x->band($y); # bit-wise and
2942 $x->bior($y); # bit-wise inclusive or
2943 $x->bxor($y); # bit-wise exclusive or
2944 $x->bnot(); # bit-wise not (two's complement)
2946 $x->bsqrt(); # calculate square-root
2947 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
2948 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2950 $x->bround($N); # accuracy: preserve $N digits
2951 $x->bfround($N); # precision: round to the $Nth digit
2953 $x->bfloor(); # return integer less or equal than $x
2954 $x->bceil(); # return integer greater or equal than $x
2956 # The following do not modify their arguments:
2958 bgcd(@values); # greatest common divisor
2959 blcm(@values); # lowest common multiplicator
2961 $x->bstr(); # return string
2962 $x->bsstr(); # return string in scientific notation
2964 $x->as_int(); # return $x as BigInt
2965 $x->exponent(); # return exponent as BigInt
2966 $x->mantissa(); # return mantissa as BigInt
2967 $x->parts(); # return (mantissa,exponent) as BigInt
2969 $x->length(); # number of digits (w/o sign and '.')
2970 ($l,$f) = $x->length(); # number of digits, and length of fraction
2972 $x->precision(); # return P of $x (or global, if P of $x undef)
2973 $x->precision($n); # set P of $x to $n
2974 $x->accuracy(); # return A of $x (or global, if A of $x undef)
2975 $x->accuracy($n); # set A $x to $n
2977 # these get/set the appropriate global value for all BigFloat objects
2978 Math::BigFloat->precision(); # Precision
2979 Math::BigFloat->accuracy(); # Accuracy
2980 Math::BigFloat->round_mode(); # rounding mode
2984 All operators (including basic math operations) are overloaded if you
2985 declare your big floating point numbers as
2987 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2989 Operations with overloaded operators preserve the arguments, which is
2990 exactly what you expect.
2992 =head2 Canonical notation
2994 Input to these routines are either BigFloat objects, or strings of the
2995 following four forms:
3009 C</^[+-]\d+E[+-]?\d+$/>
3013 C</^[+-]\d*\.\d+E[+-]?\d+$/>
3017 all with optional leading and trailing zeros and/or spaces. Additionally,
3018 numbers are allowed to have an underscore between any two digits.
3020 Empty strings as well as other illegal numbers results in 'NaN'.
3022 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
3023 are always stored in normalized form. On a string, it creates a BigFloat
3028 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
3030 The string output will always have leading and trailing zeros stripped and drop
3031 a plus sign. C<bstr()> will give you always the form with a decimal point,
3032 while C<bsstr()> (s for scientific) gives you the scientific notation.
3034 Input bstr() bsstr()
3036 ' -123 123 123' '-123123123' '-123123123E0'
3037 '00.0123' '0.0123' '123E-4'
3038 '123.45E-2' '1.2345' '12345E-4'
3039 '10E+3' '10000' '1E4'
3041 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
3042 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
3043 return either undef, <0, 0 or >0 and are suited for sort.
3045 Actual math is done by using the class defined with C<with => Class;> (which
3046 defaults to BigInts) to represent the mantissa and exponent.
3048 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
3049 represent the result when input arguments are not numbers, as well as
3050 the result of dividing by zero.
3052 =head2 C<mantissa()>, C<exponent()> and C<parts()>
3054 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
3055 as BigInts such that:
3057 $m = $x->mantissa();
3058 $e = $x->exponent();
3059 $y = $m * ( 10 ** $e );
3060 print "ok\n" if $x == $y;
3062 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
3064 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
3066 Currently the mantissa is reduced as much as possible, favouring higher
3067 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
3068 This might change in the future, so do not depend on it.
3070 =head2 Accuracy vs. Precision
3072 See also: L<Rounding|Rounding>.
3074 Math::BigFloat supports both precision (rounding to a certain place before or
3075 after the dot) and accuracy (rounding to a certain number of digits). For a
3076 full documentation, examples and tips on these topics please see the large
3077 section about rounding in L<Math::BigInt>.
3079 Since things like C<sqrt(2)> or C<1 / 3> must presented with a limited
3080 accuracy lest a operation consumes all resources, each operation produces
3081 no more than the requested number of digits.
3083 If there is no gloabl precision or accuracy set, B<and> the operation in
3084 question was not called with a requested precision or accuracy, B<and> the
3085 input $x has no accuracy or precision set, then a fallback parameter will
3086 be used. For historical reasons, it is called C<div_scale> and can be accessed
3089 $d = Math::BigFloat->div_scale(); # query
3090 Math::BigFloat->div_scale($n); # set to $n digits
3092 The default value for C<div_scale> is 40.
3094 In case the result of one operation has more digits than specified,
3095 it is rounded. The rounding mode taken is either the default mode, or the one
3096 supplied to the operation after the I<scale>:
3098 $x = Math::BigFloat->new(2);
3099 Math::BigFloat->accuracy(5); # 5 digits max
3100 $y = $x->copy()->bdiv(3); # will give 0.66667
3101 $y = $x->copy()->bdiv(3,6); # will give 0.666667
3102 $y = $x->copy()->bdiv(3,6,undef,'odd'); # will give 0.666667
3103 Math::BigFloat->round_mode('zero');
3104 $y = $x->copy()->bdiv(3,6); # will also give 0.666667
3106 Note that C<< Math::BigFloat->accuracy() >> and C<< Math::BigFloat->precision() >>
3107 set the global variables, and thus B<any> newly created number will be subject
3108 to the global rounding B<immediately>. This means that in the examples above, the
3109 C<3> as argument to C<bdiv()> will also get an accuracy of B<5>.
3111 It is less confusing to either calculate the result fully, and afterwards
3112 round it explicitly, or use the additional parameters to the math
3116 $x = Math::BigFloat->new(2);
3117 $y = $x->copy()->bdiv(3);
3118 print $y->bround(5),"\n"; # will give 0.66667
3123 $x = Math::BigFloat->new(2);
3124 $y = $x->copy()->bdiv(3,5); # will give 0.66667
3131 =item ffround ( +$scale )
3133 Rounds to the $scale'th place left from the '.', counting from the dot.
3134 The first digit is numbered 1.
3136 =item ffround ( -$scale )
3138 Rounds to the $scale'th place right from the '.', counting from the dot.
3142 Rounds to an integer.
3144 =item fround ( +$scale )
3146 Preserves accuracy to $scale digits from the left (aka significant digits)
3147 and pads the rest with zeros. If the number is between 1 and -1, the
3148 significant digits count from the first non-zero after the '.'
3150 =item fround ( -$scale ) and fround ( 0 )
3152 These are effectively no-ops.
3156 All rounding functions take as a second parameter a rounding mode from one of
3157 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
3159 The default rounding mode is 'even'. By using
3160 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
3161 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
3162 no longer supported.
3163 The second parameter to the round functions then overrides the default
3166 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
3167 'trunc' as rounding mode to make it equivalent to:
3172 You can override this by passing the desired rounding mode as parameter to
3175 $x = Math::BigFloat->new(2.5);
3176 $y = $x->as_number('odd'); # $y = 3
3182 $x->accuracy(5); # local for $x
3183 CLASS->accuracy(5); # global for all members of CLASS
3184 # Note: This also applies to new()!
3186 $A = $x->accuracy(); # read out accuracy that affects $x
3187 $A = CLASS->accuracy(); # read out global accuracy
3189 Set or get the global or local accuracy, aka how many significant digits the
3190 results have. If you set a global accuracy, then this also applies to new()!
3192 Warning! The accuracy I<sticks>, e.g. once you created a number under the
3193 influence of C<< CLASS->accuracy($A) >>, all results from math operations with
3194 that number will also be rounded.
3196 In most cases, you should probably round the results explicitly using one of
3197 L<round()>, L<bround()> or L<bfround()> or by passing the desired accuracy
3198 to the math operation as additional parameter:
3200 my $x = Math::BigInt->new(30000);
3201 my $y = Math::BigInt->new(7);
3202 print scalar $x->copy()->bdiv($y, 2); # print 4300
3203 print scalar $x->copy()->bdiv($y)->bround(2); # print 4300
3207 $x->precision(-2); # local for $x, round at the second digit right of the dot
3208 $x->precision(2); # ditto, round at the second digit left of the dot
3210 CLASS->precision(5); # Global for all members of CLASS
3211 # This also applies to new()!
3212 CLASS->precision(-5); # ditto
3214 $P = CLASS->precision(); # read out global precision
3215 $P = $x->precision(); # read out precision that affects $x
3217 Note: You probably want to use L<accuracy()> instead. With L<accuracy> you
3218 set the number of digits each result should have, with L<precision> you
3219 set the place where to round!
3221 =head1 Autocreating constants
3223 After C<use Math::BigFloat ':constant'> all the floating point constants
3224 in the given scope are converted to C<Math::BigFloat>. This conversion
3225 happens at compile time.
3229 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
3231 prints the value of C<2E-100>. Note that without conversion of
3232 constants the expression 2E-100 will be calculated as normal floating point
3235 Please note that ':constant' does not affect integer constants, nor binary
3236 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
3241 Math with the numbers is done (by default) by a module called
3242 Math::BigInt::Calc. This is equivalent to saying:
3244 use Math::BigFloat lib => 'Calc';
3246 You can change this by using:
3248 use Math::BigFloat lib => 'BitVect';
3250 The following would first try to find Math::BigInt::Foo, then
3251 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
3253 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
3255 Calc.pm uses as internal format an array of elements of some decimal base
3256 (usually 1e7, but this might be different for some systems) with the least
3257 significant digit first, while BitVect.pm uses a bit vector of base 2, most
3258 significant bit first. Other modules might use even different means of
3259 representing the numbers. See the respective module documentation for further
3262 Please note that Math::BigFloat does B<not> use the denoted library itself,
3263 but it merely passes the lib argument to Math::BigInt. So, instead of the need
3266 use Math::BigInt lib => 'GMP';
3269 you can roll it all into one line:
3271 use Math::BigFloat lib => 'GMP';
3273 It is also possible to just require Math::BigFloat:
3275 require Math::BigFloat;
3277 This will load the necessary things (like BigInt) when they are needed, and
3280 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
3281 you ever wanted to know about loading a different library.
3283 =head2 Using Math::BigInt::Lite
3285 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
3288 use Math::BigFloat with => 'Math::BigInt::Lite';
3290 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
3291 can combine these if you want. For instance, you may want to use
3292 Math::BigInt objects in your main script, too.
3296 use Math::BigFloat with => 'Math::BigInt::Lite';
3298 Of course, you can combine this with the C<lib> parameter.
3301 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
3303 There is no need for a "use Math::BigInt;" statement, even if you want to
3304 use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
3305 always loads it. But if you add it, add it B<before>:
3309 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
3311 Notice that the module with the last C<lib> will "win" and thus
3312 it's lib will be used if the lib is available:
3315 use Math::BigInt lib => 'Bar,Baz';
3316 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
3318 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
3319 words, Math::BigFloat will try to retain previously loaded libs when you
3320 don't specify it onem but if you specify one, it will try to load them.
3322 Actually, the lib loading order would be "Bar,Baz,Calc", and then
3323 "Foo,Bar,Baz,Calc", but independent of which lib exists, the result is the
3324 same as trying the latter load alone, except for the fact that one of Bar or
3325 Baz might be loaded needlessly in an intermidiate step (and thus hang around
3326 and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
3327 will still be tried to be loaded, but this is not as time/memory consuming as
3328 actually loading one of them. Still, this type of usage is not recommended due
3331 The old way (loading the lib only in BigInt) still works though:
3334 use Math::BigInt lib => 'Bar,Baz';
3337 You can even load Math::BigInt afterwards:
3341 use Math::BigInt lib => 'Bar,Baz';
3343 But this has the same problems like #5, it will first load Calc
3344 (Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
3345 Baz, depending on which of them works and is usable/loadable. Since this
3346 loads Calc unnec., it is not recommended.
3348 Since it also possible to just require Math::BigFloat, this poses the question
3349 about what libary this will use:
3351 require Math::BigFloat;
3352 my $x = Math::BigFloat->new(123); $x += 123;
3354 It will use Calc. Please note that the call to import() is still done, but
3355 only when you use for the first time some Math::BigFloat math (it is triggered
3356 via any constructor, so the first time you create a Math::BigFloat, the load
3357 will happen in the background). This means:
3359 require Math::BigFloat;
3360 Math::BigFloat->import ( lib => 'Foo,Bar' );
3362 would be the same as:
3364 use Math::BigFloat lib => 'Foo, Bar';
3366 But don't try to be clever to insert some operations in between:
3368 require Math::BigFloat;
3369 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
3370 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
3371 $x = Math::BigFloat->bone()+4; # now use Pari
3373 While this works, it loads Calc needlessly. But maybe you just wanted that?
3375 B<Examples #3 is highly recommended> for daily usage.
3379 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
3385 =item stringify, bstr()
3387 Both stringify and bstr() now drop the leading '+'. The old code would return
3388 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
3389 reasoning and details.
3393 The following will probably not print what you expect:
3395 print $c->bdiv(123.456),"\n";
3397 It prints both quotient and reminder since print works in list context. Also,
3398 bdiv() will modify $c, so be careful. You probably want to use
3400 print $c / 123.456,"\n";
3401 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
3407 The following will probably not print what you expect:
3409 my $c = Math::BigFloat->new('3.14159');
3410 print $c->brsft(3,10),"\n"; # prints 0.00314153.1415
3412 It prints both quotient and remainder, since print calls C<brsft()> in list
3413 context. Also, C<< $c->brsft() >> will modify $c, so be careful.
3414 You probably want to use
3416 print scalar $c->copy()->brsft(3,10),"\n";
3417 # or if you really want to modify $c
3418 print scalar $c->brsft(3,10),"\n";
3422 =item Modifying and =
3426 $x = Math::BigFloat->new(5);
3429 It will not do what you think, e.g. making a copy of $x. Instead it just makes
3430 a second reference to the B<same> object and stores it in $y. Thus anything
3431 that modifies $x will modify $y (except overloaded math operators), and vice
3432 versa. See L<Math::BigInt> for details and how to avoid that.
3436 C<bpow()> now modifies the first argument, unlike the old code which left
3437 it alone and only returned the result. This is to be consistent with
3438 C<badd()> etc. The first will modify $x, the second one won't:
3440 print bpow($x,$i),"\n"; # modify $x
3441 print $x->bpow($i),"\n"; # ditto
3442 print $x ** $i,"\n"; # leave $x alone
3444 =item precision() vs. accuracy()
3446 A common pitfall is to use L<precision()> when you want to round a result to
3447 a certain number of digits:
3451 Math::BigFloat->precision(4); # does not do what you think it does
3452 my $x = Math::BigFloat->new(12345); # rounds $x to "12000"!
3453 print "$x\n"; # print "12000"
3454 my $y = Math::BigFloat->new(3); # rounds $y to "0"!
3455 print "$y\n"; # print "0"
3456 $z = $x / $y; # 12000 / 0 => NaN!
3458 print $z->precision(),"\n"; # 4
3460 Replacing L<precision> with L<accuracy> is probably not what you want, either:
3464 Math::BigFloat->accuracy(4); # enables global rounding:
3465 my $x = Math::BigFloat->new(123456); # rounded immediately to "12350"
3466 print "$x\n"; # print "123500"
3467 my $y = Math::BigFloat->new(3); # rounded to "3
3468 print "$y\n"; # print "3"
3469 print $z = $x->copy()->bdiv($y),"\n"; # 41170
3470 print $z->accuracy(),"\n"; # 4
3472 What you want to use instead is:
3476 my $x = Math::BigFloat->new(123456); # no rounding
3477 print "$x\n"; # print "123456"
3478 my $y = Math::BigFloat->new(3); # no rounding
3479 print "$y\n"; # print "3"
3480 print $z = $x->copy()->bdiv($y,4),"\n"; # 41150
3481 print $z->accuracy(),"\n"; # undef
3483 In addition to computing what you expected, the last example also does B<not>
3484 "taint" the result with an accuracy or precision setting, which would
3485 influence any further operation.
3491 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
3492 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
3494 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
3495 because they solve the autoupgrading/downgrading issue, at least partly.
3497 The package at L<http://search.cpan.org/~tels/Math-BigInt> contains
3498 more documentation including a full version history, testcases, empty
3499 subclass files and benchmarks.
3503 This program is free software; you may redistribute it and/or modify it under
3504 the same terms as Perl itself.
3508 Mark Biggar, overloaded interface by Ilya Zakharevich.
3509 Completely rewritten by Tels L<http://bloodgate.com> in 2001 - 2006, and still