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 { $_[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 fneg
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 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();
137 # handle '+inf', '-inf' first
138 if ($wanted =~ /^[+-]?inf$/)
140 return $downgrade->new($wanted) if $downgrade;
142 $self->{_e} = $MBI->_zero();
144 $self->{_m} = $MBI->_zero();
145 $self->{sign} = $wanted;
146 $self->{sign} = '+inf' if $self->{sign} eq 'inf';
147 return $self->bnorm();
150 # shortcut for simple forms like '12' that neither have trailing nor leading
152 if ($wanted =~ /^([+-]?)([1-9][0-9]*[1-9])$/)
154 $self->{_e} = $MBI->_zero();
156 $self->{sign} = $1 || '+';
157 $self->{_m} = $MBI->_new($2);
158 return $self->round(@r) if !$downgrade;
161 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
167 Carp::croak ("$wanted is not a number initialized to $class");
170 return $downgrade->bnan() if $downgrade;
172 $self->{_e} = $MBI->_zero();
174 $self->{_m} = $MBI->_zero();
175 $self->{sign} = $nan;
179 # make integer from mantissa by adjusting exp, then convert to int
180 $self->{_e} = $MBI->_new($$ev); # exponent
181 $self->{_es} = $$es || '+';
182 my $mantissa = "$$miv$$mfv"; # create mant.
183 $mantissa =~ s/^0+(\d)/$1/; # strip leading zeros
184 $self->{_m} = $MBI->_new($mantissa); # create mant.
186 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
187 if (CORE::length($$mfv) != 0)
189 my $len = $MBI->_new( CORE::length($$mfv));
190 ($self->{_e}, $self->{_es}) =
191 _e_sub ($self->{_e}, $len, $self->{_es}, '+');
193 # we can only have trailing zeros on the mantissa if $$mfv eq ''
196 # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
197 # because that is faster, especially when _m is not stored in base 10.
198 my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/;
201 my $z = $MBI->_new($zeros);
202 # turn '120e2' into '12e3'
203 $MBI->_rsft ( $self->{_m}, $z, 10);
204 _e_add ( $self->{_e}, $z, $self->{_es}, '+');
207 $self->{sign} = $$mis;
209 # for something like 0Ey, set y to 1, and -0 => +0
210 # Check $$miv for beeing '0' and $$mfv eq '', because otherwise _m could not
211 # have become 0. That's faster than to call $MBI->_is_zero().
212 $self->{sign} = '+', $self->{_e} = $MBI->_one()
213 if $$miv eq '0' and $$mfv eq '';
215 return $self->round(@r) if !$downgrade;
217 # if downgrade, inf, NaN or integers go down
219 if ($downgrade && $self->{_es} eq '+')
221 if ($MBI->_is_zero( $self->{_e} ))
223 return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
225 return $downgrade->new($self->bsstr());
227 $self->bnorm()->round(@r); # first normalize, then round
235 # if two arguments, the first one is the class to "swallow" subclasses
243 return unless ref($x); # only for objects
245 my $self = {}; bless $self,$c;
247 $self->{sign} = $x->{sign};
248 $self->{_es} = $x->{_es};
249 $self->{_m} = $MBI->_copy($x->{_m});
250 $self->{_e} = $MBI->_copy($x->{_e});
251 $self->{_a} = $x->{_a} if defined $x->{_a};
252 $self->{_p} = $x->{_p} if defined $x->{_p};
258 # used by parent class bone() to initialize number to NaN
264 my $class = ref($self);
265 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
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 +-inf
282 my $class = ref($self);
283 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
286 $IMPORT=1; # call our import only once
287 $self->{_m} = $MBI->_zero();
288 $self->{_e} = $MBI->_zero();
294 # used by parent class bone() to initialize number to 1
296 $IMPORT=1; # call our import only once
297 $self->{_m} = $MBI->_one();
298 $self->{_e} = $MBI->_zero();
304 # used by parent class bone() to initialize number to 0
306 $IMPORT=1; # call our import only once
307 $self->{_m} = $MBI->_zero();
308 $self->{_e} = $MBI->_one();
314 my ($self,$class) = @_;
315 return if $class =~ /^Math::BigInt/; # we aren't one of these
316 UNIVERSAL::isa($self,$class);
321 # return (later set?) configuration data as hash ref
322 my $class = shift || 'Math::BigFloat';
324 my $cfg = $class->SUPER::config(@_);
326 # now we need only to override the ones that are different from our parent
327 $cfg->{class} = $class;
332 ##############################################################################
333 # string conversation
337 # (ref to BFLOAT or num_str ) return num_str
338 # Convert number from internal format to (non-scientific) string format.
339 # internal format is always normalized (no leading zeros, "-0" => "+0")
340 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
342 if ($x->{sign} !~ /^[+-]$/)
344 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
348 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
351 my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
354 $es = $MBI->_str($x->{_m});
355 $len = CORE::length($es);
356 my $e = $MBI->_num($x->{_e});
357 $e = -$e if $x->{_es} eq '-';
361 # if _e is bigger than a scalar, the following will blow your memory
364 my $r = abs($e) - $len;
365 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
369 substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
370 $cad = -$cad if $x->{_es} eq '-';
376 $es .= '0' x $e; $len += $e; $cad = 0;
380 $es = '-'.$es if $x->{sign} eq '-';
381 # if set accuracy or precision, pad with zeros on the right side
382 if ((defined $x->{_a}) && ($not_zero))
384 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
385 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
386 $zeros = $x->{_a} - $len if $cad != $len;
387 $es .= $dot.'0' x $zeros if $zeros > 0;
389 elsif ((($x->{_p} || 0) < 0))
391 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
392 my $zeros = -$x->{_p} + $cad;
393 $es .= $dot.'0' x $zeros if $zeros > 0;
400 # (ref to BFLOAT or num_str ) return num_str
401 # Convert number from internal format to scientific string format.
402 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
403 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
405 if ($x->{sign} !~ /^[+-]$/)
407 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
410 my $sep = 'e'.$x->{_es};
411 my $sign = $x->{sign}; $sign = '' if $sign eq '+';
412 $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
417 # Make a number from a BigFloat object
418 # simple return a string and let Perl's atoi()/atof() handle the rest
419 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
423 ##############################################################################
424 # public stuff (usually prefixed with "b")
428 # (BINT or num_str) return BINT
429 # negate number or make a negated number from string
430 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
432 return $x if $x->modify('bneg');
434 # for +0 dont negate (to have always normalized +0). Does nothing for 'NaN'
435 $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
440 # XXX TODO this must be overwritten and return NaN for non-integer values
441 # band(), bior(), bxor(), too
444 # $class->SUPER::bnot($class,@_);
449 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
452 my ($self,$x,$y) = (ref($_[0]),@_);
453 # objectify is costly, so avoid it
454 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
456 ($self,$x,$y) = objectify(2,@_);
459 return $upgrade->bcmp($x,$y) if defined $upgrade &&
460 ((!$x->isa($self)) || (!$y->isa($self)));
462 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
464 # handle +-inf and NaN
465 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
466 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
467 return +1 if $x->{sign} eq '+inf';
468 return -1 if $x->{sign} eq '-inf';
469 return -1 if $y->{sign} eq '+inf';
473 # check sign for speed first
474 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
475 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
478 my $xz = $x->is_zero();
479 my $yz = $y->is_zero();
480 return 0 if $xz && $yz; # 0 <=> 0
481 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
482 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
484 # adjust so that exponents are equal
485 my $lxm = $MBI->_len($x->{_m});
486 my $lym = $MBI->_len($y->{_m});
487 # the numify somewhat limits our length, but makes it much faster
488 my ($xes,$yes) = (1,1);
489 $xes = -1 if $x->{_es} ne '+';
490 $yes = -1 if $y->{_es} ne '+';
491 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
492 my $ly = $lym + $yes * $MBI->_num($y->{_e});
493 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
494 return $l <=> 0 if $l != 0;
496 # lengths (corrected by exponent) are equal
497 # so make mantissa equal length by padding with zero (shift left)
498 my $diff = $lxm - $lym;
499 my $xm = $x->{_m}; # not yet copy it
503 $ym = $MBI->_copy($y->{_m});
504 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
508 $xm = $MBI->_copy($x->{_m});
509 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
511 my $rc = $MBI->_acmp($xm,$ym);
512 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
518 # Compares 2 values, ignoring their signs.
519 # Returns one of undef, <0, =0, >0. (suitable for sort)
522 my ($self,$x,$y) = (ref($_[0]),@_);
523 # objectify is costly, so avoid it
524 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
526 ($self,$x,$y) = objectify(2,@_);
529 return $upgrade->bacmp($x,$y) if defined $upgrade &&
530 ((!$x->isa($self)) || (!$y->isa($self)));
532 # handle +-inf and NaN's
533 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
535 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
536 return 0 if ($x->is_inf() && $y->is_inf());
537 return 1 if ($x->is_inf() && !$y->is_inf());
542 my $xz = $x->is_zero();
543 my $yz = $y->is_zero();
544 return 0 if $xz && $yz; # 0 <=> 0
545 return -1 if $xz && !$yz; # 0 <=> +y
546 return 1 if $yz && !$xz; # +x <=> 0
548 # adjust so that exponents are equal
549 my $lxm = $MBI->_len($x->{_m});
550 my $lym = $MBI->_len($y->{_m});
551 my ($xes,$yes) = (1,1);
552 $xes = -1 if $x->{_es} ne '+';
553 $yes = -1 if $y->{_es} ne '+';
554 # the numify somewhat limits our length, but makes it much faster
555 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
556 my $ly = $lym + $yes * $MBI->_num($y->{_e});
558 return $l <=> 0 if $l != 0;
560 # lengths (corrected by exponent) are equal
561 # so make mantissa equal-length by padding with zero (shift left)
562 my $diff = $lxm - $lym;
563 my $xm = $x->{_m}; # not yet copy it
567 $ym = $MBI->_copy($y->{_m});
568 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
572 $xm = $MBI->_copy($x->{_m});
573 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
575 $MBI->_acmp($xm,$ym);
580 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
581 # return result as BFLOAT
584 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
585 # objectify is costly, so avoid it
586 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
588 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
591 # inf and NaN handling
592 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
595 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
597 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
599 # +inf++inf or -inf+-inf => same, rest is NaN
600 return $x if $x->{sign} eq $y->{sign};
603 # +-inf + something => +inf; something +-inf => +-inf
604 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
608 return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
609 ((!$x->isa($self)) || (!$y->isa($self)));
611 # speed: no add for 0+y or x+0
612 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
613 if ($x->is_zero()) # 0+y
615 # make copy, clobbering up x (modify in place!)
616 $x->{_e} = $MBI->_copy($y->{_e});
617 $x->{_es} = $y->{_es};
618 $x->{_m} = $MBI->_copy($y->{_m});
619 $x->{sign} = $y->{sign} || $nan;
620 return $x->round($a,$p,$r,$y);
623 # take lower of the two e's and adapt m1 to it to match m2
625 $e = $MBI->_zero() if !defined $e; # if no BFLOAT?
626 $e = $MBI->_copy($e); # make copy (didn't do it yet)
630 ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
632 my $add = $MBI->_copy($y->{_m});
634 if ($es eq '-') # < 0
636 $MBI->_lsft( $x->{_m}, $e, 10);
637 ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
639 elsif (!$MBI->_is_zero($e)) # > 0
641 $MBI->_lsft($add, $e, 10);
643 # else: both e are the same, so just leave them
645 if ($x->{sign} eq $y->{sign})
648 $x->{_m} = $MBI->_add($x->{_m}, $add);
652 ($x->{_m}, $x->{sign}) =
653 _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
656 # delete trailing zeros, then round
657 $x->bnorm()->round($a,$p,$r,$y);
660 # sub bsub is inherited from Math::BigInt!
664 # increment arg by one
665 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
667 if ($x->{_es} eq '-')
669 return $x->badd($self->bone(),@r); # digits after dot
672 if (!$MBI->_is_zero($x->{_e})) # _e == 0 for NaN, inf, -inf
674 # 1e2 => 100, so after the shift below _m has a '0' as last digit
675 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
676 $x->{_e} = $MBI->_zero(); # normalize
678 # we know that the last digit of $x will be '1' or '9', depending on the
682 if ($x->{sign} eq '+')
684 $MBI->_inc($x->{_m});
685 return $x->bnorm()->bround(@r);
687 elsif ($x->{sign} eq '-')
689 $MBI->_dec($x->{_m});
690 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
691 return $x->bnorm()->bround(@r);
693 # inf, nan handling etc
694 $x->badd($self->bone(),@r); # badd() does round
699 # decrement arg by one
700 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
702 if ($x->{_es} eq '-')
704 return $x->badd($self->bone('-'),@r); # digits after dot
707 if (!$MBI->_is_zero($x->{_e}))
709 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
710 $x->{_e} = $MBI->_zero(); # normalize
714 my $zero = $x->is_zero();
716 if (($x->{sign} eq '-') || $zero)
718 $MBI->_inc($x->{_m});
719 $x->{sign} = '-' if $zero; # 0 => 1 => -1
720 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
721 return $x->bnorm()->round(@r);
724 elsif ($x->{sign} eq '+')
726 $MBI->_dec($x->{_m});
727 return $x->bnorm()->round(@r);
729 # inf, nan handling etc
730 $x->badd($self->bone('-'),@r); # does round
737 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
739 # $base > 0, $base != 1; if $base == undef default to $base == e
742 # we need to limit the accuracy to protect against overflow
745 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
747 # also takes care of the "error in _find_round_parameters?" case
748 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
751 # no rounding at all, so must use fallback
752 if (scalar @params == 0)
754 # simulate old behaviour
755 $params[0] = $self->div_scale(); # and round to it as accuracy
756 $params[1] = undef; # P = undef
757 $scale = $params[0]+4; # at least four more for proper round
758 $params[2] = $r; # round mode by caller or undef
759 $fallback = 1; # to clear a/p afterwards
763 # the 4 below is empirical, and there might be cases where it is not
765 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
768 return $x->bzero(@params) if $x->is_one();
769 # base not defined => base == Euler's constant e
772 # make object, since we don't feed it through objectify() to still get the
773 # case of $base == undef
774 $base = $self->new($base) unless ref($base);
775 # $base > 0; $base != 1
776 return $x->bnan() if $base->is_zero() || $base->is_one() ||
777 $base->{sign} ne '+';
778 # if $x == $base, we know the result must be 1.0
779 if ($x->bcmp($base) == 0)
781 $x->bone('+',@params);
784 # clear a/p after round, since user did not request it
785 delete $x->{_a}; delete $x->{_p};
791 # when user set globals, they would interfere with our calculation, so
792 # disable them and later re-enable them
794 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
795 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
796 # we also need to disable any set A or P on $x (_find_round_parameters took
797 # them already into account), since these would interfere, too
798 delete $x->{_a}; delete $x->{_p};
799 # need to disable $upgrade in BigInt, to avoid deep recursion
800 local $Math::BigInt::upgrade = undef;
801 local $Math::BigFloat::downgrade = undef;
803 # upgrade $x if $x is not a BigFloat (handle BigInt input)
804 if (!$x->isa('Math::BigFloat'))
806 $x = Math::BigFloat->new($x);
812 # If the base is defined and an integer, try to calculate integer result
813 # first. This is very fast, and in case the real result was found, we can
815 if (defined $base && $base->is_int() && $x->is_int())
817 my $i = $MBI->_copy( $x->{_m} );
818 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
819 my $int = Math::BigInt->bzero();
821 $int->blog($base->as_number());
823 if ($base->as_number()->bpow($int) == $x)
825 # found result, return it
826 $x->{_m} = $int->{value};
827 $x->{_e} = $MBI->_zero();
836 # first calculate the log to base e (using reduction by 10 (and probably 2))
837 $self->_log_10($x,$scale);
839 # and if a different base was requested, convert it
842 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
843 # not ln, but some other base (don't modify $base)
844 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
848 # shortcut to not run through _find_round_parameters again
849 if (defined $params[0])
851 $x->bround($params[0],$params[2]); # then round accordingly
855 $x->bfround($params[1],$params[2]); # then round accordingly
859 # clear a/p after round, since user did not request it
860 delete $x->{_a}; delete $x->{_p};
863 $$abr = $ab; $$pbr = $pb;
870 # internal log function to calculate ln() based on Taylor series.
871 # Modifies $x in place.
872 my ($self,$x,$scale) = @_;
874 # in case of $x == 1, result is 0
875 return $x->bzero() if $x->is_one();
877 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
881 # Taylor: | u 1 u^3 1 u^5 |
882 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
883 # |_ v 3 v^3 5 v^5 _|
885 # This takes much more steps to calculate the result and is thus not used
888 # Taylor: | u 1 u^2 1 u^3 |
889 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
890 # |_ x 2 x^2 3 x^3 _|
892 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
894 $v = $x->copy(); $v->binc(); # v = x+1
895 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
896 $x->bdiv($v,$scale); # first term: u/v
899 $u *= $u; $v *= $v; # u^2, v^2
900 $below->bmul($v); # u^3, v^3
902 $factor = $self->new(3); $f = $self->new(2);
904 my $steps = 0 if DEBUG;
905 $limit = $self->new("1E-". ($scale-1));
908 # we calculate the next term, and add it to the last
909 # when the next term is below our limit, it won't affect the outcome
910 # anymore, so we stop
912 # calculating the next term simple from over/below will result in quite
913 # a time hog if the input has many digits, since over and below will
914 # accumulate more and more digits, and the result will also have many
915 # digits, but in the end it is rounded to $scale digits anyway. So if we
916 # round $over and $below first, we save a lot of time for the division
917 # (not with log(1.2345), but try log (123**123) to see what I mean. This
918 # can introduce a rounding error if the division result would be f.i.
919 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
920 # if we truncated $over and $below we might get 0.12345. Does this matter
921 # for the end result? So we give $over and $below 4 more digits to be
922 # on the safe side (unscientific error handling as usual... :+D
924 $next = $over->copy->bround($scale+4)->bdiv(
925 $below->copy->bmul($factor)->bround($scale+4),
929 ## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
931 last if $next->bacmp($limit) <= 0;
933 delete $next->{_a}; delete $next->{_p};
935 # calculate things for the next term
936 $over *= $u; $below *= $v; $factor->badd($f);
939 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
942 $x->bmul($f); # $x *= 2
943 print "took $steps steps\n" if DEBUG;
948 # Internal log function based on reducing input to the range of 0.1 .. 9.99
949 # and then "correcting" the result to the proper one. Modifies $x in place.
950 my ($self,$x,$scale) = @_;
952 # taking blog() from numbers greater than 10 takes a *very long* time, so we
953 # break the computation down into parts based on the observation that:
954 # blog(x*y) = blog(x) + blog(y)
955 # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
956 # the faster it get's, especially because 2*$x takes about 10 times as long,
957 # so by dividing $x by 10 we make it at least factor 100 faster...)
959 # The same observation is valid for numbers smaller than 0.1 (e.g. computing
960 # log(1) is fastest, and the farther away we get from 1, the longer it takes)
961 # so we also 'break' this down by multiplying $x with 10 and subtract the
962 # log(10) afterwards to get the correct result.
964 # calculate nr of digits before dot
965 my $dbd = $MBI->_num($x->{_e});
966 $dbd = -$dbd if $x->{_es} eq '-';
967 $dbd += $MBI->_len($x->{_m});
969 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
972 my $calc = 1; # do some calculation?
974 # disable the shortcut for 10, since we need log(10) and this would recurse
976 if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
978 $dbd = 0; # disable shortcut
979 # we can use the cached value in these cases
980 if ($scale <= $LOG_10_A)
982 $x->bzero(); $x->badd($LOG_10);
983 $calc = 0; # no need to calc, but round
988 # disable the shortcut for 2, since we maybe have it cached
989 if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
991 $dbd = 0; # disable shortcut
992 # we can use the cached value in these cases
993 if ($scale <= $LOG_2_A)
995 $x->bzero(); $x->badd($LOG_2);
996 $calc = 0; # no need to calc, but round
1001 # if $x = 0.1, we know the result must be 0-log(10)
1002 if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
1003 $MBI->_is_one($x->{_m}))
1005 $dbd = 0; # disable shortcut
1006 # we can use the cached value in these cases
1007 if ($scale <= $LOG_10_A)
1009 $x->bzero(); $x->bsub($LOG_10);
1010 $calc = 0; # no need to calc, but round
1014 return if $calc == 0; # already have the result
1016 # default: these correction factors are undef and thus not used
1017 my $l_10; # value of ln(10) to A of $scale
1018 my $l_2; # value of ln(2) to A of $scale
1020 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1021 # so don't do this shortcut for 1 or 0
1022 if (($dbd > 1) || ($dbd < 0))
1024 # convert our cached value to an object if not already (avoid doing this
1025 # at import() time, since not everybody needs this)
1026 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1028 #print "x = $x, dbd = $dbd, calc = $calc\n";
1029 # got more than one digit before the dot, or more than one zero after the
1031 # log(123) == log(1.23) + log(10) * 2
1032 # log(0.0123) == log(1.23) - log(10) * 2
1034 if ($scale <= $LOG_10_A)
1037 $l_10 = $LOG_10->copy(); # copy for mul
1041 # else: slower, compute it (but don't cache it, because it could be big)
1042 # also disable downgrade for this code path
1043 local $Math::BigFloat::downgrade = undef;
1044 $l_10 = $self->new(10)->blog(undef,$scale); # scale+4, actually
1046 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
1047 $l_10->bmul( $self->new($dbd)); # log(10) * (digits_before_dot-1)
1054 ($x->{_e}, $x->{_es}) =
1055 _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
1059 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1061 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1062 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1064 $HALF = $self->new($HALF) unless ref($HALF);
1066 my $twos = 0; # default: none (0 times)
1067 my $two = $self->new(2);
1068 while ($x->bacmp($HALF) <= 0)
1070 $twos--; $x->bmul($two);
1072 while ($x->bacmp($two) >= 0)
1074 $twos++; $x->bdiv($two,$scale+4); # keep all digits
1076 # $twos > 0 => did mul 2, < 0 => did div 2 (never both)
1077 # calculate correction factor based on ln(2)
1080 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1081 if ($scale <= $LOG_2_A)
1084 $l_2 = $LOG_2->copy(); # copy for mul
1088 # else: slower, compute it (but don't cache it, because it could be big)
1089 # also disable downgrade for this code path
1090 local $Math::BigFloat::downgrade = undef;
1091 $l_2 = $two->blog(undef,$scale); # scale+4, actually
1093 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
1096 $self->_log($x,$scale); # need to do the "normal" way
1097 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
1098 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
1099 # all done, $x contains now the result
1104 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1105 # does not modify arguments, but returns new object
1106 # Lowest Common Multiplicator
1108 my ($self,@arg) = objectify(0,@_);
1109 my $x = $self->new(shift @arg);
1110 while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); }
1116 # (BINT or num_str, BINT or num_str) return BINT
1117 # does not modify arguments, but returns new object
1120 $y = __PACKAGE__->new($y) if !ref($y);
1122 my $x = $y->copy()->babs(); # keep arguments
1124 return $x->bnan() if $x->{sign} !~ /^[+-]$/ # x NaN?
1125 || !$x->is_int(); # only for integers now
1129 my $t = shift; $t = $self->new($t) if !ref($t);
1130 $y = $t->copy()->babs();
1132 return $x->bnan() if $y->{sign} !~ /^[+-]$/ # y NaN?
1133 || !$y->is_int(); # only for integers now
1135 # greatest common divisor
1136 while (! $y->is_zero())
1138 ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1141 last if $x->is_one();
1146 ##############################################################################
1150 # Internal helper sub to take two positive integers and their signs and
1151 # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1152 # output ($CALC,('+'|'-'))
1153 my ($x,$y,$xs,$ys) = @_;
1155 # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1158 $x = $MBI->_add ($x, $y ); # a+b
1159 # the sign follows $xs
1163 my $a = $MBI->_acmp($x,$y);
1166 $x = $MBI->_sub ($x , $y); # abs sub
1170 $x = $MBI->_zero(); # result is 0
1175 $x = $MBI->_sub ( $y, $x, 1 ); # abs sub
1183 # Internal helper sub to take two positive integers and their signs and
1184 # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1185 # output ($CALC,('+'|'-'))
1186 my ($x,$y,$xs,$ys) = @_;
1190 _e_add($x,$y,$xs,$ys); # call add (does subtract now)
1193 ###############################################################################
1194 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1198 # return true if arg (BFLOAT or num_str) is an integer
1199 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1201 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1202 $x->{_es} eq '+'; # 1e-1 => no integer
1208 # return true if arg (BFLOAT or num_str) is zero
1209 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1211 return 1 if $x->{sign} eq '+' && $MBI->_is_zero($x->{_m});
1217 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1218 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1220 $sign = '+' if !defined $sign || $sign ne '-';
1222 if ($x->{sign} eq $sign &&
1223 $MBI->_is_zero($x->{_e}) && $MBI->_is_one($x->{_m}));
1229 # return true if arg (BFLOAT or num_str) is odd or false if even
1230 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1232 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1233 ($MBI->_is_zero($x->{_e}) && $MBI->_is_odd($x->{_m}));
1239 # return true if arg (BINT or num_str) is even or false if odd
1240 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1242 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
1243 return 1 if ($x->{_es} eq '+' # 123.45 is never
1244 && $MBI->_is_even($x->{_m})); # but 1200 is
1250 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
1251 # (BINT or num_str, BINT or num_str) return BINT
1254 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1255 # objectify is costly, so avoid it
1256 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1258 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1261 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1264 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1266 return $x->bnan() if $x->is_zero() || $y->is_zero();
1267 # result will always be +-inf:
1268 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1269 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1270 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1271 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1272 return $x->binf('-');
1275 return $x->bzero() if $x->is_zero() || $y->is_zero();
1277 return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
1278 ((!$x->isa($self)) || (!$y->isa($self)));
1280 # aEb * cEd = (a*c)E(b+d)
1281 $MBI->_mul($x->{_m},$y->{_m});
1282 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1285 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1286 return $x->bnorm()->round($a,$p,$r,$y);
1291 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1292 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1295 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1296 # objectify is costly, so avoid it
1297 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1299 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1302 return $self->_div_inf($x,$y)
1303 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1305 # x== 0 # also: or y == 1 or y == -1
1306 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1309 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1311 # we need to limit the accuracy to protect against overflow
1313 my (@params,$scale);
1314 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1316 return $x if $x->is_nan(); # error in _find_round_parameters?
1318 # no rounding at all, so must use fallback
1319 if (scalar @params == 0)
1321 # simulate old behaviour
1322 $params[0] = $self->div_scale(); # and round to it as accuracy
1323 $scale = $params[0]+4; # at least four more for proper round
1324 $params[2] = $r; # round mode by caller or undef
1325 $fallback = 1; # to clear a/p afterwards
1329 # the 4 below is empirical, and there might be cases where it is not
1331 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1334 my $rem; $rem = $self->bzero() if wantarray;
1336 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1338 my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
1339 $scale = $lx if $lx > $scale;
1340 $scale = $ly if $ly > $scale;
1341 my $diff = $ly - $lx;
1342 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
1344 # cases like $x /= $x (but not $x /= $y!) were wrong due to modifying $x
1346 require Scalar::Util;
1347 if (Scalar::Util::refaddr($x) == Scalar::Util::refaddr($y))
1349 $x->bone(); # x/x => 1, rem 0
1354 # make copy of $x in case of list context for later reminder calculation
1355 if (wantarray && !$y->is_one())
1360 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1362 # check for / +-1 ( +/- 1E0)
1365 # promote BigInts and it's subclasses (except when already a BigFloat)
1366 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1368 # calculate the result to $scale digits and then round it
1369 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1370 $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1371 $MBI->_div ($x->{_m},$y->{_m}); # a/c
1373 # correct exponent of $x
1374 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1375 # correct for 10**scale
1376 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1377 $x->bnorm(); # remove trailing 0's
1379 } # ende else $x != $y
1381 # shortcut to not run through _find_round_parameters again
1382 if (defined $params[0])
1384 delete $x->{_a}; # clear before round
1385 $x->bround($params[0],$params[2]); # then round accordingly
1389 delete $x->{_p}; # clear before round
1390 $x->bfround($params[1],$params[2]); # then round accordingly
1394 # clear a/p after round, since user did not request it
1395 delete $x->{_a}; delete $x->{_p};
1402 $rem->bmod($y,@params); # copy already done
1406 # clear a/p after round, since user did not request it
1407 delete $rem->{_a}; delete $rem->{_p};
1416 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1419 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1420 # objectify is costly, so avoid it
1421 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1423 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1426 # handle NaN, inf, -inf
1427 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1429 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1430 $x->{sign} = $re->{sign};
1431 $x->{_e} = $re->{_e};
1432 $x->{_m} = $re->{_m};
1433 return $x->round($a,$p,$r,$y);
1437 return $x->bnan() if $x->is_zero();
1440 return $x->bzero() if $y->is_one() || $x->is_zero();
1442 my $cmp = $x->bacmp($y); # equal or $x < $y?
1443 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1445 # only $y of the operands negative?
1446 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1448 $x->{sign} = $y->{sign}; # calc sign first
1449 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1451 my $ym = $MBI->_copy($y->{_m});
1454 $MBI->_lsft( $ym, $y->{_e}, 10)
1455 if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1457 # if $y has digits after dot
1458 my $shifty = 0; # correct _e of $x by this
1459 if ($y->{_es} eq '-') # has digits after dot
1461 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1462 $shifty = $MBI->_num($y->{_e}); # no more digits after dot
1463 $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1465 # $ym is now mantissa of $y based on exponent 0
1467 my $shiftx = 0; # correct _e of $x by this
1468 if ($x->{_es} eq '-') # has digits after dot
1470 # 123.4 % 20 => 1234 % 200
1471 $shiftx = $MBI->_num($x->{_e}); # no more digits after dot
1472 $MBI->_lsft($ym, $x->{_e}, 10); # 123 => 1230
1474 # 123e1 % 20 => 1230 % 20
1475 if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1477 $MBI->_lsft( $x->{_m}, $x->{_e},10); # es => '+' here
1480 $x->{_e} = $MBI->_new($shiftx);
1482 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1483 $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1485 # now mantissas are equalized, exponent of $x is adjusted, so calc result
1487 $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1489 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1492 if ($neg != 0) # one of them negative => correct in place
1495 $x->{_m} = $r->{_m};
1496 $x->{_e} = $r->{_e};
1497 $x->{_es} = $r->{_es};
1498 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1502 $x->round($a,$p,$r,$y); # round and return
1507 # calculate $y'th root of $x
1510 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1511 # objectify is costly, so avoid it
1512 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1514 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1517 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1518 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1519 $y->{sign} !~ /^\+$/;
1521 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1523 # we need to limit the accuracy to protect against overflow
1525 my (@params,$scale);
1526 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1528 return $x if $x->is_nan(); # error in _find_round_parameters?
1530 # no rounding at all, so must use fallback
1531 if (scalar @params == 0)
1533 # simulate old behaviour
1534 $params[0] = $self->div_scale(); # and round to it as accuracy
1535 $scale = $params[0]+4; # at least four more for proper round
1536 $params[2] = $r; # iound mode by caller or undef
1537 $fallback = 1; # to clear a/p afterwards
1541 # the 4 below is empirical, and there might be cases where it is not
1543 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1546 # when user set globals, they would interfere with our calculation, so
1547 # disable them and later re-enable them
1549 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1550 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1551 # we also need to disable any set A or P on $x (_find_round_parameters took
1552 # them already into account), since these would interfere, too
1553 delete $x->{_a}; delete $x->{_p};
1554 # need to disable $upgrade in BigInt, to avoid deep recursion
1555 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1557 # remember sign and make $x positive, since -4 ** (1/2) => -2
1558 my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
1561 if ($y->isa('Math::BigFloat'))
1563 $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
1567 $is_two = ($y == 2);
1570 # normal square root if $y == 2:
1573 $x->bsqrt($scale+4);
1575 elsif ($y->is_one('-'))
1578 my $u = $self->bone()->bdiv($x,$scale);
1579 # copy private parts over
1580 $x->{_m} = $u->{_m};
1581 $x->{_e} = $u->{_e};
1582 $x->{_es} = $u->{_es};
1586 # calculate the broot() as integer result first, and if it fits, return
1587 # it rightaway (but only if $x and $y are integer):
1589 my $done = 0; # not yet
1590 if ($y->is_int() && $x->is_int())
1592 my $i = $MBI->_copy( $x->{_m} );
1593 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1594 my $int = Math::BigInt->bzero();
1596 $int->broot($y->as_number());
1598 if ($int->copy()->bpow($y) == $x)
1600 # found result, return it
1601 $x->{_m} = $int->{value};
1602 $x->{_e} = $MBI->_zero();
1610 my $u = $self->bone()->bdiv($y,$scale+4);
1611 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
1612 $x->bpow($u,$scale+4); # el cheapo
1615 $x->bneg() if $sign == 1;
1617 # shortcut to not run through _find_round_parameters again
1618 if (defined $params[0])
1620 $x->bround($params[0],$params[2]); # then round accordingly
1624 $x->bfround($params[1],$params[2]); # then round accordingly
1628 # clear a/p after round, since user did not request it
1629 delete $x->{_a}; delete $x->{_p};
1632 $$abr = $ab; $$pbr = $pb;
1638 # calculate square root
1639 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1641 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
1642 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
1643 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
1645 # we need to limit the accuracy to protect against overflow
1647 my (@params,$scale);
1648 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1650 return $x if $x->is_nan(); # error in _find_round_parameters?
1652 # no rounding at all, so must use fallback
1653 if (scalar @params == 0)
1655 # simulate old behaviour
1656 $params[0] = $self->div_scale(); # and round to it as accuracy
1657 $scale = $params[0]+4; # at least four more for proper round
1658 $params[2] = $r; # round mode by caller or undef
1659 $fallback = 1; # to clear a/p afterwards
1663 # the 4 below is empirical, and there might be cases where it is not
1665 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1668 # when user set globals, they would interfere with our calculation, so
1669 # disable them and later re-enable them
1671 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1672 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1673 # we also need to disable any set A or P on $x (_find_round_parameters took
1674 # them already into account), since these would interfere, too
1675 delete $x->{_a}; delete $x->{_p};
1676 # need to disable $upgrade in BigInt, to avoid deep recursion
1677 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1679 my $i = $MBI->_copy( $x->{_m} );
1680 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1681 my $xas = Math::BigInt->bzero();
1684 my $gs = $xas->copy()->bsqrt(); # some guess
1686 if (($x->{_es} ne '-') # guess can't be accurate if there are
1687 # digits after the dot
1688 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
1690 # exact result, copy result over to keep $x
1691 $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
1693 # shortcut to not run through _find_round_parameters again
1694 if (defined $params[0])
1696 $x->bround($params[0],$params[2]); # then round accordingly
1700 $x->bfround($params[1],$params[2]); # then round accordingly
1704 # clear a/p after round, since user did not request it
1705 delete $x->{_a}; delete $x->{_p};
1707 # re-enable A and P, upgrade is taken care of by "local"
1708 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1712 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1713 # of the result by multipyling the input by 100 and then divide the integer
1714 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
1716 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
1717 my $y1 = $MBI->_copy($x->{_m});
1719 my $length = $MBI->_len($y1);
1721 # Now calculate how many digits the result of sqrt(y1) would have
1722 my $digits = int($length / 2);
1724 # But we need at least $scale digits, so calculate how many are missing
1725 my $shift = $scale - $digits;
1727 # That should never happen (we take care of integer guesses above)
1728 # $shift = 0 if $shift < 0;
1730 # Multiply in steps of 100, by shifting left two times the "missing" digits
1731 my $s2 = $shift * 2;
1733 # We now make sure that $y1 has the same odd or even number of digits than
1734 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1735 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1736 # steps of 10. The length of $x does not count, since an even or odd number
1737 # of digits before the dot is not changed by adding an even number of digits
1738 # after the dot (the result is still odd or even digits long).
1739 $s2++ if $MBI->_is_odd($x->{_e});
1741 $MBI->_lsft( $y1, $MBI->_new($s2), 10);
1743 # now take the square root and truncate to integer
1744 $y1 = $MBI->_sqrt($y1);
1746 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1747 # result, which is than later rounded to the desired scale.
1749 # calculate how many zeros $x had after the '.' (or before it, depending
1750 # on sign of $dat, the result should have half as many:
1751 my $dat = $MBI->_num($x->{_e});
1752 $dat = -$dat if $x->{_es} eq '-';
1757 # no zeros after the dot (e.g. 1.23, 0.49 etc)
1758 # preserve half as many digits before the dot than the input had
1759 # (but round this "up")
1760 $dat = int(($dat+1)/2);
1764 $dat = int(($dat)/2);
1766 $dat -= $MBI->_len($y1);
1770 $x->{_e} = $MBI->_new( $dat );
1775 $x->{_e} = $MBI->_new( $dat );
1781 # shortcut to not run through _find_round_parameters again
1782 if (defined $params[0])
1784 $x->bround($params[0],$params[2]); # then round accordingly
1788 $x->bfround($params[1],$params[2]); # then round accordingly
1792 # clear a/p after round, since user did not request it
1793 delete $x->{_a}; delete $x->{_p};
1796 $$abr = $ab; $$pbr = $pb;
1802 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1803 # compute factorial number, modifies first argument
1806 my ($self,$x,@r) = (ref($_[0]),@_);
1807 # objectify is costly, so avoid it
1808 ($self,$x,@r) = objectify(1,@_) if !ref($x);
1810 return $x if $x->{sign} eq '+inf'; # inf => inf
1812 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1813 ($x->{_es} ne '+')); # digits after dot?
1815 # use BigInt's bfac() for faster calc
1816 if (! $MBI->_is_zero($x->{_e}))
1818 $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
1819 $x->{_e} = $MBI->_zero(); # normalize
1822 $MBI->_fac($x->{_m}); # calculate factorial
1823 $x->bnorm()->round(@r); # norm again and round result
1828 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1829 my ($x,$y,$a,$p,$r) = @_;
1832 # if $y == 0.5, it is sqrt($x)
1833 $HALF = $self->new($HALF) unless ref($HALF);
1834 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp($HALF) == 0;
1837 # a ** x == e ** (x * ln a)
1841 # Taylor: | u u^2 u^3 |
1842 # x ** y = 1 + | --- + --- + ----- + ... |
1845 # we need to limit the accuracy to protect against overflow
1847 my ($scale,@params);
1848 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1850 return $x if $x->is_nan(); # error in _find_round_parameters?
1852 # no rounding at all, so must use fallback
1853 if (scalar @params == 0)
1855 # simulate old behaviour
1856 $params[0] = $self->div_scale(); # and round to it as accuracy
1857 $params[1] = undef; # disable P
1858 $scale = $params[0]+4; # at least four more for proper round
1859 $params[2] = $r; # round mode by caller or undef
1860 $fallback = 1; # to clear a/p afterwards
1864 # the 4 below is empirical, and there might be cases where it is not
1866 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1869 # when user set globals, they would interfere with our calculation, so
1870 # disable them and later re-enable them
1872 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1873 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1874 # we also need to disable any set A or P on $x (_find_round_parameters took
1875 # them already into account), since these would interfere, too
1876 delete $x->{_a}; delete $x->{_p};
1877 # need to disable $upgrade in BigInt, to avoid deep recursion
1878 local $Math::BigInt::upgrade = undef;
1880 my ($limit,$v,$u,$below,$factor,$next,$over);
1882 $u = $x->copy()->blog(undef,$scale)->bmul($y);
1883 $v = $self->bone(); # 1
1884 $factor = $self->new(2); # 2
1885 $x->bone(); # first term: 1
1887 $below = $v->copy();
1890 $limit = $self->new("1E-". ($scale-1));
1894 # we calculate the next term, and add it to the last
1895 # when the next term is below our limit, it won't affect the outcome
1896 # anymore, so we stop
1897 $next = $over->copy()->bdiv($below,$scale);
1898 last if $next->bacmp($limit) <= 0;
1900 # calculate things for the next term
1901 $over *= $u; $below *= $factor; $factor->binc();
1903 last if $x->{sign} !~ /^[-+]$/;
1908 # shortcut to not run through _find_round_parameters again
1909 if (defined $params[0])
1911 $x->bround($params[0],$params[2]); # then round accordingly
1915 $x->bfround($params[1],$params[2]); # then round accordingly
1919 # clear a/p after round, since user did not request it
1920 delete $x->{_a}; delete $x->{_p};
1923 $$abr = $ab; $$pbr = $pb;
1929 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1930 # compute power of two numbers, second arg is used as integer
1931 # modifies first argument
1934 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1935 # objectify is costly, so avoid it
1936 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1938 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1941 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1942 return $x if $x->{sign} =~ /^[+-]inf$/;
1945 return $x->bnan() if $x->{sign} eq '-' && $y->{sign} eq '-';
1947 # cache the result of is_zero
1948 my $y_is_zero = $y->is_zero();
1949 return $x->bone() if $y_is_zero;
1950 return $x if $x->is_one() || $y->is_one();
1952 my $x_is_zero = $x->is_zero();
1953 return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int(); # non-integer power
1955 my $y1 = $y->as_number()->{value}; # make MBI part
1958 if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
1960 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
1961 return $MBI->_is_odd($y1) ? $x : $x->babs(1);
1965 return $x->bone() if $y_is_zero;
1966 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1967 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
1972 $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
1974 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1975 $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
1976 $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
1978 $x->{sign} = $new_sign;
1980 if ($y->{sign} eq '-')
1982 # modify $x in place!
1983 my $z = $x->copy(); $x->bone();
1984 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
1986 $x->round($a,$p,$r,$y);
1989 ###############################################################################
1990 # rounding functions
1994 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1995 # $n == 0 means round to integer
1996 # expects and returns normalized numbers!
1997 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1999 my ($scale,$mode) = $x->_scale_p(@_);
2000 return $x if !defined $scale || $x->modify('bfround'); # no-op
2002 # never round a 0, +-inf, NaN
2005 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
2008 return $x if $x->{sign} !~ /^[+-]$/;
2010 # don't round if x already has lower precision
2011 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
2013 $x->{_p} = $scale; # remember round in any case
2014 delete $x->{_a}; # and clear A
2017 # round right from the '.'
2019 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round
2021 $scale = -$scale; # positive for simplicity
2022 my $len = $MBI->_len($x->{_m}); # length of mantissa
2024 # the following poses a restriction on _e, but if _e is bigger than a
2025 # scalar, you got other problems (memory etc) anyway
2026 my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e}))); # digits after dot
2027 my $zad = 0; # zeros after dot
2028 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
2030 # p rint "scale $scale dad $dad zad $zad len $len\n";
2031 # number bsstr len zad dad
2032 # 0.123 123e-3 3 0 3
2033 # 0.0123 123e-4 3 1 4
2036 # 1.2345 12345e-4 5 0 4
2038 # do not round after/right of the $dad
2039 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
2041 # round to zero if rounding inside the $zad, but not for last zero like:
2042 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
2043 return $x->bzero() if $scale < $zad;
2044 if ($scale == $zad) # for 0.006, scale -3 and trunc
2050 # adjust round-point to be inside mantissa
2053 $scale = $scale-$zad;
2057 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
2058 $scale = $dbd+$scale;
2064 # round left from the '.'
2066 # 123 => 100 means length(123) = 3 - $scale (2) => 1
2068 my $dbt = $MBI->_len($x->{_m});
2070 my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
2071 # should be the same, so treat it as this
2072 $scale = 1 if $scale == 0;
2073 # shortcut if already integer
2074 return $x if $scale == 1 && $dbt <= $dbd;
2075 # maximum digits before dot
2080 # not enough digits before dot, so round to zero
2083 elsif ( $scale == $dbd )
2090 $scale = $dbd - $scale;
2093 # pass sign to bround for rounding modes '+inf' and '-inf'
2094 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2095 $m->bround($scale,$mode);
2096 $x->{_m} = $m->{value}; # get our mantissa back
2102 # accuracy: preserve $N digits, and overwrite the rest with 0's
2103 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2105 if (($_[0] || 0) < 0)
2107 require Carp; Carp::croak ('bround() needs positive accuracy');
2110 my ($scale,$mode) = $x->_scale_a(@_);
2111 return $x if !defined $scale || $x->modify('bround'); # no-op
2113 # scale is now either $x->{_a}, $accuracy, or the user parameter
2114 # test whether $x already has lower accuracy, do nothing in this case
2115 # but do round if the accuracy is the same, since a math operation might
2116 # want to round a number with A=5 to 5 digits afterwards again
2117 return $x if defined $x->{_a} && $x->{_a} < $scale;
2119 # scale < 0 makes no sense
2120 # scale == 0 => keep all digits
2121 # never round a +-inf, NaN
2122 return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
2124 # 1: never round a 0
2125 # 2: if we should keep more digits than the mantissa has, do nothing
2126 if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
2128 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
2132 # pass sign to bround for '+inf' and '-inf' rounding modes
2133 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2135 $m->bround($scale,$mode); # round mantissa
2136 $x->{_m} = $m->{value}; # get our mantissa back
2137 $x->{_a} = $scale; # remember rounding
2138 delete $x->{_p}; # and clear P
2139 $x->bnorm(); # del trailing zeros gen. by bround()
2144 # return integer less or equal then $x
2145 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2147 return $x if $x->modify('bfloor');
2149 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2151 # if $x has digits after dot
2152 if ($x->{_es} eq '-')
2154 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2155 $x->{_e} = $MBI->_zero(); # trunc/norm
2156 $x->{_es} = '+'; # abs e
2157 $MBI->_inc($x->{_m}) if $x->{sign} eq '-'; # increment if negative
2159 $x->round($a,$p,$r);
2164 # return integer greater or equal then $x
2165 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2167 return $x if $x->modify('bceil');
2168 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2170 # if $x has digits after dot
2171 if ($x->{_es} eq '-')
2173 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2174 $x->{_e} = $MBI->_zero(); # trunc/norm
2175 $x->{_es} = '+'; # abs e
2176 $MBI->_inc($x->{_m}) if $x->{sign} eq '+'; # increment if positive
2178 $x->round($a,$p,$r);
2183 # shift right by $y (divide by power of $n)
2186 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2187 # objectify is costly, so avoid it
2188 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2190 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2193 return $x if $x->modify('brsft');
2194 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2196 $n = 2 if !defined $n; $n = $self->new($n);
2197 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
2202 # shift left by $y (multiply by power of $n)
2205 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2206 # objectify is costly, so avoid it
2207 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2209 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2212 return $x if $x->modify('blsft');
2213 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2215 $n = 2 if !defined $n; $n = $self->new($n);
2216 $x->bmul($n->bpow($y),$a,$p,$r,$y);
2219 ###############################################################################
2223 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
2228 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
2229 # or falling back to MBI::bxxx()
2230 my $name = $AUTOLOAD;
2232 $name =~ s/(.*):://; # split package
2233 my $c = $1 || $class;
2235 $c->import() if $IMPORT == 0;
2236 if (!method_alias($name))
2240 # delayed load of Carp and avoid recursion
2242 Carp::croak ("$c: Can't call a method without name");
2244 if (!method_hand_up($name))
2246 # delayed load of Carp and avoid recursion
2248 Carp::croak ("Can't call $c\-\>$name, not a valid method");
2250 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
2252 return &{"Math::BigInt"."::$name"}(@_);
2254 my $bname = $name; $bname =~ s/^f/b/;
2262 # return a copy of the exponent
2263 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2265 if ($x->{sign} !~ /^[+-]$/)
2267 my $s = $x->{sign}; $s =~ s/^[+-]//;
2268 return Math::BigInt->new($s); # -inf, +inf => +inf
2270 Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
2275 # return a copy of the mantissa
2276 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2278 if ($x->{sign} !~ /^[+-]$/)
2280 my $s = $x->{sign}; $s =~ s/^[+]//;
2281 return Math::BigInt->new($s); # -inf, +inf => +inf
2283 my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
2284 $m->bneg() if $x->{sign} eq '-';
2291 # return a copy of both the exponent and the mantissa
2292 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2294 if ($x->{sign} !~ /^[+-]$/)
2296 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2297 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2299 my $m = Math::BigInt->bzero();
2300 $m->{value} = $MBI->_copy($x->{_m});
2301 $m->bneg() if $x->{sign} eq '-';
2302 ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
2305 ##############################################################################
2306 # private stuff (internal use only)
2312 my $lib = ''; my @a;
2314 for ( my $i = 0; $i < $l ; $i++)
2316 if ( $_[$i] eq ':constant' )
2318 # This causes overlord er load to step in. 'binary' and 'integer'
2319 # are handled by BigInt.
2320 overload::constant float => sub { $self->new(shift); };
2322 elsif ($_[$i] eq 'upgrade')
2324 # this causes upgrading
2325 $upgrade = $_[$i+1]; # or undef to disable
2328 elsif ($_[$i] eq 'downgrade')
2330 # this causes downgrading
2331 $downgrade = $_[$i+1]; # or undef to disable
2334 elsif ($_[$i] eq 'lib')
2336 # alternative library
2337 $lib = $_[$i+1] || ''; # default Calc
2340 elsif ($_[$i] eq 'with')
2342 # alternative class for our private parts()
2343 # XXX: no longer supported
2344 # $MBI = $_[$i+1] || 'Math::BigInt';
2353 $lib =~ tr/a-zA-Z0-9,://cd; # restrict to sane characters
2354 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2355 my $mbilib = eval { Math::BigInt->config()->{lib} };
2356 if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
2358 # MBI already loaded
2359 Math::BigInt->import('lib',"$lib,$mbilib", 'objectify');
2363 # MBI not loaded, or with ne "Math::BigInt::Calc"
2364 $lib .= ",$mbilib" if defined $mbilib;
2365 $lib =~ s/^,//; # don't leave empty
2367 # replacement library can handle lib statement, but also could ignore it
2369 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2370 # used in the same script, or eval inside import(). So we require MBI:
2371 require Math::BigInt;
2372 Math::BigInt->import( lib => $lib, 'objectify' );
2376 require Carp; Carp::croak ("Couldn't load $lib: $! $@");
2378 # find out which one was actually loaded
2379 $MBI = Math::BigInt->config()->{lib};
2381 # register us with MBI to get notified of future lib changes
2382 Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
2384 # any non :constant stuff is handled by our parent, Exporter
2385 # even if @_ is empty, to give it a chance
2386 $self->SUPER::import(@a); # for subclasses
2387 $self->export_to_level(1,$self,@a); # need this, too
2392 # adjust m and e so that m is smallest possible
2393 # round number according to accuracy and precision settings
2394 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
2396 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2398 my $zeros = $MBI->_zeros($x->{_m}); # correct for trailing zeros
2401 my $z = $MBI->_new($zeros);
2402 $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
2403 if ($x->{_es} eq '-')
2405 if ($MBI->_acmp($x->{_e},$z) >= 0)
2407 $x->{_e} = $MBI->_sub ($x->{_e}, $z);
2408 $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
2412 $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
2418 $x->{_e} = $MBI->_add ($x->{_e}, $z);
2423 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
2424 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
2425 $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
2426 if $MBI->_is_zero($x->{_m});
2429 $x; # MBI bnorm is no-op, so dont call it
2432 ##############################################################################
2436 # return number as hexadecimal string (only for integers defined)
2437 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2439 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2440 return '0x0' if $x->is_zero();
2442 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2444 my $z = $MBI->_copy($x->{_m});
2445 if (! $MBI->_is_zero($x->{_e})) # > 0
2447 $MBI->_lsft( $z, $x->{_e},10);
2449 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2455 # return number as binary digit string (only for integers defined)
2456 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2458 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2459 return '0b0' if $x->is_zero();
2461 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2463 my $z = $MBI->_copy($x->{_m});
2464 if (! $MBI->_is_zero($x->{_e})) # > 0
2466 $MBI->_lsft( $z, $x->{_e},10);
2468 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2474 # return copy as a bigint representation of this BigFloat number
2475 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2477 my $z = $MBI->_copy($x->{_m});
2478 if ($x->{_es} eq '-') # < 0
2480 $MBI->_rsft( $z, $x->{_e},10);
2482 elsif (! $MBI->_is_zero($x->{_e})) # > 0
2484 $MBI->_lsft( $z, $x->{_e},10);
2486 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2493 my $class = ref($x) || $x;
2494 $x = $class->new(shift) unless ref($x);
2496 return 1 if $MBI->_is_zero($x->{_m});
2498 my $len = $MBI->_len($x->{_m});
2499 $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
2503 $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
2514 Math::BigFloat - Arbitrary size floating point math package
2521 $x = Math::BigFloat->new($str); # defaults to 0
2522 $nan = Math::BigFloat->bnan(); # create a NotANumber
2523 $zero = Math::BigFloat->bzero(); # create a +0
2524 $inf = Math::BigFloat->binf(); # create a +inf
2525 $inf = Math::BigFloat->binf('-'); # create a -inf
2526 $one = Math::BigFloat->bone(); # create a +1
2527 $one = Math::BigFloat->bone('-'); # create a -1
2530 $x->is_zero(); # true if arg is +0
2531 $x->is_nan(); # true if arg is NaN
2532 $x->is_one(); # true if arg is +1
2533 $x->is_one('-'); # true if arg is -1
2534 $x->is_odd(); # true if odd, false for even
2535 $x->is_even(); # true if even, false for odd
2536 $x->is_pos(); # true if >= 0
2537 $x->is_neg(); # true if < 0
2538 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2540 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2541 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2542 $x->sign(); # return the sign, either +,- or NaN
2543 $x->digit($n); # return the nth digit, counting from right
2544 $x->digit(-$n); # return the nth digit, counting from left
2546 # The following all modify their first argument. If you want to preserve
2547 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2548 # neccessary when mixing $a = $b assigments with non-overloaded math.
2551 $x->bzero(); # set $i to 0
2552 $x->bnan(); # set $i to NaN
2553 $x->bone(); # set $x to +1
2554 $x->bone('-'); # set $x to -1
2555 $x->binf(); # set $x to inf
2556 $x->binf('-'); # set $x to -inf
2558 $x->bneg(); # negation
2559 $x->babs(); # absolute value
2560 $x->bnorm(); # normalize (no-op)
2561 $x->bnot(); # two's complement (bit wise not)
2562 $x->binc(); # increment x by 1
2563 $x->bdec(); # decrement x by 1
2565 $x->badd($y); # addition (add $y to $x)
2566 $x->bsub($y); # subtraction (subtract $y from $x)
2567 $x->bmul($y); # multiplication (multiply $x by $y)
2568 $x->bdiv($y); # divide, set $x to quotient
2569 # return (quo,rem) or quo if scalar
2571 $x->bmod($y); # modulus ($x % $y)
2572 $x->bpow($y); # power of arguments ($x ** $y)
2573 $x->blsft($y); # left shift
2574 $x->brsft($y); # right shift
2575 # return (quo,rem) or quo if scalar
2577 $x->blog(); # logarithm of $x to base e (Euler's number)
2578 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
2580 $x->band($y); # bit-wise and
2581 $x->bior($y); # bit-wise inclusive or
2582 $x->bxor($y); # bit-wise exclusive or
2583 $x->bnot(); # bit-wise not (two's complement)
2585 $x->bsqrt(); # calculate square-root
2586 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
2587 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2589 $x->bround($N); # accuracy: preserve $N digits
2590 $x->bfround($N); # precision: round to the $Nth digit
2592 $x->bfloor(); # return integer less or equal than $x
2593 $x->bceil(); # return integer greater or equal than $x
2595 # The following do not modify their arguments:
2597 bgcd(@values); # greatest common divisor
2598 blcm(@values); # lowest common multiplicator
2600 $x->bstr(); # return string
2601 $x->bsstr(); # return string in scientific notation
2603 $x->as_int(); # return $x as BigInt
2604 $x->exponent(); # return exponent as BigInt
2605 $x->mantissa(); # return mantissa as BigInt
2606 $x->parts(); # return (mantissa,exponent) as BigInt
2608 $x->length(); # number of digits (w/o sign and '.')
2609 ($l,$f) = $x->length(); # number of digits, and length of fraction
2611 $x->precision(); # return P of $x (or global, if P of $x undef)
2612 $x->precision($n); # set P of $x to $n
2613 $x->accuracy(); # return A of $x (or global, if A of $x undef)
2614 $x->accuracy($n); # set A $x to $n
2616 # these get/set the appropriate global value for all BigFloat objects
2617 Math::BigFloat->precision(); # Precision
2618 Math::BigFloat->accuracy(); # Accuracy
2619 Math::BigFloat->round_mode(); # rounding mode
2623 All operators (inlcuding basic math operations) are overloaded if you
2624 declare your big floating point numbers as
2626 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2628 Operations with overloaded operators preserve the arguments, which is
2629 exactly what you expect.
2631 =head2 Canonical notation
2633 Input to these routines are either BigFloat objects, or strings of the
2634 following four forms:
2648 C</^[+-]\d+E[+-]?\d+$/>
2652 C</^[+-]\d*\.\d+E[+-]?\d+$/>
2656 all with optional leading and trailing zeros and/or spaces. Additonally,
2657 numbers are allowed to have an underscore between any two digits.
2659 Empty strings as well as other illegal numbers results in 'NaN'.
2661 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2662 are always stored in normalized form. On a string, it creates a BigFloat
2667 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2669 The string output will always have leading and trailing zeros stripped and drop
2670 a plus sign. C<bstr()> will give you always the form with a decimal point,
2671 while C<bsstr()> (s for scientific) gives you the scientific notation.
2673 Input bstr() bsstr()
2675 ' -123 123 123' '-123123123' '-123123123E0'
2676 '00.0123' '0.0123' '123E-4'
2677 '123.45E-2' '1.2345' '12345E-4'
2678 '10E+3' '10000' '1E4'
2680 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2681 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2682 return either undef, <0, 0 or >0 and are suited for sort.
2684 Actual math is done by using the class defined with C<with => Class;> (which
2685 defaults to BigInts) to represent the mantissa and exponent.
2687 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2688 represent the result when input arguments are not numbers, as well as
2689 the result of dividing by zero.
2691 =head2 C<mantissa()>, C<exponent()> and C<parts()>
2693 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2694 as BigInts such that:
2696 $m = $x->mantissa();
2697 $e = $x->exponent();
2698 $y = $m * ( 10 ** $e );
2699 print "ok\n" if $x == $y;
2701 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2703 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2705 Currently the mantissa is reduced as much as possible, favouring higher
2706 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2707 This might change in the future, so do not depend on it.
2709 =head2 Accuracy vs. Precision
2711 See also: L<Rounding|Rounding>.
2713 Math::BigFloat supports both precision and accuracy. For a full documentation,
2714 examples and tips on these topics please see the large section in
2717 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2718 a operation consumes all resources, each operation produces no more than
2719 the requested number of digits.
2721 Please refer to BigInt's documentation for the precedence rules of which
2722 accuracy/precision setting will be used.
2724 If there is no gloabl precision set, B<and> the operation inquestion was not
2725 called with a requested precision or accuracy, B<and> the input $x has no
2726 accuracy or precision set, then a fallback parameter will be used. For
2727 historical reasons, it is called C<div_scale> and can be accessed via:
2729 $d = Math::BigFloat->div_scale(); # query
2730 Math::BigFloat->div_scale($n); # set to $n digits
2732 The default value is 40 digits.
2734 In case the result of one operation has more precision than specified,
2735 it is rounded. The rounding mode taken is either the default mode, or the one
2736 supplied to the operation after the I<scale>:
2738 $x = Math::BigFloat->new(2);
2739 Math::BigFloat->precision(5); # 5 digits max
2740 $y = $x->copy()->bdiv(3); # will give 0.66666
2741 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2742 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
2743 Math::BigFloat->round_mode('zero');
2744 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2750 =item ffround ( +$scale )
2752 Rounds to the $scale'th place left from the '.', counting from the dot.
2753 The first digit is numbered 1.
2755 =item ffround ( -$scale )
2757 Rounds to the $scale'th place right from the '.', counting from the dot.
2761 Rounds to an integer.
2763 =item fround ( +$scale )
2765 Preserves accuracy to $scale digits from the left (aka significant digits)
2766 and pads the rest with zeros. If the number is between 1 and -1, the
2767 significant digits count from the first non-zero after the '.'
2769 =item fround ( -$scale ) and fround ( 0 )
2771 These are effectively no-ops.
2775 All rounding functions take as a second parameter a rounding mode from one of
2776 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2778 The default rounding mode is 'even'. By using
2779 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
2780 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2781 no longer supported.
2782 The second parameter to the round functions then overrides the default
2785 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
2786 'trunc' as rounding mode to make it equivalent to:
2791 You can override this by passing the desired rounding mode as parameter to
2794 $x = Math::BigFloat->new(2.5);
2795 $y = $x->as_number('odd'); # $y = 3
2801 =head1 Autocreating constants
2803 After C<use Math::BigFloat ':constant'> all the floating point constants
2804 in the given scope are converted to C<Math::BigFloat>. This conversion
2805 happens at compile time.
2809 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2811 prints the value of C<2E-100>. Note that without conversion of
2812 constants the expression 2E-100 will be calculated as normal floating point
2815 Please note that ':constant' does not affect integer constants, nor binary
2816 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2821 Math with the numbers is done (by default) by a module called
2822 Math::BigInt::Calc. This is equivalent to saying:
2824 use Math::BigFloat lib => 'Calc';
2826 You can change this by using:
2828 use Math::BigFloat lib => 'BitVect';
2830 The following would first try to find Math::BigInt::Foo, then
2831 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2833 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2835 Calc.pm uses as internal format an array of elements of some decimal base
2836 (usually 1e7, but this might be differen for some systems) with the least
2837 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2838 significant bit first. Other modules might use even different means of
2839 representing the numbers. See the respective module documentation for further
2842 Please note that Math::BigFloat does B<not> use the denoted library itself,
2843 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2846 use Math::BigInt lib => 'GMP';
2849 you can roll it all into one line:
2851 use Math::BigFloat lib => 'GMP';
2853 It is also possible to just require Math::BigFloat:
2855 require Math::BigFloat;
2857 This will load the neccessary things (like BigInt) when they are needed, and
2860 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2861 you ever wanted to know about loading a different library.
2863 =head2 Using Math::BigInt::Lite
2865 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2868 use Math::BigFloat with => 'Math::BigInt::Lite';
2870 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2871 can combine these if you want. For instance, you may want to use
2872 Math::BigInt objects in your main script, too.
2876 use Math::BigFloat with => 'Math::BigInt::Lite';
2878 Of course, you can combine this with the C<lib> parameter.
2881 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2883 There is no need for a "use Math::BigInt;" statement, even if you want to
2884 use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2885 always loads it. But if you add it, add it B<before>:
2889 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2891 Notice that the module with the last C<lib> will "win" and thus
2892 it's lib will be used if the lib is available:
2895 use Math::BigInt lib => 'Bar,Baz';
2896 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2898 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2899 words, Math::BigFloat will try to retain previously loaded libs when you
2900 don't specify it onem but if you specify one, it will try to load them.
2902 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2903 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2904 same as trying the latter load alone, except for the fact that one of Bar or
2905 Baz might be loaded needlessly in an intermidiate step (and thus hang around
2906 and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2907 will still be tried to be loaded, but this is not as time/memory consuming as
2908 actually loading one of them. Still, this type of usage is not recommended due
2911 The old way (loading the lib only in BigInt) still works though:
2914 use Math::BigInt lib => 'Bar,Baz';
2917 You can even load Math::BigInt afterwards:
2921 use Math::BigInt lib => 'Bar,Baz';
2923 But this has the same problems like #5, it will first load Calc
2924 (Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2925 Baz, depending on which of them works and is usable/loadable. Since this
2926 loads Calc unnecc., it is not recommended.
2928 Since it also possible to just require Math::BigFloat, this poses the question
2929 about what libary this will use:
2931 require Math::BigFloat;
2932 my $x = Math::BigFloat->new(123); $x += 123;
2934 It will use Calc. Please note that the call to import() is still done, but
2935 only when you use for the first time some Math::BigFloat math (it is triggered
2936 via any constructor, so the first time you create a Math::BigFloat, the load
2937 will happen in the background). This means:
2939 require Math::BigFloat;
2940 Math::BigFloat->import ( lib => 'Foo,Bar' );
2942 would be the same as:
2944 use Math::BigFloat lib => 'Foo, Bar';
2946 But don't try to be clever to insert some operations in between:
2948 require Math::BigFloat;
2949 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
2950 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
2951 $x = Math::BigFloat->bone()+4; # now use Pari
2953 While this works, it loads Calc needlessly. But maybe you just wanted that?
2955 B<Examples #3 is highly recommended> for daily usage.
2959 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
2965 =item stringify, bstr()
2967 Both stringify and bstr() now drop the leading '+'. The old code would return
2968 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2969 reasoning and details.
2973 The following will probably not do what you expect:
2975 print $c->bdiv(123.456),"\n";
2977 It prints both quotient and reminder since print works in list context. Also,
2978 bdiv() will modify $c, so be carefull. You probably want to use
2980 print $c / 123.456,"\n";
2981 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2985 =item Modifying and =
2989 $x = Math::BigFloat->new(5);
2992 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2993 a second reference to the B<same> object and stores it in $y. Thus anything
2994 that modifies $x will modify $y (except overloaded math operators), and vice
2995 versa. See L<Math::BigInt> for details and how to avoid that.
2999 C<bpow()> now modifies the first argument, unlike the old code which left
3000 it alone and only returned the result. This is to be consistent with
3001 C<badd()> etc. The first will modify $x, the second one won't:
3003 print bpow($x,$i),"\n"; # modify $x
3004 print $x->bpow($i),"\n"; # ditto
3005 print $x ** $i,"\n"; # leave $x alone
3011 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
3012 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
3014 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
3015 because they solve the autoupgrading/downgrading issue, at least partly.
3018 L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
3019 more documentation including a full version history, testcases, empty
3020 subclass files and benchmarks.
3024 This program is free software; you may redistribute it and/or modify it under
3025 the same terms as Perl itself.
3029 Mark Biggar, overloaded interface by Ilya Zakharevich.
3030 Completely rewritten by Tels L<http://bloodgate.com> in 2001 - 2004, and still