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 # already handled inf/NaN/-inf above:
1346 my $xsign = $x->{sign};
1347 $y->{sign} =~ tr/+-/-+/;
1348 my $y_not_one = !$y->is_one(); # cache this result
1349 if ($xsign ne $x->{sign})
1351 # special case of $x /= $x results in 1
1356 # correct $y's sign again
1357 $y->{sign} =~ tr/+-/-+/;
1358 # continue with normal div code:
1360 # make copy of $x in case of list context for later reminder calculation
1361 if (wantarray && $y_not_one)
1366 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1368 # check for / +-1 ( +/- 1E0)
1371 # promote BigInts and it's subclasses (except when already a BigFloat)
1372 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1374 # calculate the result to $scale digits and then round it
1375 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1376 $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1377 $MBI->_div ($x->{_m},$y->{_m}); # a/c
1379 # correct exponent of $x
1380 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1381 # correct for 10**scale
1382 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1383 $x->bnorm(); # remove trailing 0's
1385 } # ende else $x != $y
1387 # shortcut to not run through _find_round_parameters again
1388 if (defined $params[0])
1390 delete $x->{_a}; # clear before round
1391 $x->bround($params[0],$params[2]); # then round accordingly
1395 delete $x->{_p}; # clear before round
1396 $x->bfround($params[1],$params[2]); # then round accordingly
1400 # clear a/p after round, since user did not request it
1401 delete $x->{_a}; delete $x->{_p};
1408 $rem->bmod($y,@params); # copy already done
1412 # clear a/p after round, since user did not request it
1413 delete $rem->{_a}; delete $rem->{_p};
1422 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1425 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1426 # objectify is costly, so avoid it
1427 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1429 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1432 # handle NaN, inf, -inf
1433 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1435 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1436 $x->{sign} = $re->{sign};
1437 $x->{_e} = $re->{_e};
1438 $x->{_m} = $re->{_m};
1439 return $x->round($a,$p,$r,$y);
1443 return $x->bnan() if $x->is_zero();
1446 return $x->bzero() if $y->is_one() || $x->is_zero();
1448 my $cmp = $x->bacmp($y); # equal or $x < $y?
1449 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1451 # only $y of the operands negative?
1452 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1454 $x->{sign} = $y->{sign}; # calc sign first
1455 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1457 my $ym = $MBI->_copy($y->{_m});
1460 $MBI->_lsft( $ym, $y->{_e}, 10)
1461 if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1463 # if $y has digits after dot
1464 my $shifty = 0; # correct _e of $x by this
1465 if ($y->{_es} eq '-') # has digits after dot
1467 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1468 $shifty = $MBI->_num($y->{_e}); # no more digits after dot
1469 $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1471 # $ym is now mantissa of $y based on exponent 0
1473 my $shiftx = 0; # correct _e of $x by this
1474 if ($x->{_es} eq '-') # has digits after dot
1476 # 123.4 % 20 => 1234 % 200
1477 $shiftx = $MBI->_num($x->{_e}); # no more digits after dot
1478 $MBI->_lsft($ym, $x->{_e}, 10); # 123 => 1230
1480 # 123e1 % 20 => 1230 % 20
1481 if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1483 $MBI->_lsft( $x->{_m}, $x->{_e},10); # es => '+' here
1486 $x->{_e} = $MBI->_new($shiftx);
1488 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1489 $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1491 # now mantissas are equalized, exponent of $x is adjusted, so calc result
1493 $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1495 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1498 if ($neg != 0) # one of them negative => correct in place
1501 $x->{_m} = $r->{_m};
1502 $x->{_e} = $r->{_e};
1503 $x->{_es} = $r->{_es};
1504 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1508 $x->round($a,$p,$r,$y); # round and return
1513 # calculate $y'th root of $x
1516 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1517 # objectify is costly, so avoid it
1518 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1520 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1523 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1524 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1525 $y->{sign} !~ /^\+$/;
1527 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1529 # we need to limit the accuracy to protect against overflow
1531 my (@params,$scale);
1532 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1534 return $x if $x->is_nan(); # error in _find_round_parameters?
1536 # no rounding at all, so must use fallback
1537 if (scalar @params == 0)
1539 # simulate old behaviour
1540 $params[0] = $self->div_scale(); # and round to it as accuracy
1541 $scale = $params[0]+4; # at least four more for proper round
1542 $params[2] = $r; # iound mode by caller or undef
1543 $fallback = 1; # to clear a/p afterwards
1547 # the 4 below is empirical, and there might be cases where it is not
1549 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1552 # when user set globals, they would interfere with our calculation, so
1553 # disable them and later re-enable them
1555 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1556 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1557 # we also need to disable any set A or P on $x (_find_round_parameters took
1558 # them already into account), since these would interfere, too
1559 delete $x->{_a}; delete $x->{_p};
1560 # need to disable $upgrade in BigInt, to avoid deep recursion
1561 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1563 # remember sign and make $x positive, since -4 ** (1/2) => -2
1564 my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
1567 if ($y->isa('Math::BigFloat'))
1569 $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
1573 $is_two = ($y == 2);
1576 # normal square root if $y == 2:
1579 $x->bsqrt($scale+4);
1581 elsif ($y->is_one('-'))
1584 my $u = $self->bone()->bdiv($x,$scale);
1585 # copy private parts over
1586 $x->{_m} = $u->{_m};
1587 $x->{_e} = $u->{_e};
1588 $x->{_es} = $u->{_es};
1592 # calculate the broot() as integer result first, and if it fits, return
1593 # it rightaway (but only if $x and $y are integer):
1595 my $done = 0; # not yet
1596 if ($y->is_int() && $x->is_int())
1598 my $i = $MBI->_copy( $x->{_m} );
1599 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1600 my $int = Math::BigInt->bzero();
1602 $int->broot($y->as_number());
1604 if ($int->copy()->bpow($y) == $x)
1606 # found result, return it
1607 $x->{_m} = $int->{value};
1608 $x->{_e} = $MBI->_zero();
1616 my $u = $self->bone()->bdiv($y,$scale+4);
1617 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
1618 $x->bpow($u,$scale+4); # el cheapo
1621 $x->bneg() if $sign == 1;
1623 # shortcut to not run through _find_round_parameters again
1624 if (defined $params[0])
1626 $x->bround($params[0],$params[2]); # then round accordingly
1630 $x->bfround($params[1],$params[2]); # then round accordingly
1634 # clear a/p after round, since user did not request it
1635 delete $x->{_a}; delete $x->{_p};
1638 $$abr = $ab; $$pbr = $pb;
1644 # calculate square root
1645 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1647 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
1648 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
1649 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
1651 # we need to limit the accuracy to protect against overflow
1653 my (@params,$scale);
1654 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1656 return $x if $x->is_nan(); # error in _find_round_parameters?
1658 # no rounding at all, so must use fallback
1659 if (scalar @params == 0)
1661 # simulate old behaviour
1662 $params[0] = $self->div_scale(); # and round to it as accuracy
1663 $scale = $params[0]+4; # at least four more for proper round
1664 $params[2] = $r; # round mode by caller or undef
1665 $fallback = 1; # to clear a/p afterwards
1669 # the 4 below is empirical, and there might be cases where it is not
1671 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1674 # when user set globals, they would interfere with our calculation, so
1675 # disable them and later re-enable them
1677 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1678 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1679 # we also need to disable any set A or P on $x (_find_round_parameters took
1680 # them already into account), since these would interfere, too
1681 delete $x->{_a}; delete $x->{_p};
1682 # need to disable $upgrade in BigInt, to avoid deep recursion
1683 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1685 my $i = $MBI->_copy( $x->{_m} );
1686 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1687 my $xas = Math::BigInt->bzero();
1690 my $gs = $xas->copy()->bsqrt(); # some guess
1692 if (($x->{_es} ne '-') # guess can't be accurate if there are
1693 # digits after the dot
1694 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
1696 # exact result, copy result over to keep $x
1697 $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
1699 # shortcut to not run through _find_round_parameters again
1700 if (defined $params[0])
1702 $x->bround($params[0],$params[2]); # then round accordingly
1706 $x->bfround($params[1],$params[2]); # then round accordingly
1710 # clear a/p after round, since user did not request it
1711 delete $x->{_a}; delete $x->{_p};
1713 # re-enable A and P, upgrade is taken care of by "local"
1714 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1718 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1719 # of the result by multipyling the input by 100 and then divide the integer
1720 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
1722 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
1723 my $y1 = $MBI->_copy($x->{_m});
1725 my $length = $MBI->_len($y1);
1727 # Now calculate how many digits the result of sqrt(y1) would have
1728 my $digits = int($length / 2);
1730 # But we need at least $scale digits, so calculate how many are missing
1731 my $shift = $scale - $digits;
1733 # That should never happen (we take care of integer guesses above)
1734 # $shift = 0 if $shift < 0;
1736 # Multiply in steps of 100, by shifting left two times the "missing" digits
1737 my $s2 = $shift * 2;
1739 # We now make sure that $y1 has the same odd or even number of digits than
1740 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1741 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1742 # steps of 10. The length of $x does not count, since an even or odd number
1743 # of digits before the dot is not changed by adding an even number of digits
1744 # after the dot (the result is still odd or even digits long).
1745 $s2++ if $MBI->_is_odd($x->{_e});
1747 $MBI->_lsft( $y1, $MBI->_new($s2), 10);
1749 # now take the square root and truncate to integer
1750 $y1 = $MBI->_sqrt($y1);
1752 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1753 # result, which is than later rounded to the desired scale.
1755 # calculate how many zeros $x had after the '.' (or before it, depending
1756 # on sign of $dat, the result should have half as many:
1757 my $dat = $MBI->_num($x->{_e});
1758 $dat = -$dat if $x->{_es} eq '-';
1763 # no zeros after the dot (e.g. 1.23, 0.49 etc)
1764 # preserve half as many digits before the dot than the input had
1765 # (but round this "up")
1766 $dat = int(($dat+1)/2);
1770 $dat = int(($dat)/2);
1772 $dat -= $MBI->_len($y1);
1776 $x->{_e} = $MBI->_new( $dat );
1781 $x->{_e} = $MBI->_new( $dat );
1787 # shortcut to not run through _find_round_parameters again
1788 if (defined $params[0])
1790 $x->bround($params[0],$params[2]); # then round accordingly
1794 $x->bfround($params[1],$params[2]); # then round accordingly
1798 # clear a/p after round, since user did not request it
1799 delete $x->{_a}; delete $x->{_p};
1802 $$abr = $ab; $$pbr = $pb;
1808 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1809 # compute factorial number, modifies first argument
1812 my ($self,$x,@r) = (ref($_[0]),@_);
1813 # objectify is costly, so avoid it
1814 ($self,$x,@r) = objectify(1,@_) if !ref($x);
1816 return $x if $x->{sign} eq '+inf'; # inf => inf
1818 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1819 ($x->{_es} ne '+')); # digits after dot?
1821 # use BigInt's bfac() for faster calc
1822 if (! $MBI->_is_zero($x->{_e}))
1824 $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
1825 $x->{_e} = $MBI->_zero(); # normalize
1828 $MBI->_fac($x->{_m}); # calculate factorial
1829 $x->bnorm()->round(@r); # norm again and round result
1834 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1835 my ($x,$y,$a,$p,$r) = @_;
1838 # if $y == 0.5, it is sqrt($x)
1839 $HALF = $self->new($HALF) unless ref($HALF);
1840 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp($HALF) == 0;
1843 # a ** x == e ** (x * ln a)
1847 # Taylor: | u u^2 u^3 |
1848 # x ** y = 1 + | --- + --- + ----- + ... |
1851 # we need to limit the accuracy to protect against overflow
1853 my ($scale,@params);
1854 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1856 return $x if $x->is_nan(); # error in _find_round_parameters?
1858 # no rounding at all, so must use fallback
1859 if (scalar @params == 0)
1861 # simulate old behaviour
1862 $params[0] = $self->div_scale(); # and round to it as accuracy
1863 $params[1] = undef; # disable P
1864 $scale = $params[0]+4; # at least four more for proper round
1865 $params[2] = $r; # round mode by caller or undef
1866 $fallback = 1; # to clear a/p afterwards
1870 # the 4 below is empirical, and there might be cases where it is not
1872 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1875 # when user set globals, they would interfere with our calculation, so
1876 # disable them and later re-enable them
1878 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1879 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1880 # we also need to disable any set A or P on $x (_find_round_parameters took
1881 # them already into account), since these would interfere, too
1882 delete $x->{_a}; delete $x->{_p};
1883 # need to disable $upgrade in BigInt, to avoid deep recursion
1884 local $Math::BigInt::upgrade = undef;
1886 my ($limit,$v,$u,$below,$factor,$next,$over);
1888 $u = $x->copy()->blog(undef,$scale)->bmul($y);
1889 $v = $self->bone(); # 1
1890 $factor = $self->new(2); # 2
1891 $x->bone(); # first term: 1
1893 $below = $v->copy();
1896 $limit = $self->new("1E-". ($scale-1));
1900 # we calculate the next term, and add it to the last
1901 # when the next term is below our limit, it won't affect the outcome
1902 # anymore, so we stop
1903 $next = $over->copy()->bdiv($below,$scale);
1904 last if $next->bacmp($limit) <= 0;
1906 # calculate things for the next term
1907 $over *= $u; $below *= $factor; $factor->binc();
1909 last if $x->{sign} !~ /^[-+]$/;
1914 # shortcut to not run through _find_round_parameters again
1915 if (defined $params[0])
1917 $x->bround($params[0],$params[2]); # then round accordingly
1921 $x->bfround($params[1],$params[2]); # then round accordingly
1925 # clear a/p after round, since user did not request it
1926 delete $x->{_a}; delete $x->{_p};
1929 $$abr = $ab; $$pbr = $pb;
1935 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1936 # compute power of two numbers, second arg is used as integer
1937 # modifies first argument
1940 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1941 # objectify is costly, so avoid it
1942 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1944 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1947 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1948 return $x if $x->{sign} =~ /^[+-]inf$/;
1951 return $x->bnan() if $x->{sign} eq '-' && $y->{sign} eq '-';
1953 # cache the result of is_zero
1954 my $y_is_zero = $y->is_zero();
1955 return $x->bone() if $y_is_zero;
1956 return $x if $x->is_one() || $y->is_one();
1958 my $x_is_zero = $x->is_zero();
1959 return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int(); # non-integer power
1961 my $y1 = $y->as_number()->{value}; # make MBI part
1964 if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
1966 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
1967 return $MBI->_is_odd($y1) ? $x : $x->babs(1);
1971 return $x->bone() if $y_is_zero;
1972 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1973 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
1978 $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
1980 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1981 $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
1982 $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
1984 $x->{sign} = $new_sign;
1986 if ($y->{sign} eq '-')
1988 # modify $x in place!
1989 my $z = $x->copy(); $x->bone();
1990 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
1992 $x->round($a,$p,$r,$y);
1995 ###############################################################################
1996 # rounding functions
2000 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
2001 # $n == 0 means round to integer
2002 # expects and returns normalized numbers!
2003 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2005 my ($scale,$mode) = $x->_scale_p(@_);
2006 return $x if !defined $scale || $x->modify('bfround'); # no-op
2008 # never round a 0, +-inf, NaN
2011 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
2014 return $x if $x->{sign} !~ /^[+-]$/;
2016 # don't round if x already has lower precision
2017 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
2019 $x->{_p} = $scale; # remember round in any case
2020 delete $x->{_a}; # and clear A
2023 # round right from the '.'
2025 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round
2027 $scale = -$scale; # positive for simplicity
2028 my $len = $MBI->_len($x->{_m}); # length of mantissa
2030 # the following poses a restriction on _e, but if _e is bigger than a
2031 # scalar, you got other problems (memory etc) anyway
2032 my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e}))); # digits after dot
2033 my $zad = 0; # zeros after dot
2034 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
2036 # p rint "scale $scale dad $dad zad $zad len $len\n";
2037 # number bsstr len zad dad
2038 # 0.123 123e-3 3 0 3
2039 # 0.0123 123e-4 3 1 4
2042 # 1.2345 12345e-4 5 0 4
2044 # do not round after/right of the $dad
2045 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
2047 # round to zero if rounding inside the $zad, but not for last zero like:
2048 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
2049 return $x->bzero() if $scale < $zad;
2050 if ($scale == $zad) # for 0.006, scale -3 and trunc
2056 # adjust round-point to be inside mantissa
2059 $scale = $scale-$zad;
2063 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
2064 $scale = $dbd+$scale;
2070 # round left from the '.'
2072 # 123 => 100 means length(123) = 3 - $scale (2) => 1
2074 my $dbt = $MBI->_len($x->{_m});
2076 my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
2077 # should be the same, so treat it as this
2078 $scale = 1 if $scale == 0;
2079 # shortcut if already integer
2080 return $x if $scale == 1 && $dbt <= $dbd;
2081 # maximum digits before dot
2086 # not enough digits before dot, so round to zero
2089 elsif ( $scale == $dbd )
2096 $scale = $dbd - $scale;
2099 # pass sign to bround for rounding modes '+inf' and '-inf'
2100 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2101 $m->bround($scale,$mode);
2102 $x->{_m} = $m->{value}; # get our mantissa back
2108 # accuracy: preserve $N digits, and overwrite the rest with 0's
2109 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2111 if (($_[0] || 0) < 0)
2113 require Carp; Carp::croak ('bround() needs positive accuracy');
2116 my ($scale,$mode) = $x->_scale_a(@_);
2117 return $x if !defined $scale || $x->modify('bround'); # no-op
2119 # scale is now either $x->{_a}, $accuracy, or the user parameter
2120 # test whether $x already has lower accuracy, do nothing in this case
2121 # but do round if the accuracy is the same, since a math operation might
2122 # want to round a number with A=5 to 5 digits afterwards again
2123 return $x if defined $x->{_a} && $x->{_a} < $scale;
2125 # scale < 0 makes no sense
2126 # scale == 0 => keep all digits
2127 # never round a +-inf, NaN
2128 return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
2130 # 1: never round a 0
2131 # 2: if we should keep more digits than the mantissa has, do nothing
2132 if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
2134 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
2138 # pass sign to bround for '+inf' and '-inf' rounding modes
2139 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2141 $m->bround($scale,$mode); # round mantissa
2142 $x->{_m} = $m->{value}; # get our mantissa back
2143 $x->{_a} = $scale; # remember rounding
2144 delete $x->{_p}; # and clear P
2145 $x->bnorm(); # del trailing zeros gen. by bround()
2150 # return integer less or equal then $x
2151 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2153 return $x if $x->modify('bfloor');
2155 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2157 # if $x has digits after dot
2158 if ($x->{_es} eq '-')
2160 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2161 $x->{_e} = $MBI->_zero(); # trunc/norm
2162 $x->{_es} = '+'; # abs e
2163 $MBI->_inc($x->{_m}) if $x->{sign} eq '-'; # increment if negative
2165 $x->round($a,$p,$r);
2170 # return integer greater or equal then $x
2171 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2173 return $x if $x->modify('bceil');
2174 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2176 # if $x has digits after dot
2177 if ($x->{_es} eq '-')
2179 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2180 $x->{_e} = $MBI->_zero(); # trunc/norm
2181 $x->{_es} = '+'; # abs e
2182 $MBI->_inc($x->{_m}) if $x->{sign} eq '+'; # increment if positive
2184 $x->round($a,$p,$r);
2189 # shift right by $y (divide by power of $n)
2192 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2193 # objectify is costly, so avoid it
2194 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2196 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2199 return $x if $x->modify('brsft');
2200 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2202 $n = 2 if !defined $n; $n = $self->new($n);
2203 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
2208 # shift left by $y (multiply by power of $n)
2211 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2212 # objectify is costly, so avoid it
2213 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2215 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2218 return $x if $x->modify('blsft');
2219 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2221 $n = 2 if !defined $n; $n = $self->new($n);
2222 $x->bmul($n->bpow($y),$a,$p,$r,$y);
2225 ###############################################################################
2229 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
2234 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
2235 # or falling back to MBI::bxxx()
2236 my $name = $AUTOLOAD;
2238 $name =~ s/(.*):://; # split package
2239 my $c = $1 || $class;
2241 $c->import() if $IMPORT == 0;
2242 if (!method_alias($name))
2246 # delayed load of Carp and avoid recursion
2248 Carp::croak ("$c: Can't call a method without name");
2250 if (!method_hand_up($name))
2252 # delayed load of Carp and avoid recursion
2254 Carp::croak ("Can't call $c\-\>$name, not a valid method");
2256 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
2258 return &{"Math::BigInt"."::$name"}(@_);
2260 my $bname = $name; $bname =~ s/^f/b/;
2268 # return a copy of the exponent
2269 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2271 if ($x->{sign} !~ /^[+-]$/)
2273 my $s = $x->{sign}; $s =~ s/^[+-]//;
2274 return Math::BigInt->new($s); # -inf, +inf => +inf
2276 Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
2281 # return a copy of the mantissa
2282 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2284 if ($x->{sign} !~ /^[+-]$/)
2286 my $s = $x->{sign}; $s =~ s/^[+]//;
2287 return Math::BigInt->new($s); # -inf, +inf => +inf
2289 my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
2290 $m->bneg() if $x->{sign} eq '-';
2297 # return a copy of both the exponent and the mantissa
2298 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2300 if ($x->{sign} !~ /^[+-]$/)
2302 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2303 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2305 my $m = Math::BigInt->bzero();
2306 $m->{value} = $MBI->_copy($x->{_m});
2307 $m->bneg() if $x->{sign} eq '-';
2308 ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
2311 ##############################################################################
2312 # private stuff (internal use only)
2318 my $lib = ''; my @a;
2320 for ( my $i = 0; $i < $l ; $i++)
2322 if ( $_[$i] eq ':constant' )
2324 # This causes overlord er load to step in. 'binary' and 'integer'
2325 # are handled by BigInt.
2326 overload::constant float => sub { $self->new(shift); };
2328 elsif ($_[$i] eq 'upgrade')
2330 # this causes upgrading
2331 $upgrade = $_[$i+1]; # or undef to disable
2334 elsif ($_[$i] eq 'downgrade')
2336 # this causes downgrading
2337 $downgrade = $_[$i+1]; # or undef to disable
2340 elsif ($_[$i] eq 'lib')
2342 # alternative library
2343 $lib = $_[$i+1] || ''; # default Calc
2346 elsif ($_[$i] eq 'with')
2348 # alternative class for our private parts()
2349 # XXX: no longer supported
2350 # $MBI = $_[$i+1] || 'Math::BigInt';
2359 $lib =~ tr/a-zA-Z0-9,://cd; # restrict to sane characters
2360 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2361 my $mbilib = eval { Math::BigInt->config()->{lib} };
2362 if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
2364 # MBI already loaded
2365 Math::BigInt->import('lib',"$lib,$mbilib", 'objectify');
2369 # MBI not loaded, or with ne "Math::BigInt::Calc"
2370 $lib .= ",$mbilib" if defined $mbilib;
2371 $lib =~ s/^,//; # don't leave empty
2373 # replacement library can handle lib statement, but also could ignore it
2375 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2376 # used in the same script, or eval inside import(). So we require MBI:
2377 require Math::BigInt;
2378 Math::BigInt->import( lib => $lib, 'objectify' );
2382 require Carp; Carp::croak ("Couldn't load $lib: $! $@");
2384 # find out which one was actually loaded
2385 $MBI = Math::BigInt->config()->{lib};
2387 # register us with MBI to get notified of future lib changes
2388 Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
2390 # any non :constant stuff is handled by our parent, Exporter
2391 # even if @_ is empty, to give it a chance
2392 $self->SUPER::import(@a); # for subclasses
2393 $self->export_to_level(1,$self,@a); # need this, too
2398 # adjust m and e so that m is smallest possible
2399 # round number according to accuracy and precision settings
2400 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
2402 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2404 my $zeros = $MBI->_zeros($x->{_m}); # correct for trailing zeros
2407 my $z = $MBI->_new($zeros);
2408 $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
2409 if ($x->{_es} eq '-')
2411 if ($MBI->_acmp($x->{_e},$z) >= 0)
2413 $x->{_e} = $MBI->_sub ($x->{_e}, $z);
2414 $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
2418 $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
2424 $x->{_e} = $MBI->_add ($x->{_e}, $z);
2429 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
2430 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
2431 $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
2432 if $MBI->_is_zero($x->{_m});
2435 $x; # MBI bnorm is no-op, so dont call it
2438 ##############################################################################
2442 # return number as hexadecimal string (only for integers defined)
2443 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2445 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2446 return '0x0' if $x->is_zero();
2448 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2450 my $z = $MBI->_copy($x->{_m});
2451 if (! $MBI->_is_zero($x->{_e})) # > 0
2453 $MBI->_lsft( $z, $x->{_e},10);
2455 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2461 # return number as binary digit string (only for integers defined)
2462 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2464 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2465 return '0b0' if $x->is_zero();
2467 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2469 my $z = $MBI->_copy($x->{_m});
2470 if (! $MBI->_is_zero($x->{_e})) # > 0
2472 $MBI->_lsft( $z, $x->{_e},10);
2474 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2480 # return copy as a bigint representation of this BigFloat number
2481 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2483 my $z = $MBI->_copy($x->{_m});
2484 if ($x->{_es} eq '-') # < 0
2486 $MBI->_rsft( $z, $x->{_e},10);
2488 elsif (! $MBI->_is_zero($x->{_e})) # > 0
2490 $MBI->_lsft( $z, $x->{_e},10);
2492 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2499 my $class = ref($x) || $x;
2500 $x = $class->new(shift) unless ref($x);
2502 return 1 if $MBI->_is_zero($x->{_m});
2504 my $len = $MBI->_len($x->{_m});
2505 $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
2509 $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
2520 Math::BigFloat - Arbitrary size floating point math package
2527 $x = Math::BigFloat->new($str); # defaults to 0
2528 $nan = Math::BigFloat->bnan(); # create a NotANumber
2529 $zero = Math::BigFloat->bzero(); # create a +0
2530 $inf = Math::BigFloat->binf(); # create a +inf
2531 $inf = Math::BigFloat->binf('-'); # create a -inf
2532 $one = Math::BigFloat->bone(); # create a +1
2533 $one = Math::BigFloat->bone('-'); # create a -1
2536 $x->is_zero(); # true if arg is +0
2537 $x->is_nan(); # true if arg is NaN
2538 $x->is_one(); # true if arg is +1
2539 $x->is_one('-'); # true if arg is -1
2540 $x->is_odd(); # true if odd, false for even
2541 $x->is_even(); # true if even, false for odd
2542 $x->is_pos(); # true if >= 0
2543 $x->is_neg(); # true if < 0
2544 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2546 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2547 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2548 $x->sign(); # return the sign, either +,- or NaN
2549 $x->digit($n); # return the nth digit, counting from right
2550 $x->digit(-$n); # return the nth digit, counting from left
2552 # The following all modify their first argument. If you want to preserve
2553 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2554 # neccessary when mixing $a = $b assigments with non-overloaded math.
2557 $x->bzero(); # set $i to 0
2558 $x->bnan(); # set $i to NaN
2559 $x->bone(); # set $x to +1
2560 $x->bone('-'); # set $x to -1
2561 $x->binf(); # set $x to inf
2562 $x->binf('-'); # set $x to -inf
2564 $x->bneg(); # negation
2565 $x->babs(); # absolute value
2566 $x->bnorm(); # normalize (no-op)
2567 $x->bnot(); # two's complement (bit wise not)
2568 $x->binc(); # increment x by 1
2569 $x->bdec(); # decrement x by 1
2571 $x->badd($y); # addition (add $y to $x)
2572 $x->bsub($y); # subtraction (subtract $y from $x)
2573 $x->bmul($y); # multiplication (multiply $x by $y)
2574 $x->bdiv($y); # divide, set $x to quotient
2575 # return (quo,rem) or quo if scalar
2577 $x->bmod($y); # modulus ($x % $y)
2578 $x->bpow($y); # power of arguments ($x ** $y)
2579 $x->blsft($y); # left shift
2580 $x->brsft($y); # right shift
2581 # return (quo,rem) or quo if scalar
2583 $x->blog(); # logarithm of $x to base e (Euler's number)
2584 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
2586 $x->band($y); # bit-wise and
2587 $x->bior($y); # bit-wise inclusive or
2588 $x->bxor($y); # bit-wise exclusive or
2589 $x->bnot(); # bit-wise not (two's complement)
2591 $x->bsqrt(); # calculate square-root
2592 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
2593 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2595 $x->bround($N); # accuracy: preserve $N digits
2596 $x->bfround($N); # precision: round to the $Nth digit
2598 $x->bfloor(); # return integer less or equal than $x
2599 $x->bceil(); # return integer greater or equal than $x
2601 # The following do not modify their arguments:
2603 bgcd(@values); # greatest common divisor
2604 blcm(@values); # lowest common multiplicator
2606 $x->bstr(); # return string
2607 $x->bsstr(); # return string in scientific notation
2609 $x->as_int(); # return $x as BigInt
2610 $x->exponent(); # return exponent as BigInt
2611 $x->mantissa(); # return mantissa as BigInt
2612 $x->parts(); # return (mantissa,exponent) as BigInt
2614 $x->length(); # number of digits (w/o sign and '.')
2615 ($l,$f) = $x->length(); # number of digits, and length of fraction
2617 $x->precision(); # return P of $x (or global, if P of $x undef)
2618 $x->precision($n); # set P of $x to $n
2619 $x->accuracy(); # return A of $x (or global, if A of $x undef)
2620 $x->accuracy($n); # set A $x to $n
2622 # these get/set the appropriate global value for all BigFloat objects
2623 Math::BigFloat->precision(); # Precision
2624 Math::BigFloat->accuracy(); # Accuracy
2625 Math::BigFloat->round_mode(); # rounding mode
2629 All operators (inlcuding basic math operations) are overloaded if you
2630 declare your big floating point numbers as
2632 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2634 Operations with overloaded operators preserve the arguments, which is
2635 exactly what you expect.
2637 =head2 Canonical notation
2639 Input to these routines are either BigFloat objects, or strings of the
2640 following four forms:
2654 C</^[+-]\d+E[+-]?\d+$/>
2658 C</^[+-]\d*\.\d+E[+-]?\d+$/>
2662 all with optional leading and trailing zeros and/or spaces. Additonally,
2663 numbers are allowed to have an underscore between any two digits.
2665 Empty strings as well as other illegal numbers results in 'NaN'.
2667 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2668 are always stored in normalized form. On a string, it creates a BigFloat
2673 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2675 The string output will always have leading and trailing zeros stripped and drop
2676 a plus sign. C<bstr()> will give you always the form with a decimal point,
2677 while C<bsstr()> (s for scientific) gives you the scientific notation.
2679 Input bstr() bsstr()
2681 ' -123 123 123' '-123123123' '-123123123E0'
2682 '00.0123' '0.0123' '123E-4'
2683 '123.45E-2' '1.2345' '12345E-4'
2684 '10E+3' '10000' '1E4'
2686 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2687 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2688 return either undef, <0, 0 or >0 and are suited for sort.
2690 Actual math is done by using the class defined with C<with => Class;> (which
2691 defaults to BigInts) to represent the mantissa and exponent.
2693 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2694 represent the result when input arguments are not numbers, as well as
2695 the result of dividing by zero.
2697 =head2 C<mantissa()>, C<exponent()> and C<parts()>
2699 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2700 as BigInts such that:
2702 $m = $x->mantissa();
2703 $e = $x->exponent();
2704 $y = $m * ( 10 ** $e );
2705 print "ok\n" if $x == $y;
2707 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2709 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2711 Currently the mantissa is reduced as much as possible, favouring higher
2712 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2713 This might change in the future, so do not depend on it.
2715 =head2 Accuracy vs. Precision
2717 See also: L<Rounding|Rounding>.
2719 Math::BigFloat supports both precision and accuracy. For a full documentation,
2720 examples and tips on these topics please see the large section in
2723 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2724 a operation consumes all resources, each operation produces no more than
2725 the requested number of digits.
2727 Please refer to BigInt's documentation for the precedence rules of which
2728 accuracy/precision setting will be used.
2730 If there is no gloabl precision set, B<and> the operation inquestion was not
2731 called with a requested precision or accuracy, B<and> the input $x has no
2732 accuracy or precision set, then a fallback parameter will be used. For
2733 historical reasons, it is called C<div_scale> and can be accessed via:
2735 $d = Math::BigFloat->div_scale(); # query
2736 Math::BigFloat->div_scale($n); # set to $n digits
2738 The default value is 40 digits.
2740 In case the result of one operation has more precision than specified,
2741 it is rounded. The rounding mode taken is either the default mode, or the one
2742 supplied to the operation after the I<scale>:
2744 $x = Math::BigFloat->new(2);
2745 Math::BigFloat->accuracy(5); # 5 digits max
2746 $y = $x->copy()->bdiv(3); # will give 0.66667
2747 $y = $x->copy()->bdiv(3,6); # will give 0.666667
2748 $y = $x->copy()->bdiv(3,6,undef,'odd'); # will give 0.666667
2749 Math::BigFloat->round_mode('zero');
2750 $y = $x->copy()->bdiv(3,6); # will also give 0.666667
2752 Note that C<< Math::BigFloat->accuracy() >> and C<< Math::BigFloat->precision() >>
2753 set the global variables, and thus B<any> newly created number will be subject
2754 to the global rounding. This means that in the examples above, the C<3>
2755 as argument to C<bdiv()> will also get an accuracy of B<5>.
2757 It is less confusing to either calculate the result fully, and afterwards
2758 round it explicitely, or use the additional parameters to the math
2762 $x = Math::BigFloat->new(2);
2763 $y = $x->copy()->bdiv(3);
2764 print $y->bround(5),"\n"; # will give 0.66667
2769 $x = Math::BigFloat->new(2);
2770 $y = $x->copy()->bdiv(3,5); # will give 0.66667
2777 =item ffround ( +$scale )
2779 Rounds to the $scale'th place left from the '.', counting from the dot.
2780 The first digit is numbered 1.
2782 =item ffround ( -$scale )
2784 Rounds to the $scale'th place right from the '.', counting from the dot.
2788 Rounds to an integer.
2790 =item fround ( +$scale )
2792 Preserves accuracy to $scale digits from the left (aka significant digits)
2793 and pads the rest with zeros. If the number is between 1 and -1, the
2794 significant digits count from the first non-zero after the '.'
2796 =item fround ( -$scale ) and fround ( 0 )
2798 These are effectively no-ops.
2802 All rounding functions take as a second parameter a rounding mode from one of
2803 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2805 The default rounding mode is 'even'. By using
2806 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
2807 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2808 no longer supported.
2809 The second parameter to the round functions then overrides the default
2812 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
2813 'trunc' as rounding mode to make it equivalent to:
2818 You can override this by passing the desired rounding mode as parameter to
2821 $x = Math::BigFloat->new(2.5);
2822 $y = $x->as_number('odd'); # $y = 3
2828 =head1 Autocreating constants
2830 After C<use Math::BigFloat ':constant'> all the floating point constants
2831 in the given scope are converted to C<Math::BigFloat>. This conversion
2832 happens at compile time.
2836 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2838 prints the value of C<2E-100>. Note that without conversion of
2839 constants the expression 2E-100 will be calculated as normal floating point
2842 Please note that ':constant' does not affect integer constants, nor binary
2843 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2848 Math with the numbers is done (by default) by a module called
2849 Math::BigInt::Calc. This is equivalent to saying:
2851 use Math::BigFloat lib => 'Calc';
2853 You can change this by using:
2855 use Math::BigFloat lib => 'BitVect';
2857 The following would first try to find Math::BigInt::Foo, then
2858 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2860 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2862 Calc.pm uses as internal format an array of elements of some decimal base
2863 (usually 1e7, but this might be differen for some systems) with the least
2864 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2865 significant bit first. Other modules might use even different means of
2866 representing the numbers. See the respective module documentation for further
2869 Please note that Math::BigFloat does B<not> use the denoted library itself,
2870 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2873 use Math::BigInt lib => 'GMP';
2876 you can roll it all into one line:
2878 use Math::BigFloat lib => 'GMP';
2880 It is also possible to just require Math::BigFloat:
2882 require Math::BigFloat;
2884 This will load the neccessary things (like BigInt) when they are needed, and
2887 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2888 you ever wanted to know about loading a different library.
2890 =head2 Using Math::BigInt::Lite
2892 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2895 use Math::BigFloat with => 'Math::BigInt::Lite';
2897 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2898 can combine these if you want. For instance, you may want to use
2899 Math::BigInt objects in your main script, too.
2903 use Math::BigFloat with => 'Math::BigInt::Lite';
2905 Of course, you can combine this with the C<lib> parameter.
2908 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2910 There is no need for a "use Math::BigInt;" statement, even if you want to
2911 use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2912 always loads it. But if you add it, add it B<before>:
2916 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2918 Notice that the module with the last C<lib> will "win" and thus
2919 it's lib will be used if the lib is available:
2922 use Math::BigInt lib => 'Bar,Baz';
2923 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2925 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2926 words, Math::BigFloat will try to retain previously loaded libs when you
2927 don't specify it onem but if you specify one, it will try to load them.
2929 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2930 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2931 same as trying the latter load alone, except for the fact that one of Bar or
2932 Baz might be loaded needlessly in an intermidiate step (and thus hang around
2933 and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2934 will still be tried to be loaded, but this is not as time/memory consuming as
2935 actually loading one of them. Still, this type of usage is not recommended due
2938 The old way (loading the lib only in BigInt) still works though:
2941 use Math::BigInt lib => 'Bar,Baz';
2944 You can even load Math::BigInt afterwards:
2948 use Math::BigInt lib => 'Bar,Baz';
2950 But this has the same problems like #5, it will first load Calc
2951 (Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2952 Baz, depending on which of them works and is usable/loadable. Since this
2953 loads Calc unnecc., it is not recommended.
2955 Since it also possible to just require Math::BigFloat, this poses the question
2956 about what libary this will use:
2958 require Math::BigFloat;
2959 my $x = Math::BigFloat->new(123); $x += 123;
2961 It will use Calc. Please note that the call to import() is still done, but
2962 only when you use for the first time some Math::BigFloat math (it is triggered
2963 via any constructor, so the first time you create a Math::BigFloat, the load
2964 will happen in the background). This means:
2966 require Math::BigFloat;
2967 Math::BigFloat->import ( lib => 'Foo,Bar' );
2969 would be the same as:
2971 use Math::BigFloat lib => 'Foo, Bar';
2973 But don't try to be clever to insert some operations in between:
2975 require Math::BigFloat;
2976 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
2977 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
2978 $x = Math::BigFloat->bone()+4; # now use Pari
2980 While this works, it loads Calc needlessly. But maybe you just wanted that?
2982 B<Examples #3 is highly recommended> for daily usage.
2986 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
2992 =item stringify, bstr()
2994 Both stringify and bstr() now drop the leading '+'. The old code would return
2995 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2996 reasoning and details.
3000 The following will probably not do what you expect:
3002 print $c->bdiv(123.456),"\n";
3004 It prints both quotient and reminder since print works in list context. Also,
3005 bdiv() will modify $c, so be carefull. You probably want to use
3007 print $c / 123.456,"\n";
3008 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
3012 =item Modifying and =
3016 $x = Math::BigFloat->new(5);
3019 It will not do what you think, e.g. making a copy of $x. Instead it just makes
3020 a second reference to the B<same> object and stores it in $y. Thus anything
3021 that modifies $x will modify $y (except overloaded math operators), and vice
3022 versa. See L<Math::BigInt> for details and how to avoid that.
3026 C<bpow()> now modifies the first argument, unlike the old code which left
3027 it alone and only returned the result. This is to be consistent with
3028 C<badd()> etc. The first will modify $x, the second one won't:
3030 print bpow($x,$i),"\n"; # modify $x
3031 print $x->bpow($i),"\n"; # ditto
3032 print $x ** $i,"\n"; # leave $x alone
3038 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
3039 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
3041 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
3042 because they solve the autoupgrading/downgrading issue, at least partly.
3045 L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
3046 more documentation including a full version history, testcases, empty
3047 subclass files and benchmarks.
3051 This program is free software; you may redistribute it and/or modify it under
3052 the same terms as Perl itself.
3056 Mark Biggar, overloaded interface by Ilya Zakharevich.
3057 Completely rewritten by Tels L<http://bloodgate.com> in 2001 - 2004, and still