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
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);
629 # sub bsub is inherited from Math::BigInt!
633 # increment arg by one
634 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
636 if ($x->{_es} eq '-')
638 return $x->badd($self->bone(),@r); # digits after dot
641 if (!$MBI->_is_zero($x->{_e})) # _e == 0 for NaN, inf, -inf
643 # 1e2 => 100, so after the shift below _m has a '0' as last digit
644 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
645 $x->{_e} = $MBI->_zero(); # normalize
647 # we know that the last digit of $x will be '1' or '9', depending on the
651 if ($x->{sign} eq '+')
653 $MBI->_inc($x->{_m});
654 return $x->bnorm()->bround(@r);
656 elsif ($x->{sign} eq '-')
658 $MBI->_dec($x->{_m});
659 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
660 return $x->bnorm()->bround(@r);
662 # inf, nan handling etc
663 $x->badd($self->bone(),@r); # badd() does round
668 # decrement arg by one
669 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
671 if ($x->{_es} eq '-')
673 return $x->badd($self->bone('-'),@r); # digits after dot
676 if (!$MBI->_is_zero($x->{_e}))
678 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
679 $x->{_e} = $MBI->_zero(); # normalize
683 my $zero = $x->is_zero();
685 if (($x->{sign} eq '-') || $zero)
687 $MBI->_inc($x->{_m});
688 $x->{sign} = '-' if $zero; # 0 => 1 => -1
689 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
690 return $x->bnorm()->round(@r);
693 elsif ($x->{sign} eq '+')
695 $MBI->_dec($x->{_m});
696 return $x->bnorm()->round(@r);
698 # inf, nan handling etc
699 $x->badd($self->bone('-'),@r); # does round
706 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
708 # $base > 0, $base != 1; if $base == undef default to $base == e
711 # we need to limit the accuracy to protect against overflow
714 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
716 # also takes care of the "error in _find_round_parameters?" case
717 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
720 # no rounding at all, so must use fallback
721 if (scalar @params == 0)
723 # simulate old behaviour
724 $params[0] = $self->div_scale(); # and round to it as accuracy
725 $params[1] = undef; # P = undef
726 $scale = $params[0]+4; # at least four more for proper round
727 $params[2] = $r; # round mode by caller or undef
728 $fallback = 1; # to clear a/p afterwards
732 # the 4 below is empirical, and there might be cases where it is not
734 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
737 return $x->bzero(@params) if $x->is_one();
738 # base not defined => base == Euler's constant e
741 # make object, since we don't feed it through objectify() to still get the
742 # case of $base == undef
743 $base = $self->new($base) unless ref($base);
744 # $base > 0; $base != 1
745 return $x->bnan() if $base->is_zero() || $base->is_one() ||
746 $base->{sign} ne '+';
747 # if $x == $base, we know the result must be 1.0
748 return $x->bone('+',@params) if $x->bcmp($base) == 0;
751 # when user set globals, they would interfere with our calculation, so
752 # disable them and later re-enable them
754 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
755 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
756 # we also need to disable any set A or P on $x (_find_round_parameters took
757 # them already into account), since these would interfere, too
758 delete $x->{_a}; delete $x->{_p};
759 # need to disable $upgrade in BigInt, to avoid deep recursion
760 local $Math::BigInt::upgrade = undef;
761 local $Math::BigFloat::downgrade = undef;
763 # upgrade $x if $x is not a BigFloat (handle BigInt input)
764 if (!$x->isa('Math::BigFloat'))
766 $x = Math::BigFloat->new($x);
772 # If the base is defined and an integer, try to calculate integer result
773 # first. This is very fast, and in case the real result was found, we can
775 if (defined $base && $base->is_int() && $x->is_int())
777 my $i = $MBI->_copy( $x->{_m} );
778 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
779 my $int = Math::BigInt->bzero();
781 $int->blog($base->as_number());
783 if ($base->as_number()->bpow($int) == $x)
785 # found result, return it
786 $x->{_m} = $int->{value};
787 $x->{_e} = $MBI->_zero();
796 # first calculate the log to base e (using reduction by 10 (and probably 2))
797 $self->_log_10($x,$scale);
799 # and if a different base was requested, convert it
802 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
803 # not ln, but some other base (don't modify $base)
804 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
808 # shortcut to not run through _find_round_parameters again
809 if (defined $params[0])
811 $x->bround($params[0],$params[2]); # then round accordingly
815 $x->bfround($params[1],$params[2]); # then round accordingly
819 # clear a/p after round, since user did not request it
820 delete $x->{_a}; delete $x->{_p};
823 $$abr = $ab; $$pbr = $pb;
830 # internal log function to calculate ln() based on Taylor series.
831 # Modifies $x in place.
832 my ($self,$x,$scale) = @_;
834 # in case of $x == 1, result is 0
835 return $x->bzero() if $x->is_one();
837 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
841 # Taylor: | u 1 u^3 1 u^5 |
842 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
843 # |_ v 3 v^3 5 v^5 _|
845 # This takes much more steps to calculate the result and is thus not used
848 # Taylor: | u 1 u^2 1 u^3 |
849 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
850 # |_ x 2 x^2 3 x^3 _|
852 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
854 $v = $x->copy(); $v->binc(); # v = x+1
855 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
856 $x->bdiv($v,$scale); # first term: u/v
859 $u *= $u; $v *= $v; # u^2, v^2
860 $below->bmul($v); # u^3, v^3
862 $factor = $self->new(3); $f = $self->new(2);
864 my $steps = 0 if DEBUG;
865 $limit = $self->new("1E-". ($scale-1));
868 # we calculate the next term, and add it to the last
869 # when the next term is below our limit, it won't affect the outcome
870 # anymore, so we stop
872 # calculating the next term simple from over/below will result in quite
873 # a time hog if the input has many digits, since over and below will
874 # accumulate more and more digits, and the result will also have many
875 # digits, but in the end it is rounded to $scale digits anyway. So if we
876 # round $over and $below first, we save a lot of time for the division
877 # (not with log(1.2345), but try log (123**123) to see what I mean. This
878 # can introduce a rounding error if the division result would be f.i.
879 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
880 # if we truncated $over and $below we might get 0.12345. Does this matter
881 # for the end result? So we give $over and $below 4 more digits to be
882 # on the safe side (unscientific error handling as usual... :+D
884 $next = $over->copy->bround($scale+4)->bdiv(
885 $below->copy->bmul($factor)->bround($scale+4),
889 ## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
891 last if $next->bacmp($limit) <= 0;
893 delete $next->{_a}; delete $next->{_p};
895 # calculate things for the next term
896 $over *= $u; $below *= $v; $factor->badd($f);
899 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
902 $x->bmul($f); # $x *= 2
903 print "took $steps steps\n" if DEBUG;
908 # Internal log function based on reducing input to the range of 0.1 .. 9.99
909 # and then "correcting" the result to the proper one. Modifies $x in place.
910 my ($self,$x,$scale) = @_;
912 # taking blog() from numbers greater than 10 takes a *very long* time, so we
913 # break the computation down into parts based on the observation that:
914 # blog(x*y) = blog(x) + blog(y)
915 # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
916 # the faster it get's, especially because 2*$x takes about 10 times as long,
917 # so by dividing $x by 10 we make it at least factor 100 faster...)
919 # The same observation is valid for numbers smaller than 0.1 (e.g. computing
920 # log(1) is fastest, and the farther away we get from 1, the longer it takes)
921 # so we also 'break' this down by multiplying $x with 10 and subtract the
922 # log(10) afterwards to get the correct result.
924 # calculate nr of digits before dot
925 my $dbd = $MBI->_num($x->{_e});
926 $dbd = -$dbd if $x->{_es} eq '-';
927 $dbd += $MBI->_len($x->{_m});
929 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
932 my $calc = 1; # do some calculation?
934 # disable the shortcut for 10, since we need log(10) and this would recurse
936 if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
938 $dbd = 0; # disable shortcut
939 # we can use the cached value in these cases
940 if ($scale <= $LOG_10_A)
942 $x->bzero(); $x->badd($LOG_10);
943 $calc = 0; # no need to calc, but round
948 # disable the shortcut for 2, since we maybe have it cached
949 if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
951 $dbd = 0; # disable shortcut
952 # we can use the cached value in these cases
953 if ($scale <= $LOG_2_A)
955 $x->bzero(); $x->badd($LOG_2);
956 $calc = 0; # no need to calc, but round
961 # if $x = 0.1, we know the result must be 0-log(10)
962 if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
963 $MBI->_is_one($x->{_m}))
965 $dbd = 0; # disable shortcut
966 # we can use the cached value in these cases
967 if ($scale <= $LOG_10_A)
969 $x->bzero(); $x->bsub($LOG_10);
970 $calc = 0; # no need to calc, but round
974 return if $calc == 0; # already have the result
976 # default: these correction factors are undef and thus not used
977 my $l_10; # value of ln(10) to A of $scale
978 my $l_2; # value of ln(2) to A of $scale
980 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
981 # so don't do this shortcut for 1 or 0
982 if (($dbd > 1) || ($dbd < 0))
984 # convert our cached value to an object if not already (avoid doing this
985 # at import() time, since not everybody needs this)
986 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
988 #print "x = $x, dbd = $dbd, calc = $calc\n";
989 # got more than one digit before the dot, or more than one zero after the
991 # log(123) == log(1.23) + log(10) * 2
992 # log(0.0123) == log(1.23) - log(10) * 2
994 if ($scale <= $LOG_10_A)
997 $l_10 = $LOG_10->copy(); # copy for mul
1001 # else: slower, compute it (but don't cache it, because it could be big)
1002 # also disable downgrade for this code path
1003 local $Math::BigFloat::downgrade = undef;
1004 $l_10 = $self->new(10)->blog(undef,$scale); # scale+4, actually
1006 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
1007 $l_10->bmul( $self->new($dbd)); # log(10) * (digits_before_dot-1)
1014 ($x->{_e}, $x->{_es}) =
1015 _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
1019 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1021 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1022 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1024 $HALF = $self->new($HALF) unless ref($HALF);
1026 my $twos = 0; # default: none (0 times)
1027 my $two = $self->new(2);
1028 while ($x->bacmp($HALF) <= 0)
1030 $twos--; $x->bmul($two);
1032 while ($x->bacmp($two) >= 0)
1034 $twos++; $x->bdiv($two,$scale+4); # keep all digits
1036 # $twos > 0 => did mul 2, < 0 => did div 2 (never both)
1037 # calculate correction factor based on ln(2)
1040 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1041 if ($scale <= $LOG_2_A)
1044 $l_2 = $LOG_2->copy(); # copy for mul
1048 # else: slower, compute it (but don't cache it, because it could be big)
1049 # also disable downgrade for this code path
1050 local $Math::BigFloat::downgrade = undef;
1051 $l_2 = $two->blog(undef,$scale); # scale+4, actually
1053 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
1056 $self->_log($x,$scale); # need to do the "normal" way
1057 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
1058 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
1059 # all done, $x contains now the result
1064 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1065 # does not modify arguments, but returns new object
1066 # Lowest Common Multiplicator
1068 my ($self,@arg) = objectify(0,@_);
1069 my $x = $self->new(shift @arg);
1070 while (@arg) { $x = _lcm($x,shift @arg); }
1076 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
1077 # does not modify arguments, but returns new object
1078 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
1080 my ($self,@arg) = objectify(0,@_);
1081 my $x = $self->new(shift @arg);
1082 while (@arg) { $x = _gcd($x,shift @arg); }
1086 ##############################################################################
1090 # Internal helper sub to take two positive integers and their signs and
1091 # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1092 # output ($CALC,('+'|'-'))
1093 my ($x,$y,$xs,$ys) = @_;
1095 # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1098 $x = $MBI->_add ($x, $y ); # a+b
1099 # the sign follows $xs
1103 my $a = $MBI->_acmp($x,$y);
1106 $x = $MBI->_sub ($x , $y); # abs sub
1110 $x = $MBI->_zero(); # result is 0
1115 $x = $MBI->_sub ( $y, $x, 1 ); # abs sub
1123 # Internal helper sub to take two positive integers and their signs and
1124 # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1125 # output ($CALC,('+'|'-'))
1126 my ($x,$y,$xs,$ys) = @_;
1130 _e_add($x,$y,$xs,$ys); # call add (does subtract now)
1133 ###############################################################################
1134 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1138 # return true if arg (BFLOAT or num_str) is an integer
1139 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1141 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1142 $x->{_es} eq '+'; # 1e-1 => no integer
1148 # return true if arg (BFLOAT or num_str) is zero
1149 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1151 return 1 if $x->{sign} eq '+' && $MBI->_is_zero($x->{_m});
1157 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1158 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1160 $sign = '+' if !defined $sign || $sign ne '-';
1162 if ($x->{sign} eq $sign &&
1163 $MBI->_is_zero($x->{_e}) && $MBI->_is_one($x->{_m}));
1169 # return true if arg (BFLOAT or num_str) is odd or false if even
1170 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1172 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1173 ($MBI->_is_zero($x->{_e}) && $MBI->_is_odd($x->{_m}));
1179 # return true if arg (BINT or num_str) is even or false if odd
1180 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1182 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
1183 return 1 if ($x->{_es} eq '+' # 123.45 is never
1184 && $MBI->_is_even($x->{_m})); # but 1200 is
1190 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
1191 # (BINT or num_str, BINT or num_str) return BINT
1194 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1195 # objectify is costly, so avoid it
1196 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1198 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1201 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1204 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1206 return $x->bnan() if $x->is_zero() || $y->is_zero();
1207 # result will always be +-inf:
1208 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1209 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1210 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1211 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1212 return $x->binf('-');
1215 return $x->bzero() if $x->is_zero() || $y->is_zero();
1217 return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
1218 ((!$x->isa($self)) || (!$y->isa($self)));
1220 # aEb * cEd = (a*c)E(b+d)
1221 $MBI->_mul($x->{_m},$y->{_m});
1222 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1225 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1226 return $x->bnorm()->round($a,$p,$r,$y);
1231 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1232 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1235 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1236 # objectify is costly, so avoid it
1237 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1239 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1242 return $self->_div_inf($x,$y)
1243 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1245 # x== 0 # also: or y == 1 or y == -1
1246 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1249 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1251 # we need to limit the accuracy to protect against overflow
1253 my (@params,$scale);
1254 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1256 return $x if $x->is_nan(); # error in _find_round_parameters?
1258 # no rounding at all, so must use fallback
1259 if (scalar @params == 0)
1261 # simulate old behaviour
1262 $params[0] = $self->div_scale(); # and round to it as accuracy
1263 $scale = $params[0]+4; # at least four more for proper round
1264 $params[2] = $r; # round mode by caller or undef
1265 $fallback = 1; # to clear a/p afterwards
1269 # the 4 below is empirical, and there might be cases where it is not
1271 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1274 my $rem; $rem = $self->bzero() if wantarray;
1276 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1278 my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
1279 $scale = $lx if $lx > $scale;
1280 $scale = $ly if $ly > $scale;
1281 my $diff = $ly - $lx;
1282 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
1284 # cases like $x /= $x (but not $x /= $y!) were wrong due to modifying $x
1286 require Scalar::Util;
1287 if (Scalar::Util::refaddr($x) == Scalar::Util::refaddr($y))
1289 $x->bone(); # x/x => 1, rem 0
1294 # make copy of $x in case of list context for later reminder calculation
1295 if (wantarray && !$y->is_one())
1300 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1302 # check for / +-1 ( +/- 1E0)
1305 # promote BigInts and it's subclasses (except when already a BigFloat)
1306 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1308 # calculate the result to $scale digits and then round it
1309 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1310 $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1311 $MBI->_div ($x->{_m},$y->{_m}); # a/c
1313 # correct exponent of $x
1314 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1315 # correct for 10**scale
1316 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1317 $x->bnorm(); # remove trailing 0's
1319 } # ende else $x != $y
1321 # shortcut to not run through _find_round_parameters again
1322 if (defined $params[0])
1324 delete $x->{_a}; # clear before round
1325 $x->bround($params[0],$params[2]); # then round accordingly
1329 delete $x->{_p}; # clear before round
1330 $x->bfround($params[1],$params[2]); # then round accordingly
1334 # clear a/p after round, since user did not request it
1335 delete $x->{_a}; delete $x->{_p};
1342 $rem->bmod($y,@params); # copy already done
1346 # clear a/p after round, since user did not request it
1347 delete $rem->{_a}; delete $rem->{_p};
1356 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1359 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1360 # objectify is costly, so avoid it
1361 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1363 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1366 # handle NaN, inf, -inf
1367 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1369 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1370 $x->{sign} = $re->{sign};
1371 $x->{_e} = $re->{_e};
1372 $x->{_m} = $re->{_m};
1373 return $x->round($a,$p,$r,$y);
1377 return $x->bnan() if $x->is_zero();
1380 return $x->bzero() if $y->is_one() || $x->is_zero();
1382 my $cmp = $x->bacmp($y); # equal or $x < $y?
1383 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1385 # only $y of the operands negative?
1386 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1388 $x->{sign} = $y->{sign}; # calc sign first
1389 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1391 my $ym = $MBI->_copy($y->{_m});
1394 $MBI->_lsft( $ym, $y->{_e}, 10)
1395 if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1397 # if $y has digits after dot
1398 my $shifty = 0; # correct _e of $x by this
1399 if ($y->{_es} eq '-') # has digits after dot
1401 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1402 $shifty = $MBI->_num($y->{_e}); # no more digits after dot
1403 $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1405 # $ym is now mantissa of $y based on exponent 0
1407 my $shiftx = 0; # correct _e of $x by this
1408 if ($x->{_es} eq '-') # has digits after dot
1410 # 123.4 % 20 => 1234 % 200
1411 $shiftx = $MBI->_num($x->{_e}); # no more digits after dot
1412 $MBI->_lsft($ym, $x->{_e}, 10); # 123 => 1230
1414 # 123e1 % 20 => 1230 % 20
1415 if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1417 $MBI->_lsft( $x->{_m}, $x->{_e},10); # es => '+' here
1420 $x->{_e} = $MBI->_new($shiftx);
1422 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1423 $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1425 # now mantissas are equalized, exponent of $x is adjusted, so calc result
1427 $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1429 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1432 if ($neg != 0) # one of them negative => correct in place
1435 $x->{_m} = $r->{_m};
1436 $x->{_e} = $r->{_e};
1437 $x->{_es} = $r->{_es};
1438 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1442 $x->round($a,$p,$r,$y); # round and return
1447 # calculate $y'th root of $x
1450 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1451 # objectify is costly, so avoid it
1452 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1454 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1457 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1458 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1459 $y->{sign} !~ /^\+$/;
1461 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1463 # we need to limit the accuracy to protect against overflow
1465 my (@params,$scale);
1466 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1468 return $x if $x->is_nan(); # error in _find_round_parameters?
1470 # no rounding at all, so must use fallback
1471 if (scalar @params == 0)
1473 # simulate old behaviour
1474 $params[0] = $self->div_scale(); # and round to it as accuracy
1475 $scale = $params[0]+4; # at least four more for proper round
1476 $params[2] = $r; # iound mode by caller or undef
1477 $fallback = 1; # to clear a/p afterwards
1481 # the 4 below is empirical, and there might be cases where it is not
1483 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1486 # when user set globals, they would interfere with our calculation, so
1487 # disable them and later re-enable them
1489 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1490 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1491 # we also need to disable any set A or P on $x (_find_round_parameters took
1492 # them already into account), since these would interfere, too
1493 delete $x->{_a}; delete $x->{_p};
1494 # need to disable $upgrade in BigInt, to avoid deep recursion
1495 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1497 # remember sign and make $x positive, since -4 ** (1/2) => -2
1498 my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
1501 if ($y->isa('Math::BigFloat'))
1503 $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
1507 $is_two = ($y == 2);
1510 # normal square root if $y == 2:
1513 $x->bsqrt($scale+4);
1515 elsif ($y->is_one('-'))
1518 my $u = $self->bone()->bdiv($x,$scale);
1519 # copy private parts over
1520 $x->{_m} = $u->{_m};
1521 $x->{_e} = $u->{_e};
1522 $x->{_es} = $u->{_es};
1526 # calculate the broot() as integer result first, and if it fits, return
1527 # it rightaway (but only if $x and $y are integer):
1529 my $done = 0; # not yet
1530 if ($y->is_int() && $x->is_int())
1532 my $i = $MBI->_copy( $x->{_m} );
1533 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1534 my $int = Math::BigInt->bzero();
1536 $int->broot($y->as_number());
1538 if ($int->copy()->bpow($y) == $x)
1540 # found result, return it
1541 $x->{_m} = $int->{value};
1542 $x->{_e} = $MBI->_zero();
1550 my $u = $self->bone()->bdiv($y,$scale+4);
1551 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
1552 $x->bpow($u,$scale+4); # el cheapo
1555 $x->bneg() if $sign == 1;
1557 # shortcut to not run through _find_round_parameters again
1558 if (defined $params[0])
1560 $x->bround($params[0],$params[2]); # then round accordingly
1564 $x->bfround($params[1],$params[2]); # then round accordingly
1568 # clear a/p after round, since user did not request it
1569 delete $x->{_a}; delete $x->{_p};
1572 $$abr = $ab; $$pbr = $pb;
1578 # calculate square root
1579 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1581 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
1582 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
1583 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
1585 # we need to limit the accuracy to protect against overflow
1587 my (@params,$scale);
1588 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1590 return $x if $x->is_nan(); # error in _find_round_parameters?
1592 # no rounding at all, so must use fallback
1593 if (scalar @params == 0)
1595 # simulate old behaviour
1596 $params[0] = $self->div_scale(); # and round to it as accuracy
1597 $scale = $params[0]+4; # at least four more for proper round
1598 $params[2] = $r; # round mode by caller or undef
1599 $fallback = 1; # to clear a/p afterwards
1603 # the 4 below is empirical, and there might be cases where it is not
1605 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1608 # when user set globals, they would interfere with our calculation, so
1609 # disable them and later re-enable them
1611 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1612 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1613 # we also need to disable any set A or P on $x (_find_round_parameters took
1614 # them already into account), since these would interfere, too
1615 delete $x->{_a}; delete $x->{_p};
1616 # need to disable $upgrade in BigInt, to avoid deep recursion
1617 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1619 my $i = $MBI->_copy( $x->{_m} );
1620 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1621 my $xas = Math::BigInt->bzero();
1624 my $gs = $xas->copy()->bsqrt(); # some guess
1626 if (($x->{_es} ne '-') # guess can't be accurate if there are
1627 # digits after the dot
1628 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
1630 # exact result, copy result over to keep $x
1631 $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
1633 # shortcut to not run through _find_round_parameters again
1634 if (defined $params[0])
1636 $x->bround($params[0],$params[2]); # then round accordingly
1640 $x->bfround($params[1],$params[2]); # then round accordingly
1644 # clear a/p after round, since user did not request it
1645 delete $x->{_a}; delete $x->{_p};
1647 # re-enable A and P, upgrade is taken care of by "local"
1648 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1652 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1653 # of the result by multipyling the input by 100 and then divide the integer
1654 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
1656 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
1657 my $y1 = $MBI->_copy($x->{_m});
1659 my $length = $MBI->_len($y1);
1661 # Now calculate how many digits the result of sqrt(y1) would have
1662 my $digits = int($length / 2);
1664 # But we need at least $scale digits, so calculate how many are missing
1665 my $shift = $scale - $digits;
1667 # That should never happen (we take care of integer guesses above)
1668 # $shift = 0 if $shift < 0;
1670 # Multiply in steps of 100, by shifting left two times the "missing" digits
1671 my $s2 = $shift * 2;
1673 # We now make sure that $y1 has the same odd or even number of digits than
1674 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1675 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1676 # steps of 10. The length of $x does not count, since an even or odd number
1677 # of digits before the dot is not changed by adding an even number of digits
1678 # after the dot (the result is still odd or even digits long).
1679 $s2++ if $MBI->_is_odd($x->{_e});
1681 $MBI->_lsft( $y1, $MBI->_new($s2), 10);
1683 # now take the square root and truncate to integer
1684 $y1 = $MBI->_sqrt($y1);
1686 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1687 # result, which is than later rounded to the desired scale.
1689 # calculate how many zeros $x had after the '.' (or before it, depending
1690 # on sign of $dat, the result should have half as many:
1691 my $dat = $MBI->_num($x->{_e});
1692 $dat = -$dat if $x->{_es} eq '-';
1697 # no zeros after the dot (e.g. 1.23, 0.49 etc)
1698 # preserve half as many digits before the dot than the input had
1699 # (but round this "up")
1700 $dat = int(($dat+1)/2);
1704 $dat = int(($dat)/2);
1706 $dat -= $MBI->_len($y1);
1710 $x->{_e} = $MBI->_new( $dat );
1715 $x->{_e} = $MBI->_new( $dat );
1721 # shortcut to not run through _find_round_parameters again
1722 if (defined $params[0])
1724 $x->bround($params[0],$params[2]); # then round accordingly
1728 $x->bfround($params[1],$params[2]); # then round accordingly
1732 # clear a/p after round, since user did not request it
1733 delete $x->{_a}; delete $x->{_p};
1736 $$abr = $ab; $$pbr = $pb;
1742 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1743 # compute factorial number, modifies first argument
1746 my ($self,$x,@r) = (ref($_[0]),@_);
1747 # objectify is costly, so avoid it
1748 ($self,$x,@r) = objectify(1,@_) if !ref($x);
1750 return $x if $x->{sign} eq '+inf'; # inf => inf
1752 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1753 ($x->{_es} ne '+')); # digits after dot?
1755 # use BigInt's bfac() for faster calc
1756 if (! $MBI->_is_zero($x->{_e}))
1758 $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
1759 $x->{_e} = $MBI->_zero(); # normalize
1762 $MBI->_fac($x->{_m}); # calculate factorial
1763 $x->bnorm()->round(@r); # norm again and round result
1768 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1769 my ($x,$y,$a,$p,$r) = @_;
1772 # if $y == 0.5, it is sqrt($x)
1773 $HALF = $self->new($HALF) unless ref($HALF);
1774 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp($HALF) == 0;
1777 # a ** x == e ** (x * ln a)
1781 # Taylor: | u u^2 u^3 |
1782 # x ** y = 1 + | --- + --- + ----- + ... |
1785 # we need to limit the accuracy to protect against overflow
1787 my ($scale,@params);
1788 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1790 return $x if $x->is_nan(); # error in _find_round_parameters?
1792 # no rounding at all, so must use fallback
1793 if (scalar @params == 0)
1795 # simulate old behaviour
1796 $params[0] = $self->div_scale(); # and round to it as accuracy
1797 $params[1] = undef; # disable P
1798 $scale = $params[0]+4; # at least four more for proper round
1799 $params[2] = $r; # round mode by caller or undef
1800 $fallback = 1; # to clear a/p afterwards
1804 # the 4 below is empirical, and there might be cases where it is not
1806 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1809 # when user set globals, they would interfere with our calculation, so
1810 # disable them and later re-enable them
1812 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1813 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1814 # we also need to disable any set A or P on $x (_find_round_parameters took
1815 # them already into account), since these would interfere, too
1816 delete $x->{_a}; delete $x->{_p};
1817 # need to disable $upgrade in BigInt, to avoid deep recursion
1818 local $Math::BigInt::upgrade = undef;
1820 my ($limit,$v,$u,$below,$factor,$next,$over);
1822 $u = $x->copy()->blog(undef,$scale)->bmul($y);
1823 $v = $self->bone(); # 1
1824 $factor = $self->new(2); # 2
1825 $x->bone(); # first term: 1
1827 $below = $v->copy();
1830 $limit = $self->new("1E-". ($scale-1));
1834 # we calculate the next term, and add it to the last
1835 # when the next term is below our limit, it won't affect the outcome
1836 # anymore, so we stop
1837 $next = $over->copy()->bdiv($below,$scale);
1838 last if $next->bacmp($limit) <= 0;
1840 # calculate things for the next term
1841 $over *= $u; $below *= $factor; $factor->binc();
1843 last if $x->{sign} !~ /^[-+]$/;
1848 # shortcut to not run through _find_round_parameters again
1849 if (defined $params[0])
1851 $x->bround($params[0],$params[2]); # then round accordingly
1855 $x->bfround($params[1],$params[2]); # then round accordingly
1859 # clear a/p after round, since user did not request it
1860 delete $x->{_a}; delete $x->{_p};
1863 $$abr = $ab; $$pbr = $pb;
1869 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1870 # compute power of two numbers, second arg is used as integer
1871 # modifies first argument
1874 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1875 # objectify is costly, so avoid it
1876 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1878 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1881 return $x if $x->{sign} =~ /^[+-]inf$/;
1882 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1884 # cache the result of is_zero
1885 my $y_is_zero = $y->is_zero();
1886 return $x->bone() if $y_is_zero;
1887 return $x if $x->is_one() || $y->is_one();
1889 my $x_is_zero = $x->is_zero();
1890 return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int(); # non-integer power
1892 my $y1 = $y->as_number()->{value}; # make MBI part
1895 if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
1897 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
1898 return $MBI->_is_odd($y1) ? $x : $x->babs(1);
1902 return $x->bone() if $y_is_zero;
1903 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1904 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
1909 $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
1911 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1912 $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
1913 $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
1915 $x->{sign} = $new_sign;
1917 if ($y->{sign} eq '-')
1919 # modify $x in place!
1920 my $z = $x->copy(); $x->bone();
1921 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
1923 $x->round($a,$p,$r,$y);
1926 ###############################################################################
1927 # rounding functions
1931 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1932 # $n == 0 means round to integer
1933 # expects and returns normalized numbers!
1934 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1936 return $x if $x->modify('bfround');
1938 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1939 return $x if !defined $scale; # no-op
1941 # never round a 0, +-inf, NaN
1944 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1947 return $x if $x->{sign} !~ /^[+-]$/;
1949 # don't round if x already has lower precision
1950 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1952 $x->{_p} = $scale; # remember round in any case
1953 delete $x->{_a}; # and clear A
1956 # round right from the '.'
1958 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round
1960 $scale = -$scale; # positive for simplicity
1961 my $len = $MBI->_len($x->{_m}); # length of mantissa
1963 # the following poses a restriction on _e, but if _e is bigger than a
1964 # scalar, you got other problems (memory etc) anyway
1965 my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e}))); # digits after dot
1966 my $zad = 0; # zeros after dot
1967 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
1969 # p rint "scale $scale dad $dad zad $zad len $len\n";
1970 # number bsstr len zad dad
1971 # 0.123 123e-3 3 0 3
1972 # 0.0123 123e-4 3 1 4
1975 # 1.2345 12345e-4 5 0 4
1977 # do not round after/right of the $dad
1978 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1980 # round to zero if rounding inside the $zad, but not for last zero like:
1981 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1982 return $x->bzero() if $scale < $zad;
1983 if ($scale == $zad) # for 0.006, scale -3 and trunc
1989 # adjust round-point to be inside mantissa
1992 $scale = $scale-$zad;
1996 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
1997 $scale = $dbd+$scale;
2003 # round left from the '.'
2005 # 123 => 100 means length(123) = 3 - $scale (2) => 1
2007 my $dbt = $MBI->_len($x->{_m});
2009 my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
2010 # should be the same, so treat it as this
2011 $scale = 1 if $scale == 0;
2012 # shortcut if already integer
2013 return $x if $scale == 1 && $dbt <= $dbd;
2014 # maximum digits before dot
2019 # not enough digits before dot, so round to zero
2022 elsif ( $scale == $dbd )
2029 $scale = $dbd - $scale;
2032 # pass sign to bround for rounding modes '+inf' and '-inf'
2033 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2034 $m->bround($scale,$mode);
2035 $x->{_m} = $m->{value}; # get our mantissa back
2041 # accuracy: preserve $N digits, and overwrite the rest with 0's
2042 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2044 if (($_[0] || 0) < 0)
2046 require Carp; Carp::croak ('bround() needs positive accuracy');
2049 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
2050 return $x if !defined $scale; # no-op
2052 return $x if $x->modify('bround');
2054 # scale is now either $x->{_a}, $accuracy, or the user parameter
2055 # test whether $x already has lower accuracy, do nothing in this case
2056 # but do round if the accuracy is the same, since a math operation might
2057 # want to round a number with A=5 to 5 digits afterwards again
2058 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
2060 # scale < 0 makes no sense
2061 # never round a +-inf, NaN
2062 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
2064 # 1: $scale == 0 => keep all digits
2065 # 2: never round a 0
2066 # 3: if we should keep more digits than the mantissa has, do nothing
2067 if ($scale == 0 || $x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
2069 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
2073 # pass sign to bround for '+inf' and '-inf' rounding modes
2074 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2076 $m->bround($scale,$mode); # round mantissa
2077 $x->{_m} = $m->{value}; # get our mantissa back
2078 $x->{_a} = $scale; # remember rounding
2079 delete $x->{_p}; # and clear P
2080 $x->bnorm(); # del trailing zeros gen. by bround()
2085 # return integer less or equal then $x
2086 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2088 return $x if $x->modify('bfloor');
2090 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2092 # if $x has digits after dot
2093 if ($x->{_es} eq '-')
2095 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2096 $x->{_e} = $MBI->_zero(); # trunc/norm
2097 $x->{_es} = '+'; # abs e
2098 $MBI->_inc($x->{_m}) if $x->{sign} eq '-'; # increment if negative
2100 $x->round($a,$p,$r);
2105 # return integer greater or equal then $x
2106 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2108 return $x if $x->modify('bceil');
2109 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2111 # if $x has digits after dot
2112 if ($x->{_es} eq '-')
2114 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2115 $x->{_e} = $MBI->_zero(); # trunc/norm
2116 $x->{_es} = '+'; # abs e
2117 $MBI->_inc($x->{_m}) if $x->{sign} eq '+'; # increment if positive
2119 $x->round($a,$p,$r);
2124 # shift right by $y (divide by power of $n)
2127 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2128 # objectify is costly, so avoid it
2129 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2131 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2134 return $x if $x->modify('brsft');
2135 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2137 $n = 2 if !defined $n; $n = $self->new($n);
2138 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
2143 # shift left by $y (multiply by power of $n)
2146 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2147 # objectify is costly, so avoid it
2148 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2150 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2153 return $x if $x->modify('blsft');
2154 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2156 $n = 2 if !defined $n; $n = $self->new($n);
2157 $x->bmul($n->bpow($y),$a,$p,$r,$y);
2160 ###############################################################################
2164 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
2169 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
2170 # or falling back to MBI::bxxx()
2171 my $name = $AUTOLOAD;
2173 $name =~ s/(.*):://; # split package
2174 my $c = $1 || $class;
2176 $c->import() if $IMPORT == 0;
2177 if (!method_alias($name))
2181 # delayed load of Carp and avoid recursion
2183 Carp::croak ("$c: Can't call a method without name");
2185 if (!method_hand_up($name))
2187 # delayed load of Carp and avoid recursion
2189 Carp::croak ("Can't call $c\-\>$name, not a valid method");
2191 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
2193 return &{"Math::BigInt"."::$name"}(@_);
2195 my $bname = $name; $bname =~ s/^f/b/;
2203 # return a copy of the exponent
2204 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2206 if ($x->{sign} !~ /^[+-]$/)
2208 my $s = $x->{sign}; $s =~ s/^[+-]//;
2209 return Math::BigInt->new($s); # -inf, +inf => +inf
2211 Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
2216 # return a copy of the mantissa
2217 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2219 if ($x->{sign} !~ /^[+-]$/)
2221 my $s = $x->{sign}; $s =~ s/^[+]//;
2222 return Math::BigInt->new($s); # -inf, +inf => +inf
2224 my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
2225 $m->bneg() if $x->{sign} eq '-';
2232 # return a copy of both the exponent and the mantissa
2233 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2235 if ($x->{sign} !~ /^[+-]$/)
2237 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2238 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2240 my $m = Math::BigInt->bzero();
2241 $m->{value} = $MBI->_copy($x->{_m});
2242 $m->bneg() if $x->{sign} eq '-';
2243 ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
2246 ##############################################################################
2247 # private stuff (internal use only)
2253 my $lib = ''; my @a;
2255 for ( my $i = 0; $i < $l ; $i++)
2257 if ( $_[$i] eq ':constant' )
2259 # This causes overlord er load to step in. 'binary' and 'integer'
2260 # are handled by BigInt.
2261 overload::constant float => sub { $self->new(shift); };
2263 elsif ($_[$i] eq 'upgrade')
2265 # this causes upgrading
2266 $upgrade = $_[$i+1]; # or undef to disable
2269 elsif ($_[$i] eq 'downgrade')
2271 # this causes downgrading
2272 $downgrade = $_[$i+1]; # or undef to disable
2275 elsif ($_[$i] eq 'lib')
2277 # alternative library
2278 $lib = $_[$i+1] || ''; # default Calc
2281 elsif ($_[$i] eq 'with')
2283 # alternative class for our private parts()
2284 # XXX: no longer supported
2285 # $MBI = $_[$i+1] || 'Math::BigInt';
2294 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2295 my $mbilib = eval { Math::BigInt->config()->{lib} };
2296 if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
2298 # MBI already loaded
2299 Math::BigInt->import('lib',"$lib,$mbilib", 'objectify');
2303 # MBI not loaded, or with ne "Math::BigInt::Calc"
2304 $lib .= ",$mbilib" if defined $mbilib;
2305 $lib =~ s/^,//; # don't leave empty
2306 # replacement library can handle lib statement, but also could ignore it
2309 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2310 # used in the same script, or eval inside import().
2311 require Math::BigInt;
2312 Math::BigInt->import( lib => $lib, 'objectify' );
2316 my $rc = "use Math::BigInt lib => '$lib', 'objectify';";
2322 require Carp; Carp::croak ("Couldn't load $lib: $! $@");
2324 $MBI = Math::BigInt->config()->{lib};
2326 # any non :constant stuff is handled by our parent, Exporter
2327 # even if @_ is empty, to give it a chance
2328 $self->SUPER::import(@a); # for subclasses
2329 $self->export_to_level(1,$self,@a); # need this, too
2334 # adjust m and e so that m is smallest possible
2335 # round number according to accuracy and precision settings
2336 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
2338 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2340 my $zeros = $MBI->_zeros($x->{_m}); # correct for trailing zeros
2343 my $z = $MBI->_new($zeros);
2344 $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
2345 if ($x->{_es} eq '-')
2347 if ($MBI->_acmp($x->{_e},$z) >= 0)
2349 $x->{_e} = $MBI->_sub ($x->{_e}, $z);
2350 $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
2354 $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
2360 $x->{_e} = $MBI->_add ($x->{_e}, $z);
2365 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
2366 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
2367 $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
2368 if $MBI->_is_zero($x->{_m});
2371 $x; # MBI bnorm is no-op, so dont call it
2374 ##############################################################################
2378 # return number as hexadecimal string (only for integers defined)
2379 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2381 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2382 return '0x0' if $x->is_zero();
2384 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2386 my $z = $MBI->_copy($x->{_m});
2387 if (! $MBI->_is_zero($x->{_e})) # > 0
2389 $MBI->_lsft( $z, $x->{_e},10);
2391 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2397 # return number as binary digit string (only for integers defined)
2398 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2400 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2401 return '0b0' if $x->is_zero();
2403 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2405 my $z = $MBI->_copy($x->{_m});
2406 if (! $MBI->_is_zero($x->{_e})) # > 0
2408 $MBI->_lsft( $z, $x->{_e},10);
2410 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2416 # return copy as a bigint representation of this BigFloat number
2417 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2419 my $z = $MBI->_copy($x->{_m});
2420 if ($x->{_es} eq '-') # < 0
2422 $MBI->_rsft( $z, $x->{_e},10);
2424 elsif (! $MBI->_is_zero($x->{_e})) # > 0
2426 $MBI->_lsft( $z, $x->{_e},10);
2428 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2435 my $class = ref($x) || $x;
2436 $x = $class->new(shift) unless ref($x);
2438 return 1 if $MBI->_is_zero($x->{_m});
2440 my $len = $MBI->_len($x->{_m});
2441 $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
2445 $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
2456 Math::BigFloat - Arbitrary size floating point math package
2463 $x = Math::BigFloat->new($str); # defaults to 0
2464 $nan = Math::BigFloat->bnan(); # create a NotANumber
2465 $zero = Math::BigFloat->bzero(); # create a +0
2466 $inf = Math::BigFloat->binf(); # create a +inf
2467 $inf = Math::BigFloat->binf('-'); # create a -inf
2468 $one = Math::BigFloat->bone(); # create a +1
2469 $one = Math::BigFloat->bone('-'); # create a -1
2472 $x->is_zero(); # true if arg is +0
2473 $x->is_nan(); # true if arg is NaN
2474 $x->is_one(); # true if arg is +1
2475 $x->is_one('-'); # true if arg is -1
2476 $x->is_odd(); # true if odd, false for even
2477 $x->is_even(); # true if even, false for odd
2478 $x->is_pos(); # true if >= 0
2479 $x->is_neg(); # true if < 0
2480 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2482 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2483 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2484 $x->sign(); # return the sign, either +,- or NaN
2485 $x->digit($n); # return the nth digit, counting from right
2486 $x->digit(-$n); # return the nth digit, counting from left
2488 # The following all modify their first argument. If you want to preserve
2489 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2490 # neccessary when mixing $a = $b assigments with non-overloaded math.
2493 $x->bzero(); # set $i to 0
2494 $x->bnan(); # set $i to NaN
2495 $x->bone(); # set $x to +1
2496 $x->bone('-'); # set $x to -1
2497 $x->binf(); # set $x to inf
2498 $x->binf('-'); # set $x to -inf
2500 $x->bneg(); # negation
2501 $x->babs(); # absolute value
2502 $x->bnorm(); # normalize (no-op)
2503 $x->bnot(); # two's complement (bit wise not)
2504 $x->binc(); # increment x by 1
2505 $x->bdec(); # decrement x by 1
2507 $x->badd($y); # addition (add $y to $x)
2508 $x->bsub($y); # subtraction (subtract $y from $x)
2509 $x->bmul($y); # multiplication (multiply $x by $y)
2510 $x->bdiv($y); # divide, set $x to quotient
2511 # return (quo,rem) or quo if scalar
2513 $x->bmod($y); # modulus ($x % $y)
2514 $x->bpow($y); # power of arguments ($x ** $y)
2515 $x->blsft($y); # left shift
2516 $x->brsft($y); # right shift
2517 # return (quo,rem) or quo if scalar
2519 $x->blog(); # logarithm of $x to base e (Euler's number)
2520 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
2522 $x->band($y); # bit-wise and
2523 $x->bior($y); # bit-wise inclusive or
2524 $x->bxor($y); # bit-wise exclusive or
2525 $x->bnot(); # bit-wise not (two's complement)
2527 $x->bsqrt(); # calculate square-root
2528 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
2529 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2531 $x->bround($N); # accuracy: preserve $N digits
2532 $x->bfround($N); # precision: round to the $Nth digit
2534 $x->bfloor(); # return integer less or equal than $x
2535 $x->bceil(); # return integer greater or equal than $x
2537 # The following do not modify their arguments:
2539 bgcd(@values); # greatest common divisor
2540 blcm(@values); # lowest common multiplicator
2542 $x->bstr(); # return string
2543 $x->bsstr(); # return string in scientific notation
2545 $x->as_int(); # return $x as BigInt
2546 $x->exponent(); # return exponent as BigInt
2547 $x->mantissa(); # return mantissa as BigInt
2548 $x->parts(); # return (mantissa,exponent) as BigInt
2550 $x->length(); # number of digits (w/o sign and '.')
2551 ($l,$f) = $x->length(); # number of digits, and length of fraction
2553 $x->precision(); # return P of $x (or global, if P of $x undef)
2554 $x->precision($n); # set P of $x to $n
2555 $x->accuracy(); # return A of $x (or global, if A of $x undef)
2556 $x->accuracy($n); # set A $x to $n
2558 # these get/set the appropriate global value for all BigFloat objects
2559 Math::BigFloat->precision(); # Precision
2560 Math::BigFloat->accuracy(); # Accuracy
2561 Math::BigFloat->round_mode(); # rounding mode
2565 All operators (inlcuding basic math operations) are overloaded if you
2566 declare your big floating point numbers as
2568 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2570 Operations with overloaded operators preserve the arguments, which is
2571 exactly what you expect.
2573 =head2 Canonical notation
2575 Input to these routines are either BigFloat objects, or strings of the
2576 following four forms:
2590 C</^[+-]\d+E[+-]?\d+$/>
2594 C</^[+-]\d*\.\d+E[+-]?\d+$/>
2598 all with optional leading and trailing zeros and/or spaces. Additonally,
2599 numbers are allowed to have an underscore between any two digits.
2601 Empty strings as well as other illegal numbers results in 'NaN'.
2603 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2604 are always stored in normalized form. On a string, it creates a BigFloat
2609 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2611 The string output will always have leading and trailing zeros stripped and drop
2612 a plus sign. C<bstr()> will give you always the form with a decimal point,
2613 while C<bsstr()> (s for scientific) gives you the scientific notation.
2615 Input bstr() bsstr()
2617 ' -123 123 123' '-123123123' '-123123123E0'
2618 '00.0123' '0.0123' '123E-4'
2619 '123.45E-2' '1.2345' '12345E-4'
2620 '10E+3' '10000' '1E4'
2622 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2623 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2624 return either undef, <0, 0 or >0 and are suited for sort.
2626 Actual math is done by using the class defined with C<with => Class;> (which
2627 defaults to BigInts) to represent the mantissa and exponent.
2629 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2630 represent the result when input arguments are not numbers, as well as
2631 the result of dividing by zero.
2633 =head2 C<mantissa()>, C<exponent()> and C<parts()>
2635 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2636 as BigInts such that:
2638 $m = $x->mantissa();
2639 $e = $x->exponent();
2640 $y = $m * ( 10 ** $e );
2641 print "ok\n" if $x == $y;
2643 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2645 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2647 Currently the mantissa is reduced as much as possible, favouring higher
2648 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2649 This might change in the future, so do not depend on it.
2651 =head2 Accuracy vs. Precision
2653 See also: L<Rounding|Rounding>.
2655 Math::BigFloat supports both precision and accuracy. For a full documentation,
2656 examples and tips on these topics please see the large section in
2659 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2660 a operation consumes all resources, each operation produces no more than
2661 the requested number of digits.
2663 Please refer to BigInt's documentation for the precedence rules of which
2664 accuracy/precision setting will be used.
2666 If there is no gloabl precision set, B<and> the operation inquestion was not
2667 called with a requested precision or accuracy, B<and> the input $x has no
2668 accuracy or precision set, then a fallback parameter will be used. For
2669 historical reasons, it is called C<div_scale> and can be accessed via:
2671 $d = Math::BigFloat->div_scale(); # query
2672 Math::BigFloat->div_scale($n); # set to $n digits
2674 The default value is 40 digits.
2676 In case the result of one operation has more precision than specified,
2677 it is rounded. The rounding mode taken is either the default mode, or the one
2678 supplied to the operation after the I<scale>:
2680 $x = Math::BigFloat->new(2);
2681 Math::BigFloat->precision(5); # 5 digits max
2682 $y = $x->copy()->bdiv(3); # will give 0.66666
2683 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2684 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
2685 Math::BigFloat->round_mode('zero');
2686 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2692 =item ffround ( +$scale )
2694 Rounds to the $scale'th place left from the '.', counting from the dot.
2695 The first digit is numbered 1.
2697 =item ffround ( -$scale )
2699 Rounds to the $scale'th place right from the '.', counting from the dot.
2703 Rounds to an integer.
2705 =item fround ( +$scale )
2707 Preserves accuracy to $scale digits from the left (aka significant digits)
2708 and pads the rest with zeros. If the number is between 1 and -1, the
2709 significant digits count from the first non-zero after the '.'
2711 =item fround ( -$scale ) and fround ( 0 )
2713 These are effectively no-ops.
2717 All rounding functions take as a second parameter a rounding mode from one of
2718 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2720 The default rounding mode is 'even'. By using
2721 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
2722 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2723 no longer supported.
2724 The second parameter to the round functions then overrides the default
2727 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
2728 'trunc' as rounding mode to make it equivalent to:
2733 You can override this by passing the desired rounding mode as parameter to
2736 $x = Math::BigFloat->new(2.5);
2737 $y = $x->as_number('odd'); # $y = 3
2743 =head1 Autocreating constants
2745 After C<use Math::BigFloat ':constant'> all the floating point constants
2746 in the given scope are converted to C<Math::BigFloat>. This conversion
2747 happens at compile time.
2751 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2753 prints the value of C<2E-100>. Note that without conversion of
2754 constants the expression 2E-100 will be calculated as normal floating point
2757 Please note that ':constant' does not affect integer constants, nor binary
2758 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2763 Math with the numbers is done (by default) by a module called
2764 Math::BigInt::Calc. This is equivalent to saying:
2766 use Math::BigFloat lib => 'Calc';
2768 You can change this by using:
2770 use Math::BigFloat lib => 'BitVect';
2772 The following would first try to find Math::BigInt::Foo, then
2773 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2775 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2777 Calc.pm uses as internal format an array of elements of some decimal base
2778 (usually 1e7, but this might be differen for some systems) with the least
2779 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2780 significant bit first. Other modules might use even different means of
2781 representing the numbers. See the respective module documentation for further
2784 Please note that Math::BigFloat does B<not> use the denoted library itself,
2785 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2788 use Math::BigInt lib => 'GMP';
2791 you can roll it all into one line:
2793 use Math::BigFloat lib => 'GMP';
2795 It is also possible to just require Math::BigFloat:
2797 require Math::BigFloat;
2799 This will load the neccessary things (like BigInt) when they are needed, and
2802 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2803 you ever wanted to know about loading a different library.
2805 =head2 Using Math::BigInt::Lite
2807 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2810 use Math::BigFloat with => 'Math::BigInt::Lite';
2812 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2813 can combine these if you want. For instance, you may want to use
2814 Math::BigInt objects in your main script, too.
2818 use Math::BigFloat with => 'Math::BigInt::Lite';
2820 Of course, you can combine this with the C<lib> parameter.
2823 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2825 There is no need for a "use Math::BigInt;" statement, even if you want to
2826 use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2827 always loads it. But if you add it, add it B<before>:
2831 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2833 Notice that the module with the last C<lib> will "win" and thus
2834 it's lib will be used if the lib is available:
2837 use Math::BigInt lib => 'Bar,Baz';
2838 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2840 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2841 words, Math::BigFloat will try to retain previously loaded libs when you
2842 don't specify it onem but if you specify one, it will try to load them.
2844 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2845 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2846 same as trying the latter load alone, except for the fact that one of Bar or
2847 Baz might be loaded needlessly in an intermidiate step (and thus hang around
2848 and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2849 will still be tried to be loaded, but this is not as time/memory consuming as
2850 actually loading one of them. Still, this type of usage is not recommended due
2853 The old way (loading the lib only in BigInt) still works though:
2856 use Math::BigInt lib => 'Bar,Baz';
2859 You can even load Math::BigInt afterwards:
2863 use Math::BigInt lib => 'Bar,Baz';
2865 But this has the same problems like #5, it will first load Calc
2866 (Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2867 Baz, depending on which of them works and is usable/loadable. Since this
2868 loads Calc unnecc., it is not recommended.
2870 Since it also possible to just require Math::BigFloat, this poses the question
2871 about what libary this will use:
2873 require Math::BigFloat;
2874 my $x = Math::BigFloat->new(123); $x += 123;
2876 It will use Calc. Please note that the call to import() is still done, but
2877 only when you use for the first time some Math::BigFloat math (it is triggered
2878 via any constructor, so the first time you create a Math::BigFloat, the load
2879 will happen in the background). This means:
2881 require Math::BigFloat;
2882 Math::BigFloat->import ( lib => 'Foo,Bar' );
2884 would be the same as:
2886 use Math::BigFloat lib => 'Foo, Bar';
2888 But don't try to be clever to insert some operations in between:
2890 require Math::BigFloat;
2891 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
2892 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
2893 $x = Math::BigFloat->bone()+4; # now use Pari
2895 While this works, it loads Calc needlessly. But maybe you just wanted that?
2897 B<Examples #3 is highly recommended> for daily usage.
2901 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
2907 =item stringify, bstr()
2909 Both stringify and bstr() now drop the leading '+'. The old code would return
2910 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2911 reasoning and details.
2915 The following will probably not do what you expect:
2917 print $c->bdiv(123.456),"\n";
2919 It prints both quotient and reminder since print works in list context. Also,
2920 bdiv() will modify $c, so be carefull. You probably want to use
2922 print $c / 123.456,"\n";
2923 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2927 =item Modifying and =
2931 $x = Math::BigFloat->new(5);
2934 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2935 a second reference to the B<same> object and stores it in $y. Thus anything
2936 that modifies $x will modify $y (except overloaded math operators), and vice
2937 versa. See L<Math::BigInt> for details and how to avoid that.
2941 C<bpow()> now modifies the first argument, unlike the old code which left
2942 it alone and only returned the result. This is to be consistent with
2943 C<badd()> etc. The first will modify $x, the second one won't:
2945 print bpow($x,$i),"\n"; # modify $x
2946 print $x->bpow($i),"\n"; # ditto
2947 print $x ** $i,"\n"; # leave $x alone
2953 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
2954 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
2956 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
2957 because they solve the autoupgrading/downgrading issue, at least partly.
2960 L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
2961 more documentation including a full version history, testcases, empty
2962 subclass files and benchmarks.
2966 This program is free software; you may redistribute it and/or modify it under
2967 the same terms as Perl itself.
2971 Mark Biggar, overloaded interface by Ilya Zakharevich.
2972 Completely rewritten by Tels http://bloodgate.com in 2001, 2002, and still