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 and $_trap_nan are internal and should never be accessed from the 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 { $_[2] ?
29 ref($_[0])->bcmp($_[1],$_[0]) :
30 ref($_[0])->bcmp($_[0],$_[1])},
31 'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
34 ##############################################################################
35 # global constants, flags and assorted stuff
37 # the following are public, but their usage is not recommended. Use the
38 # accessor methods instead.
40 # class constants, use Class->constant_name() to access
41 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
48 # the package we are using for our private parts, defaults to:
49 # Math::BigInt->config()->{lib}
50 my $MBI = 'Math::BigInt::Calc';
52 # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
54 # the same for infinity
57 # constant for easier life
60 my $IMPORT = 0; # was import() called yet? used to make require work
62 # some digits of accuracy for blog(undef,10); which we use in blog() for speed
64 '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
65 my $LOG_10_A = length($LOG_10)-1;
68 '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
69 my $LOG_2_A = length($LOG_2)-1;
70 my $HALF = '0.5'; # made into an object if necc.
72 ##############################################################################
73 # the old code had $rnd_mode, so we need to support it, too
75 sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
76 sub FETCH { return $round_mode; }
77 sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
81 # when someone set's $rnd_mode, we catch this and check the value to see
82 # whether it is valid or not.
83 $rnd_mode = 'even'; tie $rnd_mode, 'Math::BigFloat';
86 ##############################################################################
89 # valid method aliases for AUTOLOAD
90 my %methods = map { $_ => 1 }
91 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
92 fint facmp fcmp fzero fnan finf finc fdec flog ffac
93 fceil ffloor frsft flsft fone flog froot
95 # valid method's that can be hand-ed up (for AUTOLOAD)
96 my %hand_ups = map { $_ => 1 }
97 qw / is_nan is_inf is_negative is_positive is_pos is_neg
98 accuracy precision div_scale round_mode fneg fabs fnot
99 objectify upgrade downgrade
103 sub method_alias { exists $methods{$_[0]||''}; }
104 sub method_hand_up { exists $hand_ups{$_[0]||''}; }
107 ##############################################################################
112 # create a new BigFloat object from a string or another bigfloat object.
115 # sign => sign (+/-), or "NaN"
117 my ($class,$wanted,@r) = @_;
119 # avoid numify-calls by not using || on $wanted!
120 return $class->bzero() if !defined $wanted; # default to 0
121 return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
123 $class->import() if $IMPORT == 0; # make require work
125 my $self = {}; bless $self, $class;
126 # shortcut for bigints and its subclasses
127 if ((ref($wanted)) && (ref($wanted) ne $class))
129 $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy
130 $self->{_e} = $MBI->_zero();
132 $self->{sign} = $wanted->sign();
133 return $self->bnorm();
136 # handle '+inf', '-inf' first
137 if ($wanted =~ /^[+-]?inf$/)
139 return $downgrade->new($wanted) if $downgrade;
141 $self->{_e} = $MBI->_zero();
143 $self->{_m} = $MBI->_zero();
144 $self->{sign} = $wanted;
145 $self->{sign} = '+inf' if $self->{sign} eq 'inf';
146 return $self->bnorm();
149 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
155 Carp::croak ("$wanted is not a number initialized to $class");
158 return $downgrade->bnan() if $downgrade;
160 $self->{_e} = $MBI->_zero();
162 $self->{_m} = $MBI->_zero();
163 $self->{sign} = $nan;
167 # make integer from mantissa by adjusting exp, then convert to int
168 $self->{_e} = $MBI->_new($$ev); # exponent
169 $self->{_es} = $$es || '+';
170 my $mantissa = "$$miv$$mfv"; # create mant.
171 $mantissa =~ s/^0+(\d)/$1/; # strip leading zeros
172 $self->{_m} = $MBI->_new($mantissa); # create mant.
174 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
175 if (CORE::length($$mfv) != 0)
177 my $len = $MBI->_new( CORE::length($$mfv));
178 ($self->{_e}, $self->{_es}) =
179 _e_sub ($self->{_e}, $len, $self->{_es}, '+');
181 $self->{sign} = $$mis;
183 # we can only have trailing zeros on the mantissa of $$mfv eq ''
184 if (CORE::length($$mfv) == 0)
186 my $zeros = $MBI->_zeros($self->{_m}); # correct for trailing zeros
189 my $z = $MBI->_new($zeros);
190 $MBI->_rsft ( $self->{_m}, $z, 10);
191 _e_add ( $self->{_e}, $z, $self->{_es}, '+');
194 # for something like 0Ey, set y to 1, and -0 => +0
195 $self->{sign} = '+', $self->{_e} = $MBI->_one()
196 if $MBI->_is_zero($self->{_m});
197 return $self->round(@r) if !$downgrade;
199 # if downgrade, inf, NaN or integers go down
201 if ($downgrade && $self->{_es} eq '+')
203 if ($MBI->_is_zero( $self->{_e} ))
205 return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
207 return $downgrade->new($self->bsstr());
209 $self->bnorm()->round(@r); # first normalize, then round
217 # if two arguments, the first one is the class to "swallow" subclasses
225 return unless ref($x); # only for objects
227 my $self = {}; bless $self,$c;
229 $self->{sign} = $x->{sign};
230 $self->{_es} = $x->{_es};
231 $self->{_m} = $MBI->_copy($x->{_m});
232 $self->{_e} = $MBI->_copy($x->{_e});
233 $self->{_a} = $x->{_a} if defined $x->{_a};
234 $self->{_p} = $x->{_p} if defined $x->{_p};
240 # used by parent class bone() to initialize number to NaN
246 my $class = ref($self);
247 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
250 $IMPORT=1; # call our import only once
251 $self->{_m} = $MBI->_zero();
252 $self->{_e} = $MBI->_zero();
258 # used by parent class bone() to initialize number to +-inf
264 my $class = ref($self);
265 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
268 $IMPORT=1; # call our import only once
269 $self->{_m} = $MBI->_zero();
270 $self->{_e} = $MBI->_zero();
276 # used by parent class bone() to initialize number to 1
278 $IMPORT=1; # call our import only once
279 $self->{_m} = $MBI->_one();
280 $self->{_e} = $MBI->_zero();
286 # used by parent class bone() to initialize number to 0
288 $IMPORT=1; # call our import only once
289 $self->{_m} = $MBI->_zero();
290 $self->{_e} = $MBI->_one();
296 my ($self,$class) = @_;
297 return if $class =~ /^Math::BigInt/; # we aren't one of these
298 UNIVERSAL::isa($self,$class);
303 # return (later set?) configuration data as hash ref
304 my $class = shift || 'Math::BigFloat';
306 my $cfg = $class->SUPER::config(@_);
308 # now we need only to override the ones that are different from our parent
309 $cfg->{class} = $class;
314 ##############################################################################
315 # string conversation
319 # (ref to BFLOAT or num_str ) return num_str
320 # Convert number from internal format to (non-scientific) string format.
321 # internal format is always normalized (no leading zeros, "-0" => "+0")
322 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
324 if ($x->{sign} !~ /^[+-]$/)
326 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
330 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
333 my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
336 $es = $MBI->_str($x->{_m});
337 $len = CORE::length($es);
338 my $e = $MBI->_num($x->{_e});
339 $e = -$e if $x->{_es} eq '-';
343 # if _e is bigger than a scalar, the following will blow your memory
346 my $r = abs($e) - $len;
347 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
351 substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
352 $cad = -$cad if $x->{_es} eq '-';
358 $es .= '0' x $e; $len += $e; $cad = 0;
362 $es = '-'.$es if $x->{sign} eq '-';
363 # if set accuracy or precision, pad with zeros on the right side
364 if ((defined $x->{_a}) && ($not_zero))
366 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
367 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
368 $zeros = $x->{_a} - $len if $cad != $len;
369 $es .= $dot.'0' x $zeros if $zeros > 0;
371 elsif ((($x->{_p} || 0) < 0))
373 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
374 my $zeros = -$x->{_p} + $cad;
375 $es .= $dot.'0' x $zeros if $zeros > 0;
382 # (ref to BFLOAT or num_str ) return num_str
383 # Convert number from internal format to scientific string format.
384 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
385 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
387 if ($x->{sign} !~ /^[+-]$/)
389 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
392 my $sep = 'e'.$x->{_es};
393 my $sign = $x->{sign}; $sign = '' if $sign eq '+';
394 $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
399 # Make a number from a BigFloat object
400 # simple return a string and let Perl's atoi()/atof() handle the rest
401 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
405 ##############################################################################
406 # public stuff (usually prefixed with "b")
409 # XXX TODO this must be overwritten and return NaN for non-integer values
410 # band(), bior(), bxor(), too
413 # $class->SUPER::bnot($class,@_);
418 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
421 my ($self,$x,$y) = (ref($_[0]),@_);
422 # objectify is costly, so avoid it
423 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
425 ($self,$x,$y) = objectify(2,@_);
428 return $upgrade->bcmp($x,$y) if defined $upgrade &&
429 ((!$x->isa($self)) || (!$y->isa($self)));
431 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
433 # handle +-inf and NaN
434 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
435 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
436 return +1 if $x->{sign} eq '+inf';
437 return -1 if $x->{sign} eq '-inf';
438 return -1 if $y->{sign} eq '+inf';
442 # check sign for speed first
443 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
444 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
447 my $xz = $x->is_zero();
448 my $yz = $y->is_zero();
449 return 0 if $xz && $yz; # 0 <=> 0
450 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
451 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
453 # adjust so that exponents are equal
454 my $lxm = $MBI->_len($x->{_m});
455 my $lym = $MBI->_len($y->{_m});
456 # the numify somewhat limits our length, but makes it much faster
457 my ($xes,$yes) = (1,1);
458 $xes = -1 if $x->{_es} ne '+';
459 $yes = -1 if $y->{_es} ne '+';
460 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
461 my $ly = $lym + $yes * $MBI->_num($y->{_e});
462 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
463 return $l <=> 0 if $l != 0;
465 # lengths (corrected by exponent) are equal
466 # so make mantissa equal length by padding with zero (shift left)
467 my $diff = $lxm - $lym;
468 my $xm = $x->{_m}; # not yet copy it
472 $ym = $MBI->_copy($y->{_m});
473 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
477 $xm = $MBI->_copy($x->{_m});
478 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
480 my $rc = $MBI->_acmp($xm,$ym);
481 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
487 # Compares 2 values, ignoring their signs.
488 # Returns one of undef, <0, =0, >0. (suitable for sort)
491 my ($self,$x,$y) = (ref($_[0]),@_);
492 # objectify is costly, so avoid it
493 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
495 ($self,$x,$y) = objectify(2,@_);
498 return $upgrade->bacmp($x,$y) if defined $upgrade &&
499 ((!$x->isa($self)) || (!$y->isa($self)));
501 # handle +-inf and NaN's
502 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
504 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
505 return 0 if ($x->is_inf() && $y->is_inf());
506 return 1 if ($x->is_inf() && !$y->is_inf());
511 my $xz = $x->is_zero();
512 my $yz = $y->is_zero();
513 return 0 if $xz && $yz; # 0 <=> 0
514 return -1 if $xz && !$yz; # 0 <=> +y
515 return 1 if $yz && !$xz; # +x <=> 0
517 # adjust so that exponents are equal
518 my $lxm = $MBI->_len($x->{_m});
519 my $lym = $MBI->_len($y->{_m});
520 my ($xes,$yes) = (1,1);
521 $xes = -1 if $x->{_es} ne '+';
522 $yes = -1 if $y->{_es} ne '+';
523 # the numify somewhat limits our length, but makes it much faster
524 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
525 my $ly = $lym + $yes * $MBI->_num($y->{_e});
527 return $l <=> 0 if $l != 0;
529 # lengths (corrected by exponent) are equal
530 # so make mantissa equal-length by padding with zero (shift left)
531 my $diff = $lxm - $lym;
532 my $xm = $x->{_m}; # not yet copy it
536 $ym = $MBI->_copy($y->{_m});
537 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
541 $xm = $MBI->_copy($x->{_m});
542 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
544 $MBI->_acmp($xm,$ym);
549 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
550 # return result as BFLOAT
553 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
554 # objectify is costly, so avoid it
555 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
557 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
560 # inf and NaN handling
561 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
564 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
566 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
568 # +inf++inf or -inf+-inf => same, rest is NaN
569 return $x if $x->{sign} eq $y->{sign};
572 # +-inf + something => +inf; something +-inf => +-inf
573 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
577 return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
578 ((!$x->isa($self)) || (!$y->isa($self)));
580 # speed: no add for 0+y or x+0
581 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
582 if ($x->is_zero()) # 0+y
584 # make copy, clobbering up x (modify in place!)
585 $x->{_e} = $MBI->_copy($y->{_e});
586 $x->{_es} = $y->{_es};
587 $x->{_m} = $MBI->_copy($y->{_m});
588 $x->{sign} = $y->{sign} || $nan;
589 return $x->round($a,$p,$r,$y);
592 # take lower of the two e's and adapt m1 to it to match m2
594 $e = $MBI->_zero() if !defined $e; # if no BFLOAT?
595 $e = $MBI->_copy($e); # make copy (didn't do it yet)
599 ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
601 my $add = $MBI->_copy($y->{_m});
603 if ($es eq '-') # < 0
605 $MBI->_lsft( $x->{_m}, $e, 10);
606 ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
608 elsif (!$MBI->_is_zero($e)) # > 0
610 $MBI->_lsft($add, $e, 10);
612 # else: both e are the same, so just leave them
614 if ($x->{sign} eq $y->{sign})
617 $x->{_m} = $MBI->_add($x->{_m}, $add);
621 ($x->{_m}, $x->{sign}) =
622 _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
625 # delete trailing zeros, then round
626 $x->bnorm()->round($a,$p,$r,$y);
631 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
632 # subtract second arg from first, modify first
635 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
636 # objectify is costly, so avoid it
637 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
639 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
642 if ($y->is_zero()) # still round for not adding zero
644 return $x->round($a,$p,$r);
648 $y->{sign} =~ tr/+-/-+/; # does nothing for NaN
649 $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
650 $y->{sign} =~ tr/+-/-+/; # refix $y (does nothing for NaN)
651 $x; # already rounded by badd()
656 # increment arg by one
657 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
659 if ($x->{_es} eq '-')
661 return $x->badd($self->bone(),@r); # digits after dot
664 if (!$MBI->_is_zero($x->{_e})) # _e == 0 for NaN, inf, -inf
666 # 1e2 => 100, so after the shift below _m has a '0' as last digit
667 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
668 $x->{_e} = $MBI->_zero(); # normalize
670 # we know that the last digit of $x will be '1' or '9', depending on the
674 if ($x->{sign} eq '+')
676 $MBI->_inc($x->{_m});
677 return $x->bnorm()->bround(@r);
679 elsif ($x->{sign} eq '-')
681 $MBI->_dec($x->{_m});
682 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
683 return $x->bnorm()->bround(@r);
685 # inf, nan handling etc
686 $x->badd($self->bone(),@r); # badd() does round
691 # decrement arg by one
692 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
694 if ($x->{_es} eq '-')
696 return $x->badd($self->bone('-'),@r); # digits after dot
699 if (!$MBI->_is_zero($x->{_e}))
701 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
702 $x->{_e} = $MBI->_zero(); # normalize
706 my $zero = $x->is_zero();
708 if (($x->{sign} eq '-') || $zero)
710 $MBI->_inc($x->{_m});
711 $x->{sign} = '-' if $zero; # 0 => 1 => -1
712 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
713 return $x->bnorm()->round(@r);
716 elsif ($x->{sign} eq '+')
718 $MBI->_dec($x->{_m});
719 return $x->bnorm()->round(@r);
721 # inf, nan handling etc
722 $x->badd($self->bone('-'),@r); # does round
729 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
731 # $base > 0, $base != 1; if $base == undef default to $base == e
734 # we need to limit the accuracy to protect against overflow
737 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
739 # also takes care of the "error in _find_round_parameters?" case
740 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
743 # no rounding at all, so must use fallback
744 if (scalar @params == 0)
746 # simulate old behaviour
747 $params[0] = $self->div_scale(); # and round to it as accuracy
748 $params[1] = undef; # P = undef
749 $scale = $params[0]+4; # at least four more for proper round
750 $params[2] = $r; # round mode by caller or undef
751 $fallback = 1; # to clear a/p afterwards
755 # the 4 below is empirical, and there might be cases where it is not
757 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
760 return $x->bzero(@params) if $x->is_one();
761 # base not defined => base == Euler's constant e
764 # make object, since we don't feed it through objectify() to still get the
765 # case of $base == undef
766 $base = $self->new($base) unless ref($base);
767 # $base > 0; $base != 1
768 return $x->bnan() if $base->is_zero() || $base->is_one() ||
769 $base->{sign} ne '+';
770 # if $x == $base, we know the result must be 1.0
771 return $x->bone('+',@params) if $x->bcmp($base) == 0;
774 # when user set globals, they would interfere with our calculation, so
775 # disable them and later re-enable them
777 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
778 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
779 # we also need to disable any set A or P on $x (_find_round_parameters took
780 # them already into account), since these would interfere, too
781 delete $x->{_a}; delete $x->{_p};
782 # need to disable $upgrade in BigInt, to avoid deep recursion
783 local $Math::BigInt::upgrade = undef;
784 local $Math::BigFloat::downgrade = undef;
786 # upgrade $x if $x is not a BigFloat (handle BigInt input)
787 if (!$x->isa('Math::BigFloat'))
789 $x = Math::BigFloat->new($x);
795 # If the base is defined and an integer, try to calculate integer result
796 # first. This is very fast, and in case the real result was found, we can
798 if (defined $base && $base->is_int() && $x->is_int())
800 my $i = $MBI->_copy( $x->{_m} );
801 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
802 my $int = Math::BigInt->bzero();
804 $int->blog($base->as_number());
806 if ($base->as_number()->bpow($int) == $x)
808 # found result, return it
809 $x->{_m} = $int->{value};
810 $x->{_e} = $MBI->_zero();
819 # first calculate the log to base e (using reduction by 10 (and probably 2))
820 $self->_log_10($x,$scale);
822 # and if a different base was requested, convert it
825 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
826 # not ln, but some other base (don't modify $base)
827 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
831 # shortcut to not run through _find_round_parameters again
832 if (defined $params[0])
834 $x->bround($params[0],$params[2]); # then round accordingly
838 $x->bfround($params[1],$params[2]); # then round accordingly
842 # clear a/p after round, since user did not request it
843 delete $x->{_a}; delete $x->{_p};
846 $$abr = $ab; $$pbr = $pb;
853 # internal log function to calculate ln() based on Taylor series.
854 # Modifies $x in place.
855 my ($self,$x,$scale) = @_;
857 # in case of $x == 1, result is 0
858 return $x->bzero() if $x->is_one();
860 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
864 # Taylor: | u 1 u^3 1 u^5 |
865 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
866 # |_ v 3 v^3 5 v^5 _|
868 # This takes much more steps to calculate the result and is thus not used
871 # Taylor: | u 1 u^2 1 u^3 |
872 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
873 # |_ x 2 x^2 3 x^3 _|
875 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
877 $v = $x->copy(); $v->binc(); # v = x+1
878 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
879 $x->bdiv($v,$scale); # first term: u/v
882 $u *= $u; $v *= $v; # u^2, v^2
883 $below->bmul($v); # u^3, v^3
885 $factor = $self->new(3); $f = $self->new(2);
887 my $steps = 0 if DEBUG;
888 $limit = $self->new("1E-". ($scale-1));
891 # we calculate the next term, and add it to the last
892 # when the next term is below our limit, it won't affect the outcome
893 # anymore, so we stop
895 # calculating the next term simple from over/below will result in quite
896 # a time hog if the input has many digits, since over and below will
897 # accumulate more and more digits, and the result will also have many
898 # digits, but in the end it is rounded to $scale digits anyway. So if we
899 # round $over and $below first, we save a lot of time for the division
900 # (not with log(1.2345), but try log (123**123) to see what I mean. This
901 # can introduce a rounding error if the division result would be f.i.
902 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
903 # if we truncated $over and $below we might get 0.12345. Does this matter
904 # for the end result? So we give $over and $below 4 more digits to be
905 # on the safe side (unscientific error handling as usual... :+D
907 $next = $over->copy->bround($scale+4)->bdiv(
908 $below->copy->bmul($factor)->bround($scale+4),
912 ## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
914 last if $next->bacmp($limit) <= 0;
916 delete $next->{_a}; delete $next->{_p};
918 # calculate things for the next term
919 $over *= $u; $below *= $v; $factor->badd($f);
922 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
925 $x->bmul($f); # $x *= 2
926 print "took $steps steps\n" if DEBUG;
931 # Internal log function based on reducing input to the range of 0.1 .. 9.99
932 # and then "correcting" the result to the proper one. Modifies $x in place.
933 my ($self,$x,$scale) = @_;
935 # taking blog() from numbers greater than 10 takes a *very long* time, so we
936 # break the computation down into parts based on the observation that:
937 # blog(x*y) = blog(x) + blog(y)
938 # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
939 # the faster it get's, especially because 2*$x takes about 10 times as long,
940 # so by dividing $x by 10 we make it at least factor 100 faster...)
942 # The same observation is valid for numbers smaller than 0.1 (e.g. computing
943 # log(1) is fastest, and the farther away we get from 1, the longer it takes)
944 # so we also 'break' this down by multiplying $x with 10 and subtract the
945 # log(10) afterwards to get the correct result.
947 # calculate nr of digits before dot
948 my $dbd = $MBI->_num($x->{_e});
949 $dbd = -$dbd if $x->{_es} eq '-';
950 $dbd += $MBI->_len($x->{_m});
952 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
955 my $calc = 1; # do some calculation?
957 # disable the shortcut for 10, since we need log(10) and this would recurse
959 if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
961 $dbd = 0; # disable shortcut
962 # we can use the cached value in these cases
963 if ($scale <= $LOG_10_A)
965 $x->bzero(); $x->badd($LOG_10);
966 $calc = 0; # no need to calc, but round
971 # disable the shortcut for 2, since we maybe have it cached
972 if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
974 $dbd = 0; # disable shortcut
975 # we can use the cached value in these cases
976 if ($scale <= $LOG_2_A)
978 $x->bzero(); $x->badd($LOG_2);
979 $calc = 0; # no need to calc, but round
984 # if $x = 0.1, we know the result must be 0-log(10)
985 if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
986 $MBI->_is_one($x->{_m}))
988 $dbd = 0; # disable shortcut
989 # we can use the cached value in these cases
990 if ($scale <= $LOG_10_A)
992 $x->bzero(); $x->bsub($LOG_10);
993 $calc = 0; # no need to calc, but round
997 return if $calc == 0; # already have the result
999 # default: these correction factors are undef and thus not used
1000 my $l_10; # value of ln(10) to A of $scale
1001 my $l_2; # value of ln(2) to A of $scale
1003 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1004 # so don't do this shortcut for 1 or 0
1005 if (($dbd > 1) || ($dbd < 0))
1007 # convert our cached value to an object if not already (avoid doing this
1008 # at import() time, since not everybody needs this)
1009 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1011 #print "x = $x, dbd = $dbd, calc = $calc\n";
1012 # got more than one digit before the dot, or more than one zero after the
1014 # log(123) == log(1.23) + log(10) * 2
1015 # log(0.0123) == log(1.23) - log(10) * 2
1017 if ($scale <= $LOG_10_A)
1020 $l_10 = $LOG_10->copy(); # copy for mul
1024 # else: slower, compute it (but don't cache it, because it could be big)
1025 # also disable downgrade for this code path
1026 local $Math::BigFloat::downgrade = undef;
1027 $l_10 = $self->new(10)->blog(undef,$scale); # scale+4, actually
1029 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
1030 $l_10->bmul( $self->new($dbd)); # log(10) * (digits_before_dot-1)
1037 ($x->{_e}, $x->{_es}) =
1038 _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
1042 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1044 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1045 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1047 $HALF = $self->new($HALF) unless ref($HALF);
1049 my $twos = 0; # default: none (0 times)
1050 my $two = $self->new(2);
1051 while ($x->bacmp($HALF) <= 0)
1053 $twos--; $x->bmul($two);
1055 while ($x->bacmp($two) >= 0)
1057 $twos++; $x->bdiv($two,$scale+4); # keep all digits
1059 # $twos > 0 => did mul 2, < 0 => did div 2 (never both)
1060 # calculate correction factor based on ln(2)
1063 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1064 if ($scale <= $LOG_2_A)
1067 $l_2 = $LOG_2->copy(); # copy for mul
1071 # else: slower, compute it (but don't cache it, because it could be big)
1072 # also disable downgrade for this code path
1073 local $Math::BigFloat::downgrade = undef;
1074 $l_2 = $two->blog(undef,$scale); # scale+4, actually
1076 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
1079 $self->_log($x,$scale); # need to do the "normal" way
1080 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
1081 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
1082 # all done, $x contains now the result
1087 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1088 # does not modify arguments, but returns new object
1089 # Lowest Common Multiplicator
1091 my ($self,@arg) = objectify(0,@_);
1092 my $x = $self->new(shift @arg);
1093 while (@arg) { $x = _lcm($x,shift @arg); }
1099 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
1100 # does not modify arguments, but returns new object
1101 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
1103 my ($self,@arg) = objectify(0,@_);
1104 my $x = $self->new(shift @arg);
1105 while (@arg) { $x = _gcd($x,shift @arg); }
1109 ##############################################################################
1113 # Internal helper sub to take two positive integers and their signs and
1114 # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1115 # output ($CALC,('+'|'-'))
1116 my ($x,$y,$xs,$ys) = @_;
1118 # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1121 $x = $MBI->_add ($x, $y ); # a+b
1122 # the sign follows $xs
1126 my $a = $MBI->_acmp($x,$y);
1129 $x = $MBI->_sub ($x , $y); # abs sub
1133 $x = $MBI->_zero(); # result is 0
1138 $x = $MBI->_sub ( $y, $x, 1 ); # abs sub
1146 # Internal helper sub to take two positive integers and their signs and
1147 # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1148 # output ($CALC,('+'|'-'))
1149 my ($x,$y,$xs,$ys) = @_;
1153 _e_add($x,$y,$xs,$ys); # call add (does subtract now)
1156 ###############################################################################
1157 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1161 # return true if arg (BFLOAT or num_str) is an integer
1162 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1164 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1165 $x->{_es} eq '+'; # 1e-1 => no integer
1171 # return true if arg (BFLOAT or num_str) is zero
1172 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1174 return 1 if $x->{sign} eq '+' && $MBI->_is_zero($x->{_m});
1180 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1181 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1183 $sign = '+' if !defined $sign || $sign ne '-';
1185 if ($x->{sign} eq $sign &&
1186 $MBI->_is_zero($x->{_e}) && $MBI->_is_one($x->{_m}));
1192 # return true if arg (BFLOAT or num_str) is odd or false if even
1193 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1195 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1196 ($MBI->_is_zero($x->{_e}) && $MBI->_is_odd($x->{_m}));
1202 # return true if arg (BINT or num_str) is even or false if odd
1203 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1205 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
1206 return 1 if ($x->{_es} eq '+' # 123.45 is never
1207 && $MBI->_is_even($x->{_m})); # but 1200 is
1213 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
1214 # (BINT or num_str, BINT or num_str) return BINT
1217 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1218 # objectify is costly, so avoid it
1219 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1221 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1224 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1227 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1229 return $x->bnan() if $x->is_zero() || $y->is_zero();
1230 # result will always be +-inf:
1231 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1232 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1233 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1234 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1235 return $x->binf('-');
1238 return $x->bzero() if $x->is_zero() || $y->is_zero();
1240 return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
1241 ((!$x->isa($self)) || (!$y->isa($self)));
1243 # aEb * cEd = (a*c)E(b+d)
1244 $MBI->_mul($x->{_m},$y->{_m});
1245 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1248 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1249 return $x->bnorm()->round($a,$p,$r,$y);
1254 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1255 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1258 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1259 # objectify is costly, so avoid it
1260 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1262 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1265 return $self->_div_inf($x,$y)
1266 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1268 # x== 0 # also: or y == 1 or y == -1
1269 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1272 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1274 # we need to limit the accuracy to protect against overflow
1276 my (@params,$scale);
1277 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1279 return $x if $x->is_nan(); # error in _find_round_parameters?
1281 # no rounding at all, so must use fallback
1282 if (scalar @params == 0)
1284 # simulate old behaviour
1285 $params[0] = $self->div_scale(); # and round to it as accuracy
1286 $scale = $params[0]+4; # at least four more for proper round
1287 $params[2] = $r; # round mode by caller or undef
1288 $fallback = 1; # to clear a/p afterwards
1292 # the 4 below is empirical, and there might be cases where it is not
1294 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1296 my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
1297 $scale = $lx if $lx > $scale;
1298 $scale = $ly if $ly > $scale;
1299 my $diff = $ly - $lx;
1300 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
1302 # make copy of $x in case of list context for later reminder calculation
1304 if (wantarray && !$y->is_one())
1309 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1311 # check for / +-1 ( +/- 1E0)
1314 # promote BigInts and it's subclasses (except when already a BigFloat)
1315 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1317 # calculate the result to $scale digits and then round it
1318 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1319 $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1320 $MBI->_div ($x->{_m},$y->{_m} ); # a/c
1322 ($x->{_e},$x->{_es}) =
1323 _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1324 # correct for 10**scale
1325 ($x->{_e},$x->{_es}) =
1326 _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1327 $x->bnorm(); # remove trailing 0's
1330 # shortcut to not run through _find_round_parameters again
1331 if (defined $params[0])
1333 delete $x->{_a}; # clear before round
1334 $x->bround($params[0],$params[2]); # then round accordingly
1338 delete $x->{_p}; # clear before round
1339 $x->bfround($params[1],$params[2]); # then round accordingly
1343 # clear a/p after round, since user did not request it
1344 delete $x->{_a}; delete $x->{_p};
1351 $rem->bmod($y,@params); # copy already done
1355 $rem = $self->bzero();
1359 # clear a/p after round, since user did not request it
1360 delete $rem->{_a}; delete $rem->{_p};
1369 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1372 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1373 # objectify is costly, so avoid it
1374 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1376 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1379 # handle NaN, inf, -inf
1380 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1382 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1383 $x->{sign} = $re->{sign};
1384 $x->{_e} = $re->{_e};
1385 $x->{_m} = $re->{_m};
1386 return $x->round($a,$p,$r,$y);
1390 return $x->bnan() if $x->is_zero();
1393 return $x->bzero() if $y->is_one() || $x->is_zero();
1395 my $cmp = $x->bacmp($y); # equal or $x < $y?
1396 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1398 # only $y of the operands negative?
1399 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1401 $x->{sign} = $y->{sign}; # calc sign first
1402 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1404 my $ym = $MBI->_copy($y->{_m});
1407 $MBI->_lsft( $ym, $y->{_e}, 10)
1408 if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1410 # if $y has digits after dot
1411 my $shifty = 0; # correct _e of $x by this
1412 if ($y->{_es} eq '-') # has digits after dot
1414 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1415 $shifty = $MBI->_num($y->{_e}); # no more digits after dot
1416 $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1418 # $ym is now mantissa of $y based on exponent 0
1420 my $shiftx = 0; # correct _e of $x by this
1421 if ($x->{_es} eq '-') # has digits after dot
1423 # 123.4 % 20 => 1234 % 200
1424 $shiftx = $MBI->_num($x->{_e}); # no more digits after dot
1425 $MBI->_lsft($ym, $x->{_e}, 10); # 123 => 1230
1427 # 123e1 % 20 => 1230 % 20
1428 if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1430 $MBI->_lsft( $x->{_m}, $x->{_e},10); # es => '+' here
1433 $x->{_e} = $MBI->_new($shiftx);
1435 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1436 $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1438 # now mantissas are equalized, exponent of $x is adjusted, so calc result
1440 $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1442 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1445 if ($neg != 0) # one of them negative => correct in place
1448 $x->{_m} = $r->{_m};
1449 $x->{_e} = $r->{_e};
1450 $x->{_es} = $r->{_es};
1451 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1455 $x->round($a,$p,$r,$y); # round and return
1460 # calculate $y'th root of $x
1463 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1464 # objectify is costly, so avoid it
1465 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1467 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1470 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1471 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1472 $y->{sign} !~ /^\+$/;
1474 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1476 # we need to limit the accuracy to protect against overflow
1478 my (@params,$scale);
1479 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1481 return $x if $x->is_nan(); # error in _find_round_parameters?
1483 # no rounding at all, so must use fallback
1484 if (scalar @params == 0)
1486 # simulate old behaviour
1487 $params[0] = $self->div_scale(); # and round to it as accuracy
1488 $scale = $params[0]+4; # at least four more for proper round
1489 $params[2] = $r; # iound mode by caller or undef
1490 $fallback = 1; # to clear a/p afterwards
1494 # the 4 below is empirical, and there might be cases where it is not
1496 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1499 # when user set globals, they would interfere with our calculation, so
1500 # disable them and later re-enable them
1502 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1503 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1504 # we also need to disable any set A or P on $x (_find_round_parameters took
1505 # them already into account), since these would interfere, too
1506 delete $x->{_a}; delete $x->{_p};
1507 # need to disable $upgrade in BigInt, to avoid deep recursion
1508 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1510 # remember sign and make $x positive, since -4 ** (1/2) => -2
1511 my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
1514 if ($y->isa('Math::BigFloat'))
1516 $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
1520 $is_two = ($y == 2);
1523 # normal square root if $y == 2:
1526 $x->bsqrt($scale+4);
1528 elsif ($y->is_one('-'))
1531 my $u = $self->bone()->bdiv($x,$scale);
1532 # copy private parts over
1533 $x->{_m} = $u->{_m};
1534 $x->{_e} = $u->{_e};
1535 $x->{_es} = $u->{_es};
1539 # calculate the broot() as integer result first, and if it fits, return
1540 # it rightaway (but only if $x and $y are integer):
1542 my $done = 0; # not yet
1543 if ($y->is_int() && $x->is_int())
1545 my $i = $MBI->_copy( $x->{_m} );
1546 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1547 my $int = Math::BigInt->bzero();
1549 $int->broot($y->as_number());
1551 if ($int->copy()->bpow($y) == $x)
1553 # found result, return it
1554 $x->{_m} = $int->{value};
1555 $x->{_e} = $MBI->_zero();
1563 my $u = $self->bone()->bdiv($y,$scale+4);
1564 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
1565 $x->bpow($u,$scale+4); # el cheapo
1568 $x->bneg() if $sign == 1;
1570 # shortcut to not run through _find_round_parameters again
1571 if (defined $params[0])
1573 $x->bround($params[0],$params[2]); # then round accordingly
1577 $x->bfround($params[1],$params[2]); # then round accordingly
1581 # clear a/p after round, since user did not request it
1582 delete $x->{_a}; delete $x->{_p};
1585 $$abr = $ab; $$pbr = $pb;
1591 # calculate square root
1592 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1594 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
1595 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
1596 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
1598 # we need to limit the accuracy to protect against overflow
1600 my (@params,$scale);
1601 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1603 return $x if $x->is_nan(); # error in _find_round_parameters?
1605 # no rounding at all, so must use fallback
1606 if (scalar @params == 0)
1608 # simulate old behaviour
1609 $params[0] = $self->div_scale(); # and round to it as accuracy
1610 $scale = $params[0]+4; # at least four more for proper round
1611 $params[2] = $r; # round mode by caller or undef
1612 $fallback = 1; # to clear a/p afterwards
1616 # the 4 below is empirical, and there might be cases where it is not
1618 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1621 # when user set globals, they would interfere with our calculation, so
1622 # disable them and later re-enable them
1624 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1625 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1626 # we also need to disable any set A or P on $x (_find_round_parameters took
1627 # them already into account), since these would interfere, too
1628 delete $x->{_a}; delete $x->{_p};
1629 # need to disable $upgrade in BigInt, to avoid deep recursion
1630 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1632 my $i = $MBI->_copy( $x->{_m} );
1633 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1634 my $xas = Math::BigInt->bzero();
1637 my $gs = $xas->copy()->bsqrt(); # some guess
1639 if (($x->{_es} ne '-') # guess can't be accurate if there are
1640 # digits after the dot
1641 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
1643 # exact result, copy result over to keep $x
1644 $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
1646 # shortcut to not run through _find_round_parameters again
1647 if (defined $params[0])
1649 $x->bround($params[0],$params[2]); # then round accordingly
1653 $x->bfround($params[1],$params[2]); # then round accordingly
1657 # clear a/p after round, since user did not request it
1658 delete $x->{_a}; delete $x->{_p};
1660 # re-enable A and P, upgrade is taken care of by "local"
1661 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1665 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1666 # of the result by multipyling the input by 100 and then divide the integer
1667 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
1669 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
1670 my $y1 = $MBI->_copy($x->{_m});
1672 my $length = $MBI->_len($y1);
1674 # Now calculate how many digits the result of sqrt(y1) would have
1675 my $digits = int($length / 2);
1677 # But we need at least $scale digits, so calculate how many are missing
1678 my $shift = $scale - $digits;
1680 # That should never happen (we take care of integer guesses above)
1681 # $shift = 0 if $shift < 0;
1683 # Multiply in steps of 100, by shifting left two times the "missing" digits
1684 my $s2 = $shift * 2;
1686 # We now make sure that $y1 has the same odd or even number of digits than
1687 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1688 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1689 # steps of 10. The length of $x does not count, since an even or odd number
1690 # of digits before the dot is not changed by adding an even number of digits
1691 # after the dot (the result is still odd or even digits long).
1692 $s2++ if $MBI->_is_odd($x->{_e});
1694 $MBI->_lsft( $y1, $MBI->_new($s2), 10);
1696 # now take the square root and truncate to integer
1697 $y1 = $MBI->_sqrt($y1);
1699 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1700 # result, which is than later rounded to the desired scale.
1702 # calculate how many zeros $x had after the '.' (or before it, depending
1703 # on sign of $dat, the result should have half as many:
1704 my $dat = $MBI->_num($x->{_e});
1705 $dat = -$dat if $x->{_es} eq '-';
1710 # no zeros after the dot (e.g. 1.23, 0.49 etc)
1711 # preserve half as many digits before the dot than the input had
1712 # (but round this "up")
1713 $dat = int(($dat+1)/2);
1717 $dat = int(($dat)/2);
1719 $dat -= $MBI->_len($y1);
1723 $x->{_e} = $MBI->_new( $dat );
1728 $x->{_e} = $MBI->_new( $dat );
1734 # shortcut to not run through _find_round_parameters again
1735 if (defined $params[0])
1737 $x->bround($params[0],$params[2]); # then round accordingly
1741 $x->bfround($params[1],$params[2]); # then round accordingly
1745 # clear a/p after round, since user did not request it
1746 delete $x->{_a}; delete $x->{_p};
1749 $$abr = $ab; $$pbr = $pb;
1755 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1756 # compute factorial number, modifies first argument
1759 my ($self,$x,@r) = (ref($_[0]),@_);
1760 # objectify is costly, so avoid it
1761 ($self,$x,@r) = objectify(1,@_) if !ref($x);
1763 return $x if $x->{sign} eq '+inf'; # inf => inf
1765 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1766 ($x->{_es} ne '+')); # digits after dot?
1768 # use BigInt's bfac() for faster calc
1769 if (! $MBI->_is_zero($x->{_e}))
1771 $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
1772 $x->{_e} = $MBI->_zero(); # normalize
1775 $MBI->_fac($x->{_m}); # calculate factorial
1776 $x->bnorm()->round(@r); # norm again and round result
1781 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1782 my ($x,$y,$a,$p,$r) = @_;
1785 # if $y == 0.5, it is sqrt($x)
1786 $HALF = $self->new($HALF) unless ref($HALF);
1787 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp($HALF) == 0;
1790 # a ** x == e ** (x * ln a)
1794 # Taylor: | u u^2 u^3 |
1795 # x ** y = 1 + | --- + --- + ----- + ... |
1798 # we need to limit the accuracy to protect against overflow
1800 my ($scale,@params);
1801 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1803 return $x if $x->is_nan(); # error in _find_round_parameters?
1805 # no rounding at all, so must use fallback
1806 if (scalar @params == 0)
1808 # simulate old behaviour
1809 $params[0] = $self->div_scale(); # and round to it as accuracy
1810 $params[1] = undef; # disable P
1811 $scale = $params[0]+4; # at least four more for proper round
1812 $params[2] = $r; # round mode by caller or undef
1813 $fallback = 1; # to clear a/p afterwards
1817 # the 4 below is empirical, and there might be cases where it is not
1819 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1822 # when user set globals, they would interfere with our calculation, so
1823 # disable them and later re-enable them
1825 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1826 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1827 # we also need to disable any set A or P on $x (_find_round_parameters took
1828 # them already into account), since these would interfere, too
1829 delete $x->{_a}; delete $x->{_p};
1830 # need to disable $upgrade in BigInt, to avoid deep recursion
1831 local $Math::BigInt::upgrade = undef;
1833 my ($limit,$v,$u,$below,$factor,$next,$over);
1835 $u = $x->copy()->blog(undef,$scale)->bmul($y);
1836 $v = $self->bone(); # 1
1837 $factor = $self->new(2); # 2
1838 $x->bone(); # first term: 1
1840 $below = $v->copy();
1843 $limit = $self->new("1E-". ($scale-1));
1847 # we calculate the next term, and add it to the last
1848 # when the next term is below our limit, it won't affect the outcome
1849 # anymore, so we stop
1850 $next = $over->copy()->bdiv($below,$scale);
1851 last if $next->bacmp($limit) <= 0;
1853 # calculate things for the next term
1854 $over *= $u; $below *= $factor; $factor->binc();
1856 last if $x->{sign} !~ /^[-+]$/;
1861 # shortcut to not run through _find_round_parameters again
1862 if (defined $params[0])
1864 $x->bround($params[0],$params[2]); # then round accordingly
1868 $x->bfround($params[1],$params[2]); # then round accordingly
1872 # clear a/p after round, since user did not request it
1873 delete $x->{_a}; delete $x->{_p};
1876 $$abr = $ab; $$pbr = $pb;
1882 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1883 # compute power of two numbers, second arg is used as integer
1884 # modifies first argument
1887 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1888 # objectify is costly, so avoid it
1889 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1891 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1894 return $x if $x->{sign} =~ /^[+-]inf$/;
1895 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1896 return $x->bone() if $y->is_zero();
1897 return $x if $x->is_one() || $y->is_one();
1899 return $x->_pow($y,$a,$p,$r) if !$y->is_int(); # non-integer power
1901 my $y1 = $y->as_number()->{value}; # make CALC
1904 if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
1906 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
1907 return $MBI->_is_odd($y1) ? $x : $x->babs(1);
1911 return $x->bone() if $y->is_zero();
1912 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1913 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
1918 $new_sign = $y->is_odd() ? '-' : '+' if ($x->{sign} ne '+');
1920 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1921 $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
1922 $MBI->_mul ($x->{_e}, $y1);
1924 $x->{sign} = $new_sign;
1926 if ($y->{sign} eq '-')
1928 # modify $x in place!
1929 my $z = $x->copy(); $x->bzero()->binc();
1930 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
1932 $x->round($a,$p,$r,$y);
1935 ###############################################################################
1936 # rounding functions
1940 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1941 # $n == 0 means round to integer
1942 # expects and returns normalized numbers!
1943 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1945 return $x if $x->modify('bfround');
1947 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1948 return $x if !defined $scale; # no-op
1950 # never round a 0, +-inf, NaN
1953 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1956 return $x if $x->{sign} !~ /^[+-]$/;
1958 # don't round if x already has lower precision
1959 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1961 $x->{_p} = $scale; # remember round in any case
1962 delete $x->{_a}; # and clear A
1965 # round right from the '.'
1967 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round
1969 $scale = -$scale; # positive for simplicity
1970 my $len = $MBI->_len($x->{_m}); # length of mantissa
1972 # the following poses a restriction on _e, but if _e is bigger than a
1973 # scalar, you got other problems (memory etc) anyway
1974 my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e}))); # digits after dot
1975 my $zad = 0; # zeros after dot
1976 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
1978 # p rint "scale $scale dad $dad zad $zad len $len\n";
1979 # number bsstr len zad dad
1980 # 0.123 123e-3 3 0 3
1981 # 0.0123 123e-4 3 1 4
1984 # 1.2345 12345e-4 5 0 4
1986 # do not round after/right of the $dad
1987 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1989 # round to zero if rounding inside the $zad, but not for last zero like:
1990 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1991 return $x->bzero() if $scale < $zad;
1992 if ($scale == $zad) # for 0.006, scale -3 and trunc
1998 # adjust round-point to be inside mantissa
2001 $scale = $scale-$zad;
2005 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
2006 $scale = $dbd+$scale;
2012 # round left from the '.'
2014 # 123 => 100 means length(123) = 3 - $scale (2) => 1
2016 my $dbt = $MBI->_len($x->{_m});
2018 my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
2019 # should be the same, so treat it as this
2020 $scale = 1 if $scale == 0;
2021 # shortcut if already integer
2022 return $x if $scale == 1 && $dbt <= $dbd;
2023 # maximum digits before dot
2028 # not enough digits before dot, so round to zero
2031 elsif ( $scale == $dbd )
2038 $scale = $dbd - $scale;
2041 # pass sign to bround for rounding modes '+inf' and '-inf'
2042 my $m = Math::BigInt->new( $x->{sign} . $MBI->_str($x->{_m}));
2043 $m->bround($scale,$mode);
2044 $x->{_m} = $m->{value}; # get our mantissa back
2050 # accuracy: preserve $N digits, and overwrite the rest with 0's
2051 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2053 if (($_[0] || 0) < 0)
2055 require Carp; Carp::croak ('bround() needs positive accuracy');
2058 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
2059 return $x if !defined $scale; # no-op
2061 return $x if $x->modify('bround');
2063 # scale is now either $x->{_a}, $accuracy, or the user parameter
2064 # test whether $x already has lower accuracy, do nothing in this case
2065 # but do round if the accuracy is the same, since a math operation might
2066 # want to round a number with A=5 to 5 digits afterwards again
2067 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
2069 # scale < 0 makes no sense
2070 # never round a +-inf, NaN
2071 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
2073 # 1: $scale == 0 => keep all digits
2074 # 2: never round a 0
2075 # 3: if we should keep more digits than the mantissa has, do nothing
2076 if ($scale == 0 || $x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
2078 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
2082 # pass sign to bround for '+inf' and '-inf' rounding modes
2083 my $m = Math::BigInt->new( $x->{sign} . $MBI->_str($x->{_m}));
2085 $m->bround($scale,$mode); # round mantissa
2086 $x->{_m} = $m->{value}; # get our mantissa back
2087 $x->{_a} = $scale; # remember rounding
2088 delete $x->{_p}; # and clear P
2089 $x->bnorm(); # del trailing zeros gen. by bround()
2094 # return integer less or equal then $x
2095 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2097 return $x if $x->modify('bfloor');
2099 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2101 # if $x has digits after dot
2102 if ($x->{_es} eq '-')
2104 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2105 $x->{_e} = $MBI->_zero(); # trunc/norm
2106 $x->{_es} = '+'; # abs e
2107 $MBI->_inc($x->{_m}) if $x->{sign} eq '-'; # increment if negative
2109 $x->round($a,$p,$r);
2114 # return integer greater or equal then $x
2115 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2117 return $x if $x->modify('bceil');
2118 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2120 # if $x has digits after dot
2121 if ($x->{_es} eq '-')
2123 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2124 $x->{_e} = $MBI->_zero(); # trunc/norm
2125 $x->{_es} = '+'; # abs e
2126 $MBI->_inc($x->{_m}) if $x->{sign} eq '+'; # increment if positive
2128 $x->round($a,$p,$r);
2133 # shift right by $y (divide by power of $n)
2136 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2137 # objectify is costly, so avoid it
2138 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2140 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2143 return $x if $x->modify('brsft');
2144 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2146 $n = 2 if !defined $n; $n = $self->new($n);
2147 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
2152 # shift left by $y (multiply by power of $n)
2155 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2156 # objectify is costly, so avoid it
2157 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2159 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2162 return $x if $x->modify('blsft');
2163 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2165 $n = 2 if !defined $n; $n = $self->new($n);
2166 $x->bmul($n->bpow($y),$a,$p,$r,$y);
2169 ###############################################################################
2173 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
2178 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
2179 # or falling back to MBI::bxxx()
2180 my $name = $AUTOLOAD;
2182 $name =~ s/(.*):://; # split package
2183 my $c = $1 || $class;
2185 $c->import() if $IMPORT == 0;
2186 if (!method_alias($name))
2190 # delayed load of Carp and avoid recursion
2192 Carp::croak ("$c: Can't call a method without name");
2194 if (!method_hand_up($name))
2196 # delayed load of Carp and avoid recursion
2198 Carp::croak ("Can't call $c\-\>$name, not a valid method");
2200 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
2202 return &{"Math::BigInt"."::$name"}(@_);
2204 my $bname = $name; $bname =~ s/^f/b/;
2212 # return a copy of the exponent
2213 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2215 if ($x->{sign} !~ /^[+-]$/)
2217 my $s = $x->{sign}; $s =~ s/^[+-]//;
2218 return Math::BigInt->new($s); # -inf, +inf => +inf
2220 Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
2225 # return a copy of the mantissa
2226 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2228 if ($x->{sign} !~ /^[+-]$/)
2230 my $s = $x->{sign}; $s =~ s/^[+]//;
2231 return Math::BigInt->new($s); # -inf, +inf => +inf
2233 my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
2234 $m->bneg() if $x->{sign} eq '-';
2241 # return a copy of both the exponent and the mantissa
2242 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2244 if ($x->{sign} !~ /^[+-]$/)
2246 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2247 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2249 my $m = Math::BigInt->bzero();
2250 $m->{value} = $MBI->_copy($x->{_m});
2251 $m->bneg() if $x->{sign} eq '-';
2252 ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
2255 ##############################################################################
2256 # private stuff (internal use only)
2262 my $lib = ''; my @a;
2264 for ( my $i = 0; $i < $l ; $i++)
2266 if ( $_[$i] eq ':constant' )
2268 # This causes overlord er load to step in. 'binary' and 'integer'
2269 # are handled by BigInt.
2270 overload::constant float => sub { $self->new(shift); };
2272 elsif ($_[$i] eq 'upgrade')
2274 # this causes upgrading
2275 $upgrade = $_[$i+1]; # or undef to disable
2278 elsif ($_[$i] eq 'downgrade')
2280 # this causes downgrading
2281 $downgrade = $_[$i+1]; # or undef to disable
2284 elsif ($_[$i] eq 'lib')
2286 # alternative library
2287 $lib = $_[$i+1] || ''; # default Calc
2290 elsif ($_[$i] eq 'with')
2292 # alternative class for our private parts()
2293 # XXX: no longer supported
2294 # $MBI = $_[$i+1] || 'Math::BigInt';
2303 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2304 my $mbilib = eval { Math::BigInt->config()->{lib} };
2305 if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
2307 # MBI already loaded
2308 Math::BigInt->import('lib',"$lib,$mbilib", 'objectify');
2312 # MBI not loaded, or with ne "Math::BigInt::Calc"
2313 $lib .= ",$mbilib" if defined $mbilib;
2314 $lib =~ s/^,//; # don't leave empty
2315 # replacement library can handle lib statement, but also could ignore it
2318 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2319 # used in the same script, or eval inside import().
2320 require Math::BigInt;
2321 Math::BigInt->import( lib => $lib, 'objectify' );
2325 my $rc = "use Math::BigInt lib => '$lib', 'objectify';";
2331 require Carp; Carp::croak ("Couldn't load $lib: $! $@");
2333 $MBI = Math::BigInt->config()->{lib};
2335 # any non :constant stuff is handled by our parent, Exporter
2336 # even if @_ is empty, to give it a chance
2337 $self->SUPER::import(@a); # for subclasses
2338 $self->export_to_level(1,$self,@a); # need this, too
2343 # adjust m and e so that m is smallest possible
2344 # round number according to accuracy and precision settings
2345 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
2347 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2349 my $zeros = $MBI->_zeros($x->{_m}); # correct for trailing zeros
2352 my $z = $MBI->_new($zeros);
2353 $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
2354 if ($x->{_es} eq '-')
2356 if ($MBI->_acmp($x->{_e},$z) >= 0)
2358 $x->{_e} = $MBI->_sub ($x->{_e}, $z);
2359 $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
2363 $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
2369 $x->{_e} = $MBI->_add ($x->{_e}, $z);
2374 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
2375 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
2376 $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
2377 if $MBI->_is_zero($x->{_m});
2380 $x; # MBI bnorm is no-op, so dont call it
2383 ##############################################################################
2387 # return number as hexadecimal string (only for integers defined)
2388 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2390 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2391 return '0x0' if $x->is_zero();
2393 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2395 my $z = $MBI->_copy($x->{_m});
2396 if (! $MBI->_is_zero($x->{_e})) # > 0
2398 $MBI->_lsft( $z, $x->{_e},10);
2400 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2406 # return number as binary digit string (only for integers defined)
2407 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2409 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2410 return '0b0' if $x->is_zero();
2412 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2414 my $z = $MBI->_copy($x->{_m});
2415 if (! $MBI->_is_zero($x->{_e})) # > 0
2417 $MBI->_lsft( $z, $x->{_e},10);
2419 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2425 # return copy as a bigint representation of this BigFloat number
2426 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2428 my $z = $MBI->_copy($x->{_m});
2429 if ($x->{_es} eq '-') # < 0
2431 $MBI->_rsft( $z, $x->{_e},10);
2433 elsif (! $MBI->_is_zero($x->{_e})) # > 0
2435 $MBI->_lsft( $z, $x->{_e},10);
2437 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2444 my $class = ref($x) || $x;
2445 $x = $class->new(shift) unless ref($x);
2447 return 1 if $MBI->_is_zero($x->{_m});
2449 my $len = $MBI->_len($x->{_m});
2450 $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
2454 $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
2465 Math::BigFloat - Arbitrary size floating point math package
2472 $x = Math::BigFloat->new($str); # defaults to 0
2473 $nan = Math::BigFloat->bnan(); # create a NotANumber
2474 $zero = Math::BigFloat->bzero(); # create a +0
2475 $inf = Math::BigFloat->binf(); # create a +inf
2476 $inf = Math::BigFloat->binf('-'); # create a -inf
2477 $one = Math::BigFloat->bone(); # create a +1
2478 $one = Math::BigFloat->bone('-'); # create a -1
2481 $x->is_zero(); # true if arg is +0
2482 $x->is_nan(); # true if arg is NaN
2483 $x->is_one(); # true if arg is +1
2484 $x->is_one('-'); # true if arg is -1
2485 $x->is_odd(); # true if odd, false for even
2486 $x->is_even(); # true if even, false for odd
2487 $x->is_pos(); # true if >= 0
2488 $x->is_neg(); # true if < 0
2489 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2491 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2492 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2493 $x->sign(); # return the sign, either +,- or NaN
2494 $x->digit($n); # return the nth digit, counting from right
2495 $x->digit(-$n); # return the nth digit, counting from left
2497 # The following all modify their first argument. If you want to preserve
2498 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2499 # neccessary when mixing $a = $b assigments with non-overloaded math.
2502 $x->bzero(); # set $i to 0
2503 $x->bnan(); # set $i to NaN
2504 $x->bone(); # set $x to +1
2505 $x->bone('-'); # set $x to -1
2506 $x->binf(); # set $x to inf
2507 $x->binf('-'); # set $x to -inf
2509 $x->bneg(); # negation
2510 $x->babs(); # absolute value
2511 $x->bnorm(); # normalize (no-op)
2512 $x->bnot(); # two's complement (bit wise not)
2513 $x->binc(); # increment x by 1
2514 $x->bdec(); # decrement x by 1
2516 $x->badd($y); # addition (add $y to $x)
2517 $x->bsub($y); # subtraction (subtract $y from $x)
2518 $x->bmul($y); # multiplication (multiply $x by $y)
2519 $x->bdiv($y); # divide, set $x to quotient
2520 # return (quo,rem) or quo if scalar
2522 $x->bmod($y); # modulus ($x % $y)
2523 $x->bpow($y); # power of arguments ($x ** $y)
2524 $x->blsft($y); # left shift
2525 $x->brsft($y); # right shift
2526 # return (quo,rem) or quo if scalar
2528 $x->blog(); # logarithm of $x to base e (Euler's number)
2529 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
2531 $x->band($y); # bit-wise and
2532 $x->bior($y); # bit-wise inclusive or
2533 $x->bxor($y); # bit-wise exclusive or
2534 $x->bnot(); # bit-wise not (two's complement)
2536 $x->bsqrt(); # calculate square-root
2537 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
2538 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2540 $x->bround($N); # accuracy: preserve $N digits
2541 $x->bfround($N); # precision: round to the $Nth digit
2543 $x->bfloor(); # return integer less or equal than $x
2544 $x->bceil(); # return integer greater or equal than $x
2546 # The following do not modify their arguments:
2548 bgcd(@values); # greatest common divisor
2549 blcm(@values); # lowest common multiplicator
2551 $x->bstr(); # return string
2552 $x->bsstr(); # return string in scientific notation
2554 $x->as_int(); # return $x as BigInt
2555 $x->exponent(); # return exponent as BigInt
2556 $x->mantissa(); # return mantissa as BigInt
2557 $x->parts(); # return (mantissa,exponent) as BigInt
2559 $x->length(); # number of digits (w/o sign and '.')
2560 ($l,$f) = $x->length(); # number of digits, and length of fraction
2562 $x->precision(); # return P of $x (or global, if P of $x undef)
2563 $x->precision($n); # set P of $x to $n
2564 $x->accuracy(); # return A of $x (or global, if A of $x undef)
2565 $x->accuracy($n); # set A $x to $n
2567 # these get/set the appropriate global value for all BigFloat objects
2568 Math::BigFloat->precision(); # Precision
2569 Math::BigFloat->accuracy(); # Accuracy
2570 Math::BigFloat->round_mode(); # rounding mode
2574 All operators (inlcuding basic math operations) are overloaded if you
2575 declare your big floating point numbers as
2577 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2579 Operations with overloaded operators preserve the arguments, which is
2580 exactly what you expect.
2582 =head2 Canonical notation
2584 Input to these routines are either BigFloat objects, or strings of the
2585 following four forms:
2599 C</^[+-]\d+E[+-]?\d+$/>
2603 C</^[+-]\d*\.\d+E[+-]?\d+$/>
2607 all with optional leading and trailing zeros and/or spaces. Additonally,
2608 numbers are allowed to have an underscore between any two digits.
2610 Empty strings as well as other illegal numbers results in 'NaN'.
2612 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2613 are always stored in normalized form. On a string, it creates a BigFloat
2618 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2620 The string output will always have leading and trailing zeros stripped and drop
2621 a plus sign. C<bstr()> will give you always the form with a decimal point,
2622 while C<bsstr()> (s for scientific) gives you the scientific notation.
2624 Input bstr() bsstr()
2626 ' -123 123 123' '-123123123' '-123123123E0'
2627 '00.0123' '0.0123' '123E-4'
2628 '123.45E-2' '1.2345' '12345E-4'
2629 '10E+3' '10000' '1E4'
2631 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2632 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2633 return either undef, <0, 0 or >0 and are suited for sort.
2635 Actual math is done by using the class defined with C<with => Class;> (which
2636 defaults to BigInts) to represent the mantissa and exponent.
2638 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2639 represent the result when input arguments are not numbers, as well as
2640 the result of dividing by zero.
2642 =head2 C<mantissa()>, C<exponent()> and C<parts()>
2644 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2645 as BigInts such that:
2647 $m = $x->mantissa();
2648 $e = $x->exponent();
2649 $y = $m * ( 10 ** $e );
2650 print "ok\n" if $x == $y;
2652 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2654 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2656 Currently the mantissa is reduced as much as possible, favouring higher
2657 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2658 This might change in the future, so do not depend on it.
2660 =head2 Accuracy vs. Precision
2662 See also: L<Rounding|Rounding>.
2664 Math::BigFloat supports both precision and accuracy. For a full documentation,
2665 examples and tips on these topics please see the large section in
2668 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2669 a operation consumes all resources, each operation produces no more than
2670 the requested number of digits.
2672 Please refer to BigInt's documentation for the precedence rules of which
2673 accuracy/precision setting will be used.
2675 If there is no gloabl precision set, B<and> the operation inquestion was not
2676 called with a requested precision or accuracy, B<and> the input $x has no
2677 accuracy or precision set, then a fallback parameter will be used. For
2678 historical reasons, it is called C<div_scale> and can be accessed via:
2680 $d = Math::BigFloat->div_scale(); # query
2681 Math::BigFloat->div_scale($n); # set to $n digits
2683 The default value is 40 digits.
2685 In case the result of one operation has more precision than specified,
2686 it is rounded. The rounding mode taken is either the default mode, or the one
2687 supplied to the operation after the I<scale>:
2689 $x = Math::BigFloat->new(2);
2690 Math::BigFloat->precision(5); # 5 digits max
2691 $y = $x->copy()->bdiv(3); # will give 0.66666
2692 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2693 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
2694 Math::BigFloat->round_mode('zero');
2695 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2701 =item ffround ( +$scale )
2703 Rounds to the $scale'th place left from the '.', counting from the dot.
2704 The first digit is numbered 1.
2706 =item ffround ( -$scale )
2708 Rounds to the $scale'th place right from the '.', counting from the dot.
2712 Rounds to an integer.
2714 =item fround ( +$scale )
2716 Preserves accuracy to $scale digits from the left (aka significant digits)
2717 and pads the rest with zeros. If the number is between 1 and -1, the
2718 significant digits count from the first non-zero after the '.'
2720 =item fround ( -$scale ) and fround ( 0 )
2722 These are effectively no-ops.
2726 All rounding functions take as a second parameter a rounding mode from one of
2727 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2729 The default rounding mode is 'even'. By using
2730 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
2731 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2732 no longer supported.
2733 The second parameter to the round functions then overrides the default
2736 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
2737 'trunc' as rounding mode to make it equivalent to:
2742 You can override this by passing the desired rounding mode as parameter to
2745 $x = Math::BigFloat->new(2.5);
2746 $y = $x->as_number('odd'); # $y = 3
2752 =head1 Autocreating constants
2754 After C<use Math::BigFloat ':constant'> all the floating point constants
2755 in the given scope are converted to C<Math::BigFloat>. This conversion
2756 happens at compile time.
2760 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2762 prints the value of C<2E-100>. Note that without conversion of
2763 constants the expression 2E-100 will be calculated as normal floating point
2766 Please note that ':constant' does not affect integer constants, nor binary
2767 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2772 Math with the numbers is done (by default) by a module called
2773 Math::BigInt::Calc. This is equivalent to saying:
2775 use Math::BigFloat lib => 'Calc';
2777 You can change this by using:
2779 use Math::BigFloat lib => 'BitVect';
2781 The following would first try to find Math::BigInt::Foo, then
2782 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2784 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2786 Calc.pm uses as internal format an array of elements of some decimal base
2787 (usually 1e7, but this might be differen for some systems) with the least
2788 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2789 significant bit first. Other modules might use even different means of
2790 representing the numbers. See the respective module documentation for further
2793 Please note that Math::BigFloat does B<not> use the denoted library itself,
2794 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2797 use Math::BigInt lib => 'GMP';
2800 you can roll it all into one line:
2802 use Math::BigFloat lib => 'GMP';
2804 It is also possible to just require Math::BigFloat:
2806 require Math::BigFloat;
2808 This will load the neccessary things (like BigInt) when they are needed, and
2811 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2812 you ever wanted to know about loading a different library.
2814 =head2 Using Math::BigInt::Lite
2816 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2819 use Math::BigFloat with => 'Math::BigInt::Lite';
2821 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2822 can combine these if you want. For instance, you may want to use
2823 Math::BigInt objects in your main script, too.
2827 use Math::BigFloat with => 'Math::BigInt::Lite';
2829 Of course, you can combine this with the C<lib> parameter.
2832 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2834 There is no need for a "use Math::BigInt;" statement, even if you want to
2835 use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2836 always loads it. But if you add it, add it B<before>:
2840 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2842 Notice that the module with the last C<lib> will "win" and thus
2843 it's lib will be used if the lib is available:
2846 use Math::BigInt lib => 'Bar,Baz';
2847 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2849 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2850 words, Math::BigFloat will try to retain previously loaded libs when you
2851 don't specify it onem but if you specify one, it will try to load them.
2853 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2854 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2855 same as trying the latter load alone, except for the fact that one of Bar or
2856 Baz might be loaded needlessly in an intermidiate step (and thus hang around
2857 and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2858 will still be tried to be loaded, but this is not as time/memory consuming as
2859 actually loading one of them. Still, this type of usage is not recommended due
2862 The old way (loading the lib only in BigInt) still works though:
2865 use Math::BigInt lib => 'Bar,Baz';
2868 You can even load Math::BigInt afterwards:
2872 use Math::BigInt lib => 'Bar,Baz';
2874 But this has the same problems like #5, it will first load Calc
2875 (Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2876 Baz, depending on which of them works and is usable/loadable. Since this
2877 loads Calc unnecc., it is not recommended.
2879 Since it also possible to just require Math::BigFloat, this poses the question
2880 about what libary this will use:
2882 require Math::BigFloat;
2883 my $x = Math::BigFloat->new(123); $x += 123;
2885 It will use Calc. Please note that the call to import() is still done, but
2886 only when you use for the first time some Math::BigFloat math (it is triggered
2887 via any constructor, so the first time you create a Math::BigFloat, the load
2888 will happen in the background). This means:
2890 require Math::BigFloat;
2891 Math::BigFloat->import ( lib => 'Foo,Bar' );
2893 would be the same as:
2895 use Math::BigFloat lib => 'Foo, Bar';
2897 But don't try to be clever to insert some operations in between:
2899 require Math::BigFloat;
2900 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
2901 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
2902 $x = Math::BigFloat->bone()+4; # now use Pari
2904 While this works, it loads Calc needlessly. But maybe you just wanted that?
2906 B<Examples #3 is highly recommended> for daily usage.
2910 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
2916 =item stringify, bstr()
2918 Both stringify and bstr() now drop the leading '+'. The old code would return
2919 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2920 reasoning and details.
2924 The following will probably not do what you expect:
2926 print $c->bdiv(123.456),"\n";
2928 It prints both quotient and reminder since print works in list context. Also,
2929 bdiv() will modify $c, so be carefull. You probably want to use
2931 print $c / 123.456,"\n";
2932 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2936 =item Modifying and =
2940 $x = Math::BigFloat->new(5);
2943 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2944 a second reference to the B<same> object and stores it in $y. Thus anything
2945 that modifies $x will modify $y (except overloaded math operators), and vice
2946 versa. See L<Math::BigInt> for details and how to avoid that.
2950 C<bpow()> now modifies the first argument, unlike the old code which left
2951 it alone and only returned the result. This is to be consistent with
2952 C<badd()> etc. The first will modify $x, the second one won't:
2954 print bpow($x,$i),"\n"; # modify $x
2955 print $x->bpow($i),"\n"; # ditto
2956 print $x ** $i,"\n"; # leave $x alone
2962 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
2963 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
2965 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
2966 because they solve the autoupgrading/downgrading issue, at least partly.
2969 L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
2970 more documentation including a full version history, testcases, empty
2971 subclass files and benchmarks.
2975 This program is free software; you may redistribute it and/or modify it under
2976 the same terms as Perl itself.
2980 Mark Biggar, overloaded interface by Ilya Zakharevich.
2981 Completely rewritten by Tels http://bloodgate.com in 2001, 2002, and still