1 package Math::BigFloat;
4 # Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
7 # The following hash values are internally used:
8 # _e : exponent (ref to $CALC object)
9 # _m : mantissa (ref to $CALC object)
11 # sign : +,-,+inf,-inf, or "NaN" if not a number
19 @ISA = qw(Exporter Math::BigInt);
22 # $_trap_inf/$_trap_nan are internal and should never be accessed from outside
23 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode
24 $upgrade $downgrade $_trap_nan $_trap_inf/;
25 my $class = "Math::BigFloat";
28 '<=>' => sub { my $rc = $_[2] ?
29 ref($_[0])->bcmp($_[1],$_[0]) :
30 ref($_[0])->bcmp($_[0],$_[1]);
31 $rc = 1 unless defined $rc;
34 # we need '>=' to get things like "1 >= NaN" right:
35 '>=' => sub { my $rc = $_[2] ?
36 ref($_[0])->bcmp($_[1],$_[0]) :
37 ref($_[0])->bcmp($_[0],$_[1]);
38 # if there was a NaN involved, return false
39 return '' unless defined $rc;
42 'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
45 ##############################################################################
46 # global constants, flags and assorted stuff
48 # the following are public, but their usage is not recommended. Use the
49 # accessor methods instead.
51 # class constants, use Class->constant_name() to access
52 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
59 # the package we are using for our private parts, defaults to:
60 # Math::BigInt->config()->{lib}
61 my $MBI = 'Math::BigInt::FastCalc';
63 # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
65 # the same for infinity
68 # constant for easier life
71 my $IMPORT = 0; # was import() called yet? used to make require work
73 # some digits of accuracy for blog(undef,10); which we use in blog() for speed
75 '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
76 my $LOG_10_A = length($LOG_10)-1;
79 '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
80 my $LOG_2_A = length($LOG_2)-1;
81 my $HALF = '0.5'; # made into an object if nec.
83 ##############################################################################
84 # the old code had $rnd_mode, so we need to support it, too
86 sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
87 sub FETCH { return $round_mode; }
88 sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
92 # when someone sets $rnd_mode, we catch this and check the value to see
93 # whether it is valid or not.
94 $rnd_mode = 'even'; tie $rnd_mode, 'Math::BigFloat';
96 # we need both of them in this package:
97 *as_int = \&as_number;
100 ##############################################################################
103 # valid method aliases for AUTOLOAD
104 my %methods = map { $_ => 1 }
105 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
106 fint facmp fcmp fzero fnan finf finc fdec ffac fneg
107 fceil ffloor frsft flsft fone flog froot fexp
109 # valid methods that can be handed up (for AUTOLOAD)
110 my %hand_ups = map { $_ => 1 }
111 qw / is_nan is_inf is_negative is_positive is_pos is_neg
112 accuracy precision div_scale round_mode fabs fnot
113 objectify upgrade downgrade
118 sub _method_alias { exists $methods{$_[0]||''}; }
119 sub _method_hand_up { exists $hand_ups{$_[0]||''}; }
122 ##############################################################################
127 # create a new BigFloat object from a string or another bigfloat object.
130 # sign => sign (+/-), or "NaN"
132 my ($class,$wanted,@r) = @_;
134 # avoid numify-calls by not using || on $wanted!
135 return $class->bzero() if !defined $wanted; # default to 0
136 return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
138 $class->import() if $IMPORT == 0; # make require work
140 my $self = {}; bless $self, $class;
141 # shortcut for bigints and its subclasses
142 if ((ref($wanted)) && UNIVERSAL::can( $wanted, "as_number"))
144 $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy
145 $self->{_e} = $MBI->_zero();
147 $self->{sign} = $wanted->sign();
148 return $self->bnorm();
150 # else: got a string or something maskerading as number (with overload)
152 # handle '+inf', '-inf' first
153 if ($wanted =~ /^[+-]?inf\z/)
155 return $downgrade->new($wanted) if $downgrade;
157 $self->{sign} = $wanted; # set a default sign for bstr()
158 return $self->binf($wanted);
161 # shortcut for simple forms like '12' that neither have trailing nor leading
163 if ($wanted =~ /^([+-]?)([1-9][0-9]*[1-9])$/)
165 $self->{_e} = $MBI->_zero();
167 $self->{sign} = $1 || '+';
168 $self->{_m} = $MBI->_new($2);
169 return $self->round(@r) if !$downgrade;
172 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
178 Carp::croak ("$wanted is not a number initialized to $class");
181 return $downgrade->bnan() if $downgrade;
183 $self->{_e} = $MBI->_zero();
185 $self->{_m} = $MBI->_zero();
186 $self->{sign} = $nan;
190 # make integer from mantissa by adjusting exp, then convert to int
191 $self->{_e} = $MBI->_new($$ev); # exponent
192 $self->{_es} = $$es || '+';
193 my $mantissa = "$$miv$$mfv"; # create mant.
194 $mantissa =~ s/^0+(\d)/$1/; # strip leading zeros
195 $self->{_m} = $MBI->_new($mantissa); # create mant.
197 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
198 if (CORE::length($$mfv) != 0)
200 my $len = $MBI->_new( CORE::length($$mfv));
201 ($self->{_e}, $self->{_es}) =
202 _e_sub ($self->{_e}, $len, $self->{_es}, '+');
204 # we can only have trailing zeros on the mantissa if $$mfv eq ''
207 # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
208 # because that is faster, especially when _m is not stored in base 10.
209 my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/;
212 my $z = $MBI->_new($zeros);
213 # turn '120e2' into '12e3'
214 $MBI->_rsft ( $self->{_m}, $z, 10);
215 ($self->{_e}, $self->{_es}) =
216 _e_add ( $self->{_e}, $z, $self->{_es}, '+');
219 $self->{sign} = $$mis;
221 # for something like 0Ey, set y to 1, and -0 => +0
222 # Check $$miv for being '0' and $$mfv eq '', because otherwise _m could not
223 # have become 0. That's faster than to call $MBI->_is_zero().
224 $self->{sign} = '+', $self->{_e} = $MBI->_one()
225 if $$miv eq '0' and $$mfv eq '';
227 return $self->round(@r) if !$downgrade;
229 # if downgrade, inf, NaN or integers go down
231 if ($downgrade && $self->{_es} eq '+')
233 if ($MBI->_is_zero( $self->{_e} ))
235 return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
237 return $downgrade->new($self->bsstr());
239 $self->bnorm()->round(@r); # first normalize, then round
247 # if two arguments, the first one is the class to "swallow" subclasses
255 return unless ref($x); # only for objects
257 my $self = {}; bless $self,$c;
259 $self->{sign} = $x->{sign};
260 $self->{_es} = $x->{_es};
261 $self->{_m} = $MBI->_copy($x->{_m});
262 $self->{_e} = $MBI->_copy($x->{_e});
263 $self->{_a} = $x->{_a} if defined $x->{_a};
264 $self->{_p} = $x->{_p} if defined $x->{_p};
270 # used by parent class bone() to initialize number to NaN
276 my $class = ref($self);
277 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
280 $IMPORT=1; # call our import only once
281 $self->{_m} = $MBI->_zero();
282 $self->{_e} = $MBI->_zero();
288 # used by parent class bone() to initialize number to +-inf
294 my $class = ref($self);
295 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
298 $IMPORT=1; # call our import only once
299 $self->{_m} = $MBI->_zero();
300 $self->{_e} = $MBI->_zero();
306 # used by parent class bone() to initialize number to 1
308 $IMPORT=1; # call our import only once
309 $self->{_m} = $MBI->_one();
310 $self->{_e} = $MBI->_zero();
316 # used by parent class bone() to initialize number to 0
318 $IMPORT=1; # call our import only once
319 $self->{_m} = $MBI->_zero();
320 $self->{_e} = $MBI->_one();
326 my ($self,$class) = @_;
327 return if $class =~ /^Math::BigInt/; # we aren't one of these
328 UNIVERSAL::isa($self,$class);
333 # return (later set?) configuration data as hash ref
334 my $class = shift || 'Math::BigFloat';
336 if (@_ == 1 && ref($_[0]) ne 'HASH')
338 my $cfg = $class->SUPER::config();
339 return $cfg->{$_[0]};
342 my $cfg = $class->SUPER::config(@_);
344 # now we need only to override the ones that are different from our parent
345 $cfg->{class} = $class;
350 ##############################################################################
351 # string conversation
355 # (ref to BFLOAT or num_str ) return num_str
356 # Convert number from internal format to (non-scientific) string format.
357 # internal format is always normalized (no leading zeros, "-0" => "+0")
358 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
360 if ($x->{sign} !~ /^[+-]$/)
362 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
366 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
369 my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
372 $es = $MBI->_str($x->{_m});
373 $len = CORE::length($es);
374 my $e = $MBI->_num($x->{_e});
375 $e = -$e if $x->{_es} eq '-';
379 # if _e is bigger than a scalar, the following will blow your memory
382 my $r = abs($e) - $len;
383 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
387 substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
388 $cad = -$cad if $x->{_es} eq '-';
394 $es .= '0' x $e; $len += $e; $cad = 0;
398 $es = '-'.$es if $x->{sign} eq '-';
399 # if set accuracy or precision, pad with zeros on the right side
400 if ((defined $x->{_a}) && ($not_zero))
402 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
403 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
404 $zeros = $x->{_a} - $len if $cad != $len;
405 $es .= $dot.'0' x $zeros if $zeros > 0;
407 elsif ((($x->{_p} || 0) < 0))
409 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
410 my $zeros = -$x->{_p} + $cad;
411 $es .= $dot.'0' x $zeros if $zeros > 0;
418 # (ref to BFLOAT or num_str ) return num_str
419 # Convert number from internal format to scientific string format.
420 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
421 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
423 if ($x->{sign} !~ /^[+-]$/)
425 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
428 my $sep = 'e'.$x->{_es};
429 my $sign = $x->{sign}; $sign = '' if $sign eq '+';
430 $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
435 # Make a number from a BigFloat object
436 # simple return a string and let Perl's atoi()/atof() handle the rest
437 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
441 ##############################################################################
442 # public stuff (usually prefixed with "b")
446 # (BINT or num_str) return BINT
447 # negate number or make a negated number from string
448 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
450 return $x if $x->modify('bneg');
452 # for +0 dont negate (to have always normalized +0). Does nothing for 'NaN'
453 $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
458 # XXX TODO this must be overwritten and return NaN for non-integer values
459 # band(), bior(), bxor(), too
462 # $class->SUPER::bnot($class,@_);
467 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
470 my ($self,$x,$y) = (ref($_[0]),@_);
471 # objectify is costly, so avoid it
472 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
474 ($self,$x,$y) = objectify(2,@_);
477 return $upgrade->bcmp($x,$y) if defined $upgrade &&
478 ((!$x->isa($self)) || (!$y->isa($self)));
480 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
482 # handle +-inf and NaN
483 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
484 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
485 return +1 if $x->{sign} eq '+inf';
486 return -1 if $x->{sign} eq '-inf';
487 return -1 if $y->{sign} eq '+inf';
491 # check sign for speed first
492 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
493 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
496 my $xz = $x->is_zero();
497 my $yz = $y->is_zero();
498 return 0 if $xz && $yz; # 0 <=> 0
499 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
500 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
502 # adjust so that exponents are equal
503 my $lxm = $MBI->_len($x->{_m});
504 my $lym = $MBI->_len($y->{_m});
505 # the numify somewhat limits our length, but makes it much faster
506 my ($xes,$yes) = (1,1);
507 $xes = -1 if $x->{_es} ne '+';
508 $yes = -1 if $y->{_es} ne '+';
509 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
510 my $ly = $lym + $yes * $MBI->_num($y->{_e});
511 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
512 return $l <=> 0 if $l != 0;
514 # lengths (corrected by exponent) are equal
515 # so make mantissa equal length by padding with zero (shift left)
516 my $diff = $lxm - $lym;
517 my $xm = $x->{_m}; # not yet copy it
521 $ym = $MBI->_copy($y->{_m});
522 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
526 $xm = $MBI->_copy($x->{_m});
527 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
529 my $rc = $MBI->_acmp($xm,$ym);
530 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
536 # Compares 2 values, ignoring their signs.
537 # Returns one of undef, <0, =0, >0. (suitable for sort)
540 my ($self,$x,$y) = (ref($_[0]),@_);
541 # objectify is costly, so avoid it
542 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
544 ($self,$x,$y) = objectify(2,@_);
547 return $upgrade->bacmp($x,$y) if defined $upgrade &&
548 ((!$x->isa($self)) || (!$y->isa($self)));
550 # handle +-inf and NaN's
551 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
553 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
554 return 0 if ($x->is_inf() && $y->is_inf());
555 return 1 if ($x->is_inf() && !$y->is_inf());
560 my $xz = $x->is_zero();
561 my $yz = $y->is_zero();
562 return 0 if $xz && $yz; # 0 <=> 0
563 return -1 if $xz && !$yz; # 0 <=> +y
564 return 1 if $yz && !$xz; # +x <=> 0
566 # adjust so that exponents are equal
567 my $lxm = $MBI->_len($x->{_m});
568 my $lym = $MBI->_len($y->{_m});
569 my ($xes,$yes) = (1,1);
570 $xes = -1 if $x->{_es} ne '+';
571 $yes = -1 if $y->{_es} ne '+';
572 # the numify somewhat limits our length, but makes it much faster
573 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
574 my $ly = $lym + $yes * $MBI->_num($y->{_e});
576 return $l <=> 0 if $l != 0;
578 # lengths (corrected by exponent) are equal
579 # so make mantissa equal-length by padding with zero (shift left)
580 my $diff = $lxm - $lym;
581 my $xm = $x->{_m}; # not yet copy it
585 $ym = $MBI->_copy($y->{_m});
586 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
590 $xm = $MBI->_copy($x->{_m});
591 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
593 $MBI->_acmp($xm,$ym);
598 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
599 # return result as BFLOAT
602 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
603 # objectify is costly, so avoid it
604 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
606 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
609 return $x if $x->modify('badd');
611 # inf and NaN handling
612 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
615 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
617 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
619 # +inf++inf or -inf+-inf => same, rest is NaN
620 return $x if $x->{sign} eq $y->{sign};
623 # +-inf + something => +inf; something +-inf => +-inf
624 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
628 return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
629 ((!$x->isa($self)) || (!$y->isa($self)));
631 # speed: no add for 0+y or x+0
632 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
633 if ($x->is_zero()) # 0+y
635 # make copy, clobbering up x (modify in place!)
636 $x->{_e} = $MBI->_copy($y->{_e});
637 $x->{_es} = $y->{_es};
638 $x->{_m} = $MBI->_copy($y->{_m});
639 $x->{sign} = $y->{sign} || $nan;
640 return $x->round($a,$p,$r,$y);
643 # take lower of the two e's and adapt m1 to it to match m2
645 $e = $MBI->_zero() if !defined $e; # if no BFLOAT?
646 $e = $MBI->_copy($e); # make copy (didn't do it yet)
650 ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
652 my $add = $MBI->_copy($y->{_m});
654 if ($es eq '-') # < 0
656 $MBI->_lsft( $x->{_m}, $e, 10);
657 ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
659 elsif (!$MBI->_is_zero($e)) # > 0
661 $MBI->_lsft($add, $e, 10);
663 # else: both e are the same, so just leave them
665 if ($x->{sign} eq $y->{sign})
668 $x->{_m} = $MBI->_add($x->{_m}, $add);
672 ($x->{_m}, $x->{sign}) =
673 _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
676 # delete trailing zeros, then round
677 $x->bnorm()->round($a,$p,$r,$y);
680 # sub bsub is inherited from Math::BigInt!
684 # increment arg by one
685 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
687 return $x if $x->modify('binc');
689 if ($x->{_es} eq '-')
691 return $x->badd($self->bone(),@r); # digits after dot
694 if (!$MBI->_is_zero($x->{_e})) # _e == 0 for NaN, inf, -inf
696 # 1e2 => 100, so after the shift below _m has a '0' as last digit
697 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
698 $x->{_e} = $MBI->_zero(); # normalize
700 # we know that the last digit of $x will be '1' or '9', depending on the
704 if ($x->{sign} eq '+')
706 $MBI->_inc($x->{_m});
707 return $x->bnorm()->bround(@r);
709 elsif ($x->{sign} eq '-')
711 $MBI->_dec($x->{_m});
712 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
713 return $x->bnorm()->bround(@r);
715 # inf, nan handling etc
716 $x->badd($self->bone(),@r); # badd() does round
721 # decrement arg by one
722 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
724 return $x if $x->modify('bdec');
726 if ($x->{_es} eq '-')
728 return $x->badd($self->bone('-'),@r); # digits after dot
731 if (!$MBI->_is_zero($x->{_e}))
733 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
734 $x->{_e} = $MBI->_zero(); # normalize
738 my $zero = $x->is_zero();
740 if (($x->{sign} eq '-') || $zero)
742 $MBI->_inc($x->{_m});
743 $x->{sign} = '-' if $zero; # 0 => 1 => -1
744 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
745 return $x->bnorm()->round(@r);
748 elsif ($x->{sign} eq '+')
750 $MBI->_dec($x->{_m});
751 return $x->bnorm()->round(@r);
753 # inf, nan handling etc
754 $x->badd($self->bone('-'),@r); # does round
761 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
763 return $x if $x->modify('blog');
765 # $base > 0, $base != 1; if $base == undef default to $base == e
768 # we need to limit the accuracy to protect against overflow
771 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
773 # also takes care of the "error in _find_round_parameters?" case
774 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
776 # no rounding at all, so must use fallback
777 if (scalar @params == 0)
779 # simulate old behaviour
780 $params[0] = $self->div_scale(); # and round to it as accuracy
781 $params[1] = undef; # P = undef
782 $scale = $params[0]+4; # at least four more for proper round
783 $params[2] = $r; # round mode by caller or undef
784 $fallback = 1; # to clear a/p afterwards
788 # the 4 below is empirical, and there might be cases where it is not
790 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
793 return $x->bzero(@params) if $x->is_one();
794 # base not defined => base == Euler's number e
797 # make object, since we don't feed it through objectify() to still get the
798 # case of $base == undef
799 $base = $self->new($base) unless ref($base);
800 # $base > 0; $base != 1
801 return $x->bnan() if $base->is_zero() || $base->is_one() ||
802 $base->{sign} ne '+';
803 # if $x == $base, we know the result must be 1.0
804 if ($x->bcmp($base) == 0)
806 $x->bone('+',@params);
809 # clear a/p after round, since user did not request it
810 delete $x->{_a}; delete $x->{_p};
816 # when user set globals, they would interfere with our calculation, so
817 # disable them and later re-enable them
819 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
820 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
821 # we also need to disable any set A or P on $x (_find_round_parameters took
822 # them already into account), since these would interfere, too
823 delete $x->{_a}; delete $x->{_p};
824 # need to disable $upgrade in BigInt, to avoid deep recursion
825 local $Math::BigInt::upgrade = undef;
826 local $Math::BigFloat::downgrade = undef;
828 # upgrade $x if $x is not a BigFloat (handle BigInt input)
830 if (!$x->isa('Math::BigFloat'))
832 $x = Math::BigFloat->new($x);
838 # If the base is defined and an integer, try to calculate integer result
839 # first. This is very fast, and in case the real result was found, we can
841 if (defined $base && $base->is_int() && $x->is_int())
843 my $i = $MBI->_copy( $x->{_m} );
844 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
845 my $int = Math::BigInt->bzero();
847 $int->blog($base->as_number());
849 if ($base->as_number()->bpow($int) == $x)
851 # found result, return it
852 $x->{_m} = $int->{value};
853 $x->{_e} = $MBI->_zero();
862 # base is undef, so base should be e (Euler's number), so first calculate the
863 # log to base e (using reduction by 10 (and probably 2)):
864 $self->_log_10($x,$scale);
866 # and if a different base was requested, convert it
869 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
870 # not ln, but some other base (don't modify $base)
871 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
875 # shortcut to not run through _find_round_parameters again
876 if (defined $params[0])
878 $x->bround($params[0],$params[2]); # then round accordingly
882 $x->bfround($params[1],$params[2]); # then round accordingly
886 # clear a/p after round, since user did not request it
887 delete $x->{_a}; delete $x->{_p};
890 $$abr = $ab; $$pbr = $pb;
897 # Given D (digits in decimal), compute N so that N! (N factorial) is
898 # at least D digits long. D should be at least 50.
901 # two constants for the Ramanujan estimate of ln(N!)
902 my $lg2 = log(2 * 3.14159265) / 2;
905 # D = 50 => N => 42, so L = 40 and R = 50
906 my $l = 40; my $r = $d;
908 # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
909 $l = $l->numify if ref($l);
910 $r = $r->numify if ref($r);
911 $lg2 = $lg2->numify if ref($lg2);
912 $lg10 = $lg10->numify if ref($lg10);
914 # binary search for the right value (could this be written as the reverse of lg(n!)?)
917 my $n = int(($r - $l) / 2) + $l;
919 int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
920 $ramanujan > $d ? $r = $n : $l = $n;
927 # Calculate n over k (binomial coefficient or "choose" function) as integer.
929 my ($self,$x,$y,@r) = (ref($_[0]),@_);
931 # objectify is costly, so avoid it
932 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
934 ($self,$x,$y,@r) = objectify(2,@_);
937 return $x if $x->modify('bnok');
939 return $x->bnan() if $x->is_nan() || $y->is_nan();
940 return $x->binf() if $x->is_inf();
942 my $u = $x->as_int();
943 $u->bnok($y->as_int());
945 $x->{_m} = $u->{value};
946 $x->{_e} = $MBI->_zero();
954 # Calculate e ** X (Euler's number to the power of X)
955 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
957 return $x if $x->modify('bexp');
959 return $x->binf() if $x->{sign} eq '+inf';
960 return $x->bzero() if $x->{sign} eq '-inf';
962 # we need to limit the accuracy to protect against overflow
965 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
967 # also takes care of the "error in _find_round_parameters?" case
968 return $x if $x->{sign} eq 'NaN';
970 # no rounding at all, so must use fallback
971 if (scalar @params == 0)
973 # simulate old behaviour
974 $params[0] = $self->div_scale(); # and round to it as accuracy
975 $params[1] = undef; # P = undef
976 $scale = $params[0]+4; # at least four more for proper round
977 $params[2] = $r; # round mode by caller or undef
978 $fallback = 1; # to clear a/p afterwards
982 # the 4 below is empirical, and there might be cases where it's not enough...
983 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
986 return $x->bone(@params) if $x->is_zero();
988 if (!$x->isa('Math::BigFloat'))
990 $x = Math::BigFloat->new($x);
994 # when user set globals, they would interfere with our calculation, so
995 # disable them and later re-enable them
997 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
998 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
999 # we also need to disable any set A or P on $x (_find_round_parameters took
1000 # them already into account), since these would interfere, too
1001 delete $x->{_a}; delete $x->{_p};
1002 # need to disable $upgrade in BigInt, to avoid deep recursion
1003 local $Math::BigInt::upgrade = undef;
1004 local $Math::BigFloat::downgrade = undef;
1006 my $x_org = $x->copy();
1008 # We use the following Taylor series:
1011 # e = 1 + --- + --- + --- + --- ...
1014 # The difference for each term is X and N, which would result in:
1015 # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term
1017 # But it is faster to compute exp(1) and then raising it to the
1018 # given power, esp. if $x is really big and an integer because:
1020 # * The numerator is always 1, making the computation faster
1021 # * the series converges faster in the case of x == 1
1022 # * We can also easily check when we have reached our limit: when the
1023 # term to be added is smaller than "1E$scale", we can stop - f.i.
1024 # scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
1025 # * we can compute the *exact* result by simulating bigrat math:
1027 # 1 1 gcd(3,4) = 1 1*24 + 1*6 5
1028 # - + - = ---------- = --
1031 # We do not compute the gcd() here, but simple do:
1033 # - + - = --------- = --
1037 # a c a*d + c*b and note that c is always 1 and d = (b*f)
1041 # This leads to: which can be reduced by b to:
1042 # a 1 a*b*f + b a*f + 1
1043 # - + - = --------- = -------
1046 # The first terms in the series are:
1048 # 1 1 1 1 1 1 1 1 13700
1049 # -- + -- + -- + -- + -- + --- + --- + ---- = -----
1050 # 1 1 2 6 24 120 720 5040 5040
1052 # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!
1056 # set $x directly from a cached string form
1057 $x->{_m} = $MBI->_new(
1058 "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
1061 $x->{_e} = $MBI->_new(79);
1065 # compute A and B so that e = A / B.
1067 # After some terms we end up with this, so we use it as a starting point:
1068 my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
1069 my $F = $MBI->_new(42); my $step = 42;
1071 # Compute how many steps we need to take to get $A and $B sufficiently big
1072 my $steps = _len_to_steps($scale - 4);
1073 # print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
1074 while ($step++ <= $steps)
1076 # calculate $a * $f + 1
1077 $A = $MBI->_mul($A, $F);
1078 $A = $MBI->_inc($A);
1080 $F = $MBI->_inc($F);
1082 # compute $B as factorial of $steps (this is faster than doing it manually)
1083 my $B = $MBI->_fac($MBI->_new($steps));
1085 # print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
1087 # compute A/B with $scale digits in the result (truncate, not round)
1088 $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
1089 $A = $MBI->_div( $A, $B );
1094 $x->{_e} = $MBI->_new($scale);
1097 # $x contains now an estimate of e, with some surplus digits, so we can round
1098 if (!$x_org->is_one())
1100 # raise $x to the wanted power and round it in one step:
1101 $x->bpow($x_org, @params);
1105 # else just round the already computed result
1106 delete $x->{_a}; delete $x->{_p};
1107 # shortcut to not run through _find_round_parameters again
1108 if (defined $params[0])
1110 $x->bround($params[0],$params[2]); # then round accordingly
1114 $x->bfround($params[1],$params[2]); # then round accordingly
1119 # clear a/p after round, since user did not request it
1120 delete $x->{_a}; delete $x->{_p};
1123 $$abr = $ab; $$pbr = $pb;
1125 $x; # return modified $x
1130 # internal log function to calculate ln() based on Taylor series.
1131 # Modifies $x in place.
1132 my ($self,$x,$scale) = @_;
1134 # in case of $x == 1, result is 0
1135 return $x->bzero() if $x->is_one();
1137 # XXX TODO: rewrite this in a similiar manner to bexp()
1139 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
1143 # Taylor: | u 1 u^3 1 u^5 |
1144 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
1145 # |_ v 3 v^3 5 v^5 _|
1147 # This takes much more steps to calculate the result and is thus not used
1150 # Taylor: | u 1 u^2 1 u^3 |
1151 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
1152 # |_ x 2 x^2 3 x^3 _|
1154 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
1156 $v = $x->copy(); $v->binc(); # v = x+1
1157 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
1158 $x->bdiv($v,$scale); # first term: u/v
1159 $below = $v->copy();
1161 $u *= $u; $v *= $v; # u^2, v^2
1162 $below->bmul($v); # u^3, v^3
1164 $factor = $self->new(3); $f = $self->new(2);
1166 my $steps = 0 if DEBUG;
1167 $limit = $self->new("1E-". ($scale-1));
1170 # we calculate the next term, and add it to the last
1171 # when the next term is below our limit, it won't affect the outcome
1172 # anymore, so we stop
1174 # calculating the next term simple from over/below will result in quite
1175 # a time hog if the input has many digits, since over and below will
1176 # accumulate more and more digits, and the result will also have many
1177 # digits, but in the end it is rounded to $scale digits anyway. So if we
1178 # round $over and $below first, we save a lot of time for the division
1179 # (not with log(1.2345), but try log (123**123) to see what I mean. This
1180 # can introduce a rounding error if the division result would be f.i.
1181 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
1182 # if we truncated $over and $below we might get 0.12345. Does this matter
1183 # for the end result? So we give $over and $below 4 more digits to be
1184 # on the safe side (unscientific error handling as usual... :+D
1186 $next = $over->copy->bround($scale+4)->bdiv(
1187 $below->copy->bmul($factor)->bround($scale+4),
1191 ## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
1193 last if $next->bacmp($limit) <= 0;
1195 delete $next->{_a}; delete $next->{_p};
1197 # calculate things for the next term
1198 $over *= $u; $below *= $v; $factor->badd($f);
1201 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
1204 print "took $steps steps\n" if DEBUG;
1205 $x->bmul($f); # $x *= 2
1210 # Internal log function based on reducing input to the range of 0.1 .. 9.99
1211 # and then "correcting" the result to the proper one. Modifies $x in place.
1212 my ($self,$x,$scale) = @_;
1214 # Taking blog() from numbers greater than 10 takes a *very long* time, so we
1215 # break the computation down into parts based on the observation that:
1216 # blog(X*Y) = blog(X) + blog(Y)
1217 # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
1218 # $x is the faster it gets. Since 2*$x takes about 10 times as
1219 # long, we make it faster by about a factor of 100 by dividing $x by 10.
1221 # The same observation is valid for numbers smaller than 0.1, e.g. computing
1222 # log(1) is fastest, and the further away we get from 1, the longer it takes.
1223 # So we also 'break' this down by multiplying $x with 10 and subtract the
1224 # log(10) afterwards to get the correct result.
1226 # To get $x even closer to 1, we also divide by 2 and then use log(2) to
1227 # correct for this. For instance if $x is 2.4, we use the formula:
1228 # blog(2.4 * 2) == blog (1.2) + blog(2)
1229 # and thus calculate only blog(1.2) and blog(2), which is faster in total
1230 # than calculating blog(2.4).
1232 # In addition, the values for blog(2) and blog(10) are cached.
1234 # Calculate nr of digits before dot:
1235 my $dbd = $MBI->_num($x->{_e});
1236 $dbd = -$dbd if $x->{_es} eq '-';
1237 $dbd += $MBI->_len($x->{_m});
1239 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
1240 # infinite recursion
1242 my $calc = 1; # do some calculation?
1244 # disable the shortcut for 10, since we need log(10) and this would recurse
1246 if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
1248 $dbd = 0; # disable shortcut
1249 # we can use the cached value in these cases
1250 if ($scale <= $LOG_10_A)
1252 $x->bzero(); $x->badd($LOG_10); # modify $x in place
1253 $calc = 0; # no need to calc, but round
1255 # if we can't use the shortcut, we continue normally
1259 # disable the shortcut for 2, since we maybe have it cached
1260 if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
1262 $dbd = 0; # disable shortcut
1263 # we can use the cached value in these cases
1264 if ($scale <= $LOG_2_A)
1266 $x->bzero(); $x->badd($LOG_2); # modify $x in place
1267 $calc = 0; # no need to calc, but round
1269 # if we can't use the shortcut, we continue normally
1273 # if $x = 0.1, we know the result must be 0-log(10)
1274 if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
1275 $MBI->_is_one($x->{_m}))
1277 $dbd = 0; # disable shortcut
1278 # we can use the cached value in these cases
1279 if ($scale <= $LOG_10_A)
1281 $x->bzero(); $x->bsub($LOG_10);
1282 $calc = 0; # no need to calc, but round
1286 return if $calc == 0; # already have the result
1288 # default: these correction factors are undef and thus not used
1289 my $l_10; # value of ln(10) to A of $scale
1290 my $l_2; # value of ln(2) to A of $scale
1292 my $two = $self->new(2);
1294 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1295 # so don't do this shortcut for 1 or 0
1296 if (($dbd > 1) || ($dbd < 0))
1298 # convert our cached value to an object if not already (avoid doing this
1299 # at import() time, since not everybody needs this)
1300 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1302 #print "x = $x, dbd = $dbd, calc = $calc\n";
1303 # got more than one digit before the dot, or more than one zero after the
1305 # log(123) == log(1.23) + log(10) * 2
1306 # log(0.0123) == log(1.23) - log(10) * 2
1308 if ($scale <= $LOG_10_A)
1311 $l_10 = $LOG_10->copy(); # copy for mul
1315 # else: slower, compute and cache result
1316 # also disable downgrade for this code path
1317 local $Math::BigFloat::downgrade = undef;
1319 # shorten the time to calculate log(10) based on the following:
1320 # log(1.25 * 8) = log(1.25) + log(8)
1321 # = log(1.25) + log(2) + log(2) + log(2)
1323 # first get $l_2 (and possible compute and cache log(2))
1324 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1325 if ($scale <= $LOG_2_A)
1328 $l_2 = $LOG_2->copy(); # copy() for the mul below
1332 # else: slower, compute and cache result
1333 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1334 $LOG_2 = $l_2->copy(); # cache the result for later
1335 # the copy() is for mul below
1339 # now calculate log(1.25):
1340 $l_10 = $self->new('1.25'); $self->_log($l_10, $scale); # scale+4, actually
1342 # log(1.25) + log(2) + log(2) + log(2):
1346 $LOG_10 = $l_10->copy(); # cache the result for later
1347 # the copy() is for mul below
1350 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
1351 $l_10->bmul( $self->new($dbd)); # log(10) * (digits_before_dot-1)
1358 ($x->{_e}, $x->{_es}) =
1359 _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
1363 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1365 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1366 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1368 $HALF = $self->new($HALF) unless ref($HALF);
1370 my $twos = 0; # default: none (0 times)
1371 while ($x->bacmp($HALF) <= 0) # X <= 0.5
1373 $twos--; $x->bmul($two);
1375 while ($x->bacmp($two) >= 0) # X >= 2
1377 $twos++; $x->bdiv($two,$scale+4); # keep all digits
1379 # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
1380 # So calculate correction factor based on ln(2):
1383 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1384 if ($scale <= $LOG_2_A)
1387 $l_2 = $LOG_2->copy(); # copy() for the mul below
1391 # else: slower, compute and cache result
1392 # also disable downgrade for this code path
1393 local $Math::BigFloat::downgrade = undef;
1394 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1395 $LOG_2 = $l_2->copy(); # cache the result for later
1396 # the copy() is for mul below
1399 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
1402 $self->_log($x,$scale); # need to do the "normal" way
1403 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
1404 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
1406 # all done, $x contains now the result
1412 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1413 # does not modify arguments, but returns new object
1414 # Lowest Common Multiplicator
1416 my ($self,@arg) = objectify(0,@_);
1417 my $x = $self->new(shift @arg);
1418 while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); }
1424 # (BINT or num_str, BINT or num_str) return BINT
1425 # does not modify arguments, but returns new object
1428 $y = __PACKAGE__->new($y) if !ref($y);
1430 my $x = $y->copy()->babs(); # keep arguments
1432 return $x->bnan() if $x->{sign} !~ /^[+-]$/ # x NaN?
1433 || !$x->is_int(); # only for integers now
1437 my $t = shift; $t = $self->new($t) if !ref($t);
1438 $y = $t->copy()->babs();
1440 return $x->bnan() if $y->{sign} !~ /^[+-]$/ # y NaN?
1441 || !$y->is_int(); # only for integers now
1443 # greatest common divisor
1444 while (! $y->is_zero())
1446 ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1449 last if $x->is_one();
1454 ##############################################################################
1458 # Internal helper sub to take two positive integers and their signs and
1459 # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1460 # output ($CALC,('+'|'-'))
1461 my ($x,$y,$xs,$ys) = @_;
1463 # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1466 $x = $MBI->_add ($x, $y ); # a+b
1467 # the sign follows $xs
1471 my $a = $MBI->_acmp($x,$y);
1474 $x = $MBI->_sub ($x , $y); # abs sub
1478 $x = $MBI->_zero(); # result is 0
1483 $x = $MBI->_sub ( $y, $x, 1 ); # abs sub
1491 # Internal helper sub to take two positive integers and their signs and
1492 # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1493 # output ($CALC,('+'|'-'))
1494 my ($x,$y,$xs,$ys) = @_;
1498 _e_add($x,$y,$xs,$ys); # call add (does subtract now)
1501 ###############################################################################
1502 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1506 # return true if arg (BFLOAT or num_str) is an integer
1507 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1509 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1510 $x->{_es} eq '+'; # 1e-1 => no integer
1516 # return true if arg (BFLOAT or num_str) is zero
1517 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1519 return 1 if $x->{sign} eq '+' && $MBI->_is_zero($x->{_m});
1525 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1526 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1528 $sign = '+' if !defined $sign || $sign ne '-';
1530 if ($x->{sign} eq $sign &&
1531 $MBI->_is_zero($x->{_e}) && $MBI->_is_one($x->{_m}));
1537 # return true if arg (BFLOAT or num_str) is odd or false if even
1538 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1540 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1541 ($MBI->_is_zero($x->{_e}) && $MBI->_is_odd($x->{_m}));
1547 # return true if arg (BINT or num_str) is even or false if odd
1548 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1550 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
1551 return 1 if ($x->{_es} eq '+' # 123.45 is never
1552 && $MBI->_is_even($x->{_m})); # but 1200 is
1558 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
1559 # (BINT or num_str, BINT or num_str) return BINT
1562 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1563 # objectify is costly, so avoid it
1564 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1566 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1569 return $x if $x->modify('bmul');
1571 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1574 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1576 return $x->bnan() if $x->is_zero() || $y->is_zero();
1577 # result will always be +-inf:
1578 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1579 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1580 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1581 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1582 return $x->binf('-');
1585 return $x->bzero() if $x->is_zero() || $y->is_zero();
1587 return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
1588 ((!$x->isa($self)) || (!$y->isa($self)));
1590 # aEb * cEd = (a*c)E(b+d)
1591 $MBI->_mul($x->{_m},$y->{_m});
1592 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1595 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1596 return $x->bnorm()->round($a,$p,$r,$y);
1601 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1602 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1605 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1606 # objectify is costly, so avoid it
1607 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1609 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1612 return $x if $x->modify('bdiv');
1614 return $self->_div_inf($x,$y)
1615 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1617 # x== 0 # also: or y == 1 or y == -1
1618 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1621 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1623 # we need to limit the accuracy to protect against overflow
1625 my (@params,$scale);
1626 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1628 return $x if $x->is_nan(); # error in _find_round_parameters?
1630 # no rounding at all, so must use fallback
1631 if (scalar @params == 0)
1633 # simulate old behaviour
1634 $params[0] = $self->div_scale(); # and round to it as accuracy
1635 $scale = $params[0]+4; # at least four more for proper round
1636 $params[2] = $r; # round mode by caller or undef
1637 $fallback = 1; # to clear a/p afterwards
1641 # the 4 below is empirical, and there might be cases where it is not
1643 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1646 my $rem; $rem = $self->bzero() if wantarray;
1648 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1650 my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
1651 $scale = $lx if $lx > $scale;
1652 $scale = $ly if $ly > $scale;
1653 my $diff = $ly - $lx;
1654 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
1656 # already handled inf/NaN/-inf above:
1658 # check that $y is not 1 nor -1 and cache the result:
1659 my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));
1661 # flipping the sign of $y will also flip the sign of $x for the special
1662 # case of $x->bsub($x); so we can catch it below:
1663 my $xsign = $x->{sign};
1664 $y->{sign} =~ tr/+-/-+/;
1666 if ($xsign ne $x->{sign})
1668 # special case of $x /= $x results in 1
1669 $x->bone(); # "fixes" also sign of $y, since $x is $y
1673 # correct $y's sign again
1674 $y->{sign} =~ tr/+-/-+/;
1675 # continue with normal div code:
1677 # make copy of $x in case of list context for later reminder calculation
1678 if (wantarray && $y_not_one)
1683 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1685 # check for / +-1 ( +/- 1E0)
1688 # promote BigInts and it's subclasses (except when already a BigFloat)
1689 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1691 # calculate the result to $scale digits and then round it
1692 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1693 $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1694 $MBI->_div ($x->{_m},$y->{_m}); # a/c
1696 # correct exponent of $x
1697 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1698 # correct for 10**scale
1699 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1700 $x->bnorm(); # remove trailing 0's
1702 } # ende else $x != $y
1704 # shortcut to not run through _find_round_parameters again
1705 if (defined $params[0])
1707 delete $x->{_a}; # clear before round
1708 $x->bround($params[0],$params[2]); # then round accordingly
1712 delete $x->{_p}; # clear before round
1713 $x->bfround($params[1],$params[2]); # then round accordingly
1717 # clear a/p after round, since user did not request it
1718 delete $x->{_a}; delete $x->{_p};
1725 $rem->bmod($y,@params); # copy already done
1729 # clear a/p after round, since user did not request it
1730 delete $rem->{_a}; delete $rem->{_p};
1739 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1742 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1743 # objectify is costly, so avoid it
1744 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1746 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1749 return $x if $x->modify('bmod');
1751 # handle NaN, inf, -inf
1752 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1754 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1755 $x->{sign} = $re->{sign};
1756 $x->{_e} = $re->{_e};
1757 $x->{_m} = $re->{_m};
1758 return $x->round($a,$p,$r,$y);
1762 return $x->bnan() if $x->is_zero();
1766 return $x->bzero() if $x->is_zero()
1768 # check that $y == +1 or $y == -1:
1769 ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));
1771 my $cmp = $x->bacmp($y); # equal or $x < $y?
1772 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1774 # only $y of the operands negative?
1775 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1777 $x->{sign} = $y->{sign}; # calc sign first
1778 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1780 my $ym = $MBI->_copy($y->{_m});
1783 $MBI->_lsft( $ym, $y->{_e}, 10)
1784 if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1786 # if $y has digits after dot
1787 my $shifty = 0; # correct _e of $x by this
1788 if ($y->{_es} eq '-') # has digits after dot
1790 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1791 $shifty = $MBI->_num($y->{_e}); # no more digits after dot
1792 $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1794 # $ym is now mantissa of $y based on exponent 0
1796 my $shiftx = 0; # correct _e of $x by this
1797 if ($x->{_es} eq '-') # has digits after dot
1799 # 123.4 % 20 => 1234 % 200
1800 $shiftx = $MBI->_num($x->{_e}); # no more digits after dot
1801 $MBI->_lsft($ym, $x->{_e}, 10); # 123 => 1230
1803 # 123e1 % 20 => 1230 % 20
1804 if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1806 $MBI->_lsft( $x->{_m}, $x->{_e},10); # es => '+' here
1809 $x->{_e} = $MBI->_new($shiftx);
1811 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1812 $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1814 # now mantissas are equalized, exponent of $x is adjusted, so calc result
1816 $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1818 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1821 if ($neg != 0) # one of them negative => correct in place
1824 $x->{_m} = $r->{_m};
1825 $x->{_e} = $r->{_e};
1826 $x->{_es} = $r->{_es};
1827 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1831 $x->round($a,$p,$r,$y); # round and return
1836 # calculate $y'th root of $x
1839 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1840 # objectify is costly, so avoid it
1841 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1843 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1846 return $x if $x->modify('broot');
1848 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1849 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1850 $y->{sign} !~ /^\+$/;
1852 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1854 # we need to limit the accuracy to protect against overflow
1856 my (@params,$scale);
1857 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1859 return $x if $x->is_nan(); # error in _find_round_parameters?
1861 # no rounding at all, so must use fallback
1862 if (scalar @params == 0)
1864 # simulate old behaviour
1865 $params[0] = $self->div_scale(); # and round to it as accuracy
1866 $scale = $params[0]+4; # at least four more for proper round
1867 $params[2] = $r; # iound mode by caller or undef
1868 $fallback = 1; # to clear a/p afterwards
1872 # the 4 below is empirical, and there might be cases where it is not
1874 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1877 # when user set globals, they would interfere with our calculation, so
1878 # disable them and later re-enable them
1880 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1881 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1882 # we also need to disable any set A or P on $x (_find_round_parameters took
1883 # them already into account), since these would interfere, too
1884 delete $x->{_a}; delete $x->{_p};
1885 # need to disable $upgrade in BigInt, to avoid deep recursion
1886 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1888 # remember sign and make $x positive, since -4 ** (1/2) => -2
1889 my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
1892 if ($y->isa('Math::BigFloat'))
1894 $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
1898 $is_two = ($y == 2);
1901 # normal square root if $y == 2:
1904 $x->bsqrt($scale+4);
1906 elsif ($y->is_one('-'))
1909 my $u = $self->bone()->bdiv($x,$scale);
1910 # copy private parts over
1911 $x->{_m} = $u->{_m};
1912 $x->{_e} = $u->{_e};
1913 $x->{_es} = $u->{_es};
1917 # calculate the broot() as integer result first, and if it fits, return
1918 # it rightaway (but only if $x and $y are integer):
1920 my $done = 0; # not yet
1921 if ($y->is_int() && $x->is_int())
1923 my $i = $MBI->_copy( $x->{_m} );
1924 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1925 my $int = Math::BigInt->bzero();
1927 $int->broot($y->as_number());
1929 if ($int->copy()->bpow($y) == $x)
1931 # found result, return it
1932 $x->{_m} = $int->{value};
1933 $x->{_e} = $MBI->_zero();
1941 my $u = $self->bone()->bdiv($y,$scale+4);
1942 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
1943 $x->bpow($u,$scale+4); # el cheapo
1946 $x->bneg() if $sign == 1;
1948 # shortcut to not run through _find_round_parameters again
1949 if (defined $params[0])
1951 $x->bround($params[0],$params[2]); # then round accordingly
1955 $x->bfround($params[1],$params[2]); # then round accordingly
1959 # clear a/p after round, since user did not request it
1960 delete $x->{_a}; delete $x->{_p};
1963 $$abr = $ab; $$pbr = $pb;
1969 # calculate square root
1970 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1972 return $x if $x->modify('bsqrt');
1974 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
1975 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
1976 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
1978 # we need to limit the accuracy to protect against overflow
1980 my (@params,$scale);
1981 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1983 return $x if $x->is_nan(); # error in _find_round_parameters?
1985 # no rounding at all, so must use fallback
1986 if (scalar @params == 0)
1988 # simulate old behaviour
1989 $params[0] = $self->div_scale(); # and round to it as accuracy
1990 $scale = $params[0]+4; # at least four more for proper round
1991 $params[2] = $r; # round mode by caller or undef
1992 $fallback = 1; # to clear a/p afterwards
1996 # the 4 below is empirical, and there might be cases where it is not
1998 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2001 # when user set globals, they would interfere with our calculation, so
2002 # disable them and later re-enable them
2004 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2005 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2006 # we also need to disable any set A or P on $x (_find_round_parameters took
2007 # them already into account), since these would interfere, too
2008 delete $x->{_a}; delete $x->{_p};
2009 # need to disable $upgrade in BigInt, to avoid deep recursion
2010 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2012 my $i = $MBI->_copy( $x->{_m} );
2013 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2014 my $xas = Math::BigInt->bzero();
2017 my $gs = $xas->copy()->bsqrt(); # some guess
2019 if (($x->{_es} ne '-') # guess can't be accurate if there are
2020 # digits after the dot
2021 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
2023 # exact result, copy result over to keep $x
2024 $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
2026 # shortcut to not run through _find_round_parameters again
2027 if (defined $params[0])
2029 $x->bround($params[0],$params[2]); # then round accordingly
2033 $x->bfround($params[1],$params[2]); # then round accordingly
2037 # clear a/p after round, since user did not request it
2038 delete $x->{_a}; delete $x->{_p};
2040 # re-enable A and P, upgrade is taken care of by "local"
2041 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
2045 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
2046 # of the result by multipyling the input by 100 and then divide the integer
2047 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
2049 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
2050 my $y1 = $MBI->_copy($x->{_m});
2052 my $length = $MBI->_len($y1);
2054 # Now calculate how many digits the result of sqrt(y1) would have
2055 my $digits = int($length / 2);
2057 # But we need at least $scale digits, so calculate how many are missing
2058 my $shift = $scale - $digits;
2060 # That should never happen (we take care of integer guesses above)
2061 # $shift = 0 if $shift < 0;
2063 # Multiply in steps of 100, by shifting left two times the "missing" digits
2064 my $s2 = $shift * 2;
2066 # We now make sure that $y1 has the same odd or even number of digits than
2067 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
2068 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
2069 # steps of 10. The length of $x does not count, since an even or odd number
2070 # of digits before the dot is not changed by adding an even number of digits
2071 # after the dot (the result is still odd or even digits long).
2072 $s2++ if $MBI->_is_odd($x->{_e});
2074 $MBI->_lsft( $y1, $MBI->_new($s2), 10);
2076 # now take the square root and truncate to integer
2077 $y1 = $MBI->_sqrt($y1);
2079 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
2080 # result, which is than later rounded to the desired scale.
2082 # calculate how many zeros $x had after the '.' (or before it, depending
2083 # on sign of $dat, the result should have half as many:
2084 my $dat = $MBI->_num($x->{_e});
2085 $dat = -$dat if $x->{_es} eq '-';
2090 # no zeros after the dot (e.g. 1.23, 0.49 etc)
2091 # preserve half as many digits before the dot than the input had
2092 # (but round this "up")
2093 $dat = int(($dat+1)/2);
2097 $dat = int(($dat)/2);
2099 $dat -= $MBI->_len($y1);
2103 $x->{_e} = $MBI->_new( $dat );
2108 $x->{_e} = $MBI->_new( $dat );
2114 # shortcut to not run through _find_round_parameters again
2115 if (defined $params[0])
2117 $x->bround($params[0],$params[2]); # then round accordingly
2121 $x->bfround($params[1],$params[2]); # then round accordingly
2125 # clear a/p after round, since user did not request it
2126 delete $x->{_a}; delete $x->{_p};
2129 $$abr = $ab; $$pbr = $pb;
2135 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2136 # compute factorial number, modifies first argument
2139 my ($self,$x,@r) = (ref($_[0]),@_);
2140 # objectify is costly, so avoid it
2141 ($self,$x,@r) = objectify(1,@_) if !ref($x);
2144 return $x if $x->modify('bfac') || $x->{sign} eq '+inf';
2147 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
2148 ($x->{_es} ne '+')); # digits after dot?
2150 # use BigInt's bfac() for faster calc
2151 if (! $MBI->_is_zero($x->{_e}))
2153 $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
2154 $x->{_e} = $MBI->_zero(); # normalize
2157 $MBI->_fac($x->{_m}); # calculate factorial
2158 $x->bnorm()->round(@r); # norm again and round result
2163 # Calculate a power where $y is a non-integer, like 2 ** 0.5
2164 my ($x,$y,$a,$p,$r) = @_;
2167 # if $y == 0.5, it is sqrt($x)
2168 $HALF = $self->new($HALF) unless ref($HALF);
2169 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp($HALF) == 0;
2172 # a ** x == e ** (x * ln a)
2176 # Taylor: | u u^2 u^3 |
2177 # x ** y = 1 + | --- + --- + ----- + ... |
2180 # we need to limit the accuracy to protect against overflow
2182 my ($scale,@params);
2183 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2185 return $x if $x->is_nan(); # error in _find_round_parameters?
2187 # no rounding at all, so must use fallback
2188 if (scalar @params == 0)
2190 # simulate old behaviour
2191 $params[0] = $self->div_scale(); # and round to it as accuracy
2192 $params[1] = undef; # disable P
2193 $scale = $params[0]+4; # at least four more for proper round
2194 $params[2] = $r; # round mode by caller or undef
2195 $fallback = 1; # to clear a/p afterwards
2199 # the 4 below is empirical, and there might be cases where it is not
2201 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2204 # when user set globals, they would interfere with our calculation, so
2205 # disable them and later re-enable them
2207 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2208 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2209 # we also need to disable any set A or P on $x (_find_round_parameters took
2210 # them already into account), since these would interfere, too
2211 delete $x->{_a}; delete $x->{_p};
2212 # need to disable $upgrade in BigInt, to avoid deep recursion
2213 local $Math::BigInt::upgrade = undef;
2215 my ($limit,$v,$u,$below,$factor,$next,$over);
2217 $u = $x->copy()->blog(undef,$scale)->bmul($y);
2218 $v = $self->bone(); # 1
2219 $factor = $self->new(2); # 2
2220 $x->bone(); # first term: 1
2222 $below = $v->copy();
2225 $limit = $self->new("1E-". ($scale-1));
2229 # we calculate the next term, and add it to the last
2230 # when the next term is below our limit, it won't affect the outcome
2231 # anymore, so we stop:
2232 $next = $over->copy()->bdiv($below,$scale);
2233 last if $next->bacmp($limit) <= 0;
2235 # calculate things for the next term
2236 $over *= $u; $below *= $factor; $factor->binc();
2238 last if $x->{sign} !~ /^[-+]$/;
2243 # shortcut to not run through _find_round_parameters again
2244 if (defined $params[0])
2246 $x->bround($params[0],$params[2]); # then round accordingly
2250 $x->bfround($params[1],$params[2]); # then round accordingly
2254 # clear a/p after round, since user did not request it
2255 delete $x->{_a}; delete $x->{_p};
2258 $$abr = $ab; $$pbr = $pb;
2264 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2265 # compute power of two numbers, second arg is used as integer
2266 # modifies first argument
2269 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2270 # objectify is costly, so avoid it
2271 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2273 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2276 return $x if $x->modify('bpow');
2278 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
2279 return $x if $x->{sign} =~ /^[+-]inf$/;
2281 # cache the result of is_zero
2282 my $y_is_zero = $y->is_zero();
2283 return $x->bone() if $y_is_zero;
2284 return $x if $x->is_one() || $y->is_one();
2286 my $x_is_zero = $x->is_zero();
2287 return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int(); # non-integer power
2289 my $y1 = $y->as_number()->{value}; # make MBI part
2292 if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
2294 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
2295 return $MBI->_is_odd($y1) ? $x : $x->babs(1);
2299 return $x->bone() if $y_is_zero;
2300 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
2301 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
2306 $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
2308 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
2309 $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
2310 $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
2312 $x->{sign} = $new_sign;
2314 if ($y->{sign} eq '-')
2316 # modify $x in place!
2317 my $z = $x->copy(); $x->bone();
2318 return scalar $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
2320 $x->round($a,$p,$r,$y);
2323 ###############################################################################
2324 # rounding functions
2328 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
2329 # $n == 0 means round to integer
2330 # expects and returns normalized numbers!
2331 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2333 my ($scale,$mode) = $x->_scale_p(@_);
2334 return $x if !defined $scale || $x->modify('bfround'); # no-op
2336 # never round a 0, +-inf, NaN
2339 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
2342 return $x if $x->{sign} !~ /^[+-]$/;
2344 # don't round if x already has lower precision
2345 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
2347 $x->{_p} = $scale; # remember round in any case
2348 delete $x->{_a}; # and clear A
2351 # round right from the '.'
2353 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round
2355 $scale = -$scale; # positive for simplicity
2356 my $len = $MBI->_len($x->{_m}); # length of mantissa
2358 # the following poses a restriction on _e, but if _e is bigger than a
2359 # scalar, you got other problems (memory etc) anyway
2360 my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e}))); # digits after dot
2361 my $zad = 0; # zeros after dot
2362 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
2364 # p rint "scale $scale dad $dad zad $zad len $len\n";
2365 # number bsstr len zad dad
2366 # 0.123 123e-3 3 0 3
2367 # 0.0123 123e-4 3 1 4
2370 # 1.2345 12345e-4 5 0 4
2372 # do not round after/right of the $dad
2373 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
2375 # round to zero if rounding inside the $zad, but not for last zero like:
2376 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
2377 return $x->bzero() if $scale < $zad;
2378 if ($scale == $zad) # for 0.006, scale -3 and trunc
2384 # adjust round-point to be inside mantissa
2387 $scale = $scale-$zad;
2391 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
2392 $scale = $dbd+$scale;
2398 # round left from the '.'
2400 # 123 => 100 means length(123) = 3 - $scale (2) => 1
2402 my $dbt = $MBI->_len($x->{_m});
2404 my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
2405 # should be the same, so treat it as this
2406 $scale = 1 if $scale == 0;
2407 # shortcut if already integer
2408 return $x if $scale == 1 && $dbt <= $dbd;
2409 # maximum digits before dot
2414 # not enough digits before dot, so round to zero
2417 elsif ( $scale == $dbd )
2424 $scale = $dbd - $scale;
2427 # pass sign to bround for rounding modes '+inf' and '-inf'
2428 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2429 $m->bround($scale,$mode);
2430 $x->{_m} = $m->{value}; # get our mantissa back
2436 # accuracy: preserve $N digits, and overwrite the rest with 0's
2437 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2439 if (($_[0] || 0) < 0)
2441 require Carp; Carp::croak ('bround() needs positive accuracy');
2444 my ($scale,$mode) = $x->_scale_a(@_);
2445 return $x if !defined $scale || $x->modify('bround'); # no-op
2447 # scale is now either $x->{_a}, $accuracy, or the user parameter
2448 # test whether $x already has lower accuracy, do nothing in this case
2449 # but do round if the accuracy is the same, since a math operation might
2450 # want to round a number with A=5 to 5 digits afterwards again
2451 return $x if defined $x->{_a} && $x->{_a} < $scale;
2453 # scale < 0 makes no sense
2454 # scale == 0 => keep all digits
2455 # never round a +-inf, NaN
2456 return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
2458 # 1: never round a 0
2459 # 2: if we should keep more digits than the mantissa has, do nothing
2460 if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
2462 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
2466 # pass sign to bround for '+inf' and '-inf' rounding modes
2467 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2469 $m->bround($scale,$mode); # round mantissa
2470 $x->{_m} = $m->{value}; # get our mantissa back
2471 $x->{_a} = $scale; # remember rounding
2472 delete $x->{_p}; # and clear P
2473 $x->bnorm(); # del trailing zeros gen. by bround()
2478 # return integer less or equal then $x
2479 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2481 return $x if $x->modify('bfloor');
2483 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2485 # if $x has digits after dot
2486 if ($x->{_es} eq '-')
2488 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2489 $x->{_e} = $MBI->_zero(); # trunc/norm
2490 $x->{_es} = '+'; # abs e
2491 $MBI->_inc($x->{_m}) if $x->{sign} eq '-'; # increment if negative
2493 $x->round($a,$p,$r);
2498 # return integer greater or equal then $x
2499 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2501 return $x if $x->modify('bceil');
2502 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2504 # if $x has digits after dot
2505 if ($x->{_es} eq '-')
2507 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2508 $x->{_e} = $MBI->_zero(); # trunc/norm
2509 $x->{_es} = '+'; # abs e
2510 $MBI->_inc($x->{_m}) if $x->{sign} eq '+'; # increment if positive
2512 $x->round($a,$p,$r);
2517 # shift right by $y (divide by power of $n)
2520 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2521 # objectify is costly, so avoid it
2522 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2524 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2527 return $x if $x->modify('brsft');
2528 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2530 $n = 2 if !defined $n; $n = $self->new($n);
2533 return $x->blsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
2535 # the following call to bdiv() will return either quo or (quo,reminder):
2536 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
2541 # shift left by $y (multiply by power of $n)
2544 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2545 # objectify is costly, so avoid it
2546 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2548 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2551 return $x if $x->modify('blsft');
2552 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2554 $n = 2 if !defined $n; $n = $self->new($n);
2557 return $x->brsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
2559 $x->bmul($n->bpow($y),$a,$p,$r,$y);
2562 ###############################################################################
2566 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
2571 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
2572 # or falling back to MBI::bxxx()
2573 my $name = $AUTOLOAD;
2575 $name =~ s/(.*):://; # split package
2576 my $c = $1 || $class;
2578 $c->import() if $IMPORT == 0;
2579 if (!_method_alias($name))
2583 # delayed load of Carp and avoid recursion
2585 Carp::croak ("$c: Can't call a method without name");
2587 if (!_method_hand_up($name))
2589 # delayed load of Carp and avoid recursion
2591 Carp::croak ("Can't call $c\-\>$name, not a valid method");
2593 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
2595 return &{"Math::BigInt"."::$name"}(@_);
2597 my $bname = $name; $bname =~ s/^f/b/;
2605 # return a copy of the exponent
2606 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2608 if ($x->{sign} !~ /^[+-]$/)
2610 my $s = $x->{sign}; $s =~ s/^[+-]//;
2611 return Math::BigInt->new($s); # -inf, +inf => +inf
2613 Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
2618 # return a copy of the mantissa
2619 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2621 if ($x->{sign} !~ /^[+-]$/)
2623 my $s = $x->{sign}; $s =~ s/^[+]//;
2624 return Math::BigInt->new($s); # -inf, +inf => +inf
2626 my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
2627 $m->bneg() if $x->{sign} eq '-';
2634 # return a copy of both the exponent and the mantissa
2635 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2637 if ($x->{sign} !~ /^[+-]$/)
2639 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2640 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2642 my $m = Math::BigInt->bzero();
2643 $m->{value} = $MBI->_copy($x->{_m});
2644 $m->bneg() if $x->{sign} eq '-';
2645 ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
2648 ##############################################################################
2649 # private stuff (internal use only)
2655 my $lib = ''; my @a;
2656 my $lib_kind = 'try';
2658 for ( my $i = 0; $i < $l ; $i++)
2660 if ( $_[$i] eq ':constant' )
2662 # This causes overlord er load to step in. 'binary' and 'integer'
2663 # are handled by BigInt.
2664 overload::constant float => sub { $self->new(shift); };
2666 elsif ($_[$i] eq 'upgrade')
2668 # this causes upgrading
2669 $upgrade = $_[$i+1]; # or undef to disable
2672 elsif ($_[$i] eq 'downgrade')
2674 # this causes downgrading
2675 $downgrade = $_[$i+1]; # or undef to disable
2678 elsif ($_[$i] =~ /^(lib|try|only)\z/)
2680 # alternative library
2681 $lib = $_[$i+1] || ''; # default Calc
2682 $lib_kind = $1; # lib, try or only
2685 elsif ($_[$i] eq 'with')
2687 # alternative class for our private parts()
2688 # XXX: no longer supported
2689 # $MBI = $_[$i+1] || 'Math::BigInt';
2698 $lib =~ tr/a-zA-Z0-9,://cd; # restrict to sane characters
2699 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2700 my $mbilib = eval { Math::BigInt->config()->{lib} };
2701 if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
2703 # MBI already loaded
2704 Math::BigInt->import( $lib_kind, "$lib,$mbilib", 'objectify');
2708 # MBI not loaded, or with ne "Math::BigInt::Calc"
2709 $lib .= ",$mbilib" if defined $mbilib;
2710 $lib =~ s/^,//; # don't leave empty
2712 # replacement library can handle lib statement, but also could ignore it
2714 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2715 # used in the same script, or eval inside import(). So we require MBI:
2716 require Math::BigInt;
2717 Math::BigInt->import( $lib_kind => $lib, 'objectify' );
2721 require Carp; Carp::croak ("Couldn't load $lib: $! $@");
2723 # find out which one was actually loaded
2724 $MBI = Math::BigInt->config()->{lib};
2726 # register us with MBI to get notified of future lib changes
2727 Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
2729 # any non :constant stuff is handled by our parent, Exporter
2730 # even if @_ is empty, to give it a chance
2731 $self->SUPER::import(@a); # for subclasses
2732 $self->export_to_level(1,$self,@a); # need this, too
2737 # adjust m and e so that m is smallest possible
2738 # round number according to accuracy and precision settings
2739 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
2741 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2743 my $zeros = $MBI->_zeros($x->{_m}); # correct for trailing zeros
2746 my $z = $MBI->_new($zeros);
2747 $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
2748 if ($x->{_es} eq '-')
2750 if ($MBI->_acmp($x->{_e},$z) >= 0)
2752 $x->{_e} = $MBI->_sub ($x->{_e}, $z);
2753 $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
2757 $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
2763 $x->{_e} = $MBI->_add ($x->{_e}, $z);
2768 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
2769 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
2770 $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
2771 if $MBI->_is_zero($x->{_m});
2774 $x; # MBI bnorm is no-op, so dont call it
2777 ##############################################################################
2781 # return number as hexadecimal string (only for integers defined)
2782 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2784 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2785 return '0x0' if $x->is_zero();
2787 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2789 my $z = $MBI->_copy($x->{_m});
2790 if (! $MBI->_is_zero($x->{_e})) # > 0
2792 $MBI->_lsft( $z, $x->{_e},10);
2794 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2800 # return number as binary digit string (only for integers defined)
2801 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2803 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2804 return '0b0' if $x->is_zero();
2806 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2808 my $z = $MBI->_copy($x->{_m});
2809 if (! $MBI->_is_zero($x->{_e})) # > 0
2811 $MBI->_lsft( $z, $x->{_e},10);
2813 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2819 # return number as octal digit string (only for integers defined)
2820 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2822 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2823 return '0' if $x->is_zero();
2825 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
2827 my $z = $MBI->_copy($x->{_m});
2828 if (! $MBI->_is_zero($x->{_e})) # > 0
2830 $MBI->_lsft( $z, $x->{_e},10);
2832 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2838 # return copy as a bigint representation of this BigFloat number
2839 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2841 return $x if $x->modify('as_number');
2843 my $z = $MBI->_copy($x->{_m});
2844 if ($x->{_es} eq '-') # < 0
2846 $MBI->_rsft( $z, $x->{_e},10);
2848 elsif (! $MBI->_is_zero($x->{_e})) # > 0
2850 $MBI->_lsft( $z, $x->{_e},10);
2852 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2859 my $class = ref($x) || $x;
2860 $x = $class->new(shift) unless ref($x);
2862 return 1 if $MBI->_is_zero($x->{_m});
2864 my $len = $MBI->_len($x->{_m});
2865 $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
2869 $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
2880 Math::BigFloat - Arbitrary size floating point math package
2887 $x = Math::BigFloat->new($str); # defaults to 0
2888 $nan = Math::BigFloat->bnan(); # create a NotANumber
2889 $zero = Math::BigFloat->bzero(); # create a +0
2890 $inf = Math::BigFloat->binf(); # create a +inf
2891 $inf = Math::BigFloat->binf('-'); # create a -inf
2892 $one = Math::BigFloat->bone(); # create a +1
2893 $one = Math::BigFloat->bone('-'); # create a -1
2896 $x->is_zero(); # true if arg is +0
2897 $x->is_nan(); # true if arg is NaN
2898 $x->is_one(); # true if arg is +1
2899 $x->is_one('-'); # true if arg is -1
2900 $x->is_odd(); # true if odd, false for even
2901 $x->is_even(); # true if even, false for odd
2902 $x->is_pos(); # true if >= 0
2903 $x->is_neg(); # true if < 0
2904 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2906 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2907 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2908 $x->sign(); # return the sign, either +,- or NaN
2909 $x->digit($n); # return the nth digit, counting from right
2910 $x->digit(-$n); # return the nth digit, counting from left
2912 # The following all modify their first argument. If you want to preserve
2913 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2914 # necessary when mixing $a = $b assignments with non-overloaded math.
2917 $x->bzero(); # set $i to 0
2918 $x->bnan(); # set $i to NaN
2919 $x->bone(); # set $x to +1
2920 $x->bone('-'); # set $x to -1
2921 $x->binf(); # set $x to inf
2922 $x->binf('-'); # set $x to -inf
2924 $x->bneg(); # negation
2925 $x->babs(); # absolute value
2926 $x->bnorm(); # normalize (no-op)
2927 $x->bnot(); # two's complement (bit wise not)
2928 $x->binc(); # increment x by 1
2929 $x->bdec(); # decrement x by 1
2931 $x->badd($y); # addition (add $y to $x)
2932 $x->bsub($y); # subtraction (subtract $y from $x)
2933 $x->bmul($y); # multiplication (multiply $x by $y)
2934 $x->bdiv($y); # divide, set $x to quotient
2935 # return (quo,rem) or quo if scalar
2937 $x->bmod($y); # modulus ($x % $y)
2938 $x->bpow($y); # power of arguments ($x ** $y)
2939 $x->blsft($y, $n); # left shift by $y places in base $n
2940 $x->brsft($y, $n); # right shift by $y places in base $n
2941 # returns (quo,rem) or quo if in scalar context
2943 $x->blog(); # logarithm of $x to base e (Euler's number)
2944 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
2945 $x->bexp(); # calculate e ** $x where e is Euler's number
2947 $x->band($y); # bit-wise and
2948 $x->bior($y); # bit-wise inclusive or
2949 $x->bxor($y); # bit-wise exclusive or
2950 $x->bnot(); # bit-wise not (two's complement)
2952 $x->bsqrt(); # calculate square-root
2953 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
2954 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2956 $x->bround($N); # accuracy: preserve $N digits
2957 $x->bfround($N); # precision: round to the $Nth digit
2959 $x->bfloor(); # return integer less or equal than $x
2960 $x->bceil(); # return integer greater or equal than $x
2962 # The following do not modify their arguments:
2964 bgcd(@values); # greatest common divisor
2965 blcm(@values); # lowest common multiplicator
2967 $x->bstr(); # return string
2968 $x->bsstr(); # return string in scientific notation
2970 $x->as_int(); # return $x as BigInt
2971 $x->exponent(); # return exponent as BigInt
2972 $x->mantissa(); # return mantissa as BigInt
2973 $x->parts(); # return (mantissa,exponent) as BigInt
2975 $x->length(); # number of digits (w/o sign and '.')
2976 ($l,$f) = $x->length(); # number of digits, and length of fraction
2978 $x->precision(); # return P of $x (or global, if P of $x undef)
2979 $x->precision($n); # set P of $x to $n
2980 $x->accuracy(); # return A of $x (or global, if A of $x undef)
2981 $x->accuracy($n); # set A $x to $n
2983 # these get/set the appropriate global value for all BigFloat objects
2984 Math::BigFloat->precision(); # Precision
2985 Math::BigFloat->accuracy(); # Accuracy
2986 Math::BigFloat->round_mode(); # rounding mode
2990 All operators (including basic math operations) are overloaded if you
2991 declare your big floating point numbers as
2993 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2995 Operations with overloaded operators preserve the arguments, which is
2996 exactly what you expect.
2998 =head2 Canonical notation
3000 Input to these routines are either BigFloat objects, or strings of the
3001 following four forms:
3015 C</^[+-]\d+E[+-]?\d+$/>
3019 C</^[+-]\d*\.\d+E[+-]?\d+$/>
3023 all with optional leading and trailing zeros and/or spaces. Additionally,
3024 numbers are allowed to have an underscore between any two digits.
3026 Empty strings as well as other illegal numbers results in 'NaN'.
3028 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
3029 are always stored in normalized form. On a string, it creates a BigFloat
3034 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
3036 The string output will always have leading and trailing zeros stripped and drop
3037 a plus sign. C<bstr()> will give you always the form with a decimal point,
3038 while C<bsstr()> (s for scientific) gives you the scientific notation.
3040 Input bstr() bsstr()
3042 ' -123 123 123' '-123123123' '-123123123E0'
3043 '00.0123' '0.0123' '123E-4'
3044 '123.45E-2' '1.2345' '12345E-4'
3045 '10E+3' '10000' '1E4'
3047 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
3048 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
3049 return either undef, <0, 0 or >0 and are suited for sort.
3051 Actual math is done by using the class defined with C<with => Class;> (which
3052 defaults to BigInts) to represent the mantissa and exponent.
3054 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
3055 represent the result when input arguments are not numbers, as well as
3056 the result of dividing by zero.
3058 =head2 C<mantissa()>, C<exponent()> and C<parts()>
3060 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
3061 as BigInts such that:
3063 $m = $x->mantissa();
3064 $e = $x->exponent();
3065 $y = $m * ( 10 ** $e );
3066 print "ok\n" if $x == $y;
3068 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
3070 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
3072 Currently the mantissa is reduced as much as possible, favouring higher
3073 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
3074 This might change in the future, so do not depend on it.
3076 =head2 Accuracy vs. Precision
3078 See also: L<Rounding|Rounding>.
3080 Math::BigFloat supports both precision (rounding to a certain place before or
3081 after the dot) and accuracy (rounding to a certain number of digits). For a
3082 full documentation, examples and tips on these topics please see the large
3083 section about rounding in L<Math::BigInt>.
3085 Since things like C<sqrt(2)> or C<1 / 3> must presented with a limited
3086 accuracy lest a operation consumes all resources, each operation produces
3087 no more than the requested number of digits.
3089 If there is no gloabl precision or accuracy set, B<and> the operation in
3090 question was not called with a requested precision or accuracy, B<and> the
3091 input $x has no accuracy or precision set, then a fallback parameter will
3092 be used. For historical reasons, it is called C<div_scale> and can be accessed
3095 $d = Math::BigFloat->div_scale(); # query
3096 Math::BigFloat->div_scale($n); # set to $n digits
3098 The default value for C<div_scale> is 40.
3100 In case the result of one operation has more digits than specified,
3101 it is rounded. The rounding mode taken is either the default mode, or the one
3102 supplied to the operation after the I<scale>:
3104 $x = Math::BigFloat->new(2);
3105 Math::BigFloat->accuracy(5); # 5 digits max
3106 $y = $x->copy()->bdiv(3); # will give 0.66667
3107 $y = $x->copy()->bdiv(3,6); # will give 0.666667
3108 $y = $x->copy()->bdiv(3,6,undef,'odd'); # will give 0.666667
3109 Math::BigFloat->round_mode('zero');
3110 $y = $x->copy()->bdiv(3,6); # will also give 0.666667
3112 Note that C<< Math::BigFloat->accuracy() >> and C<< Math::BigFloat->precision() >>
3113 set the global variables, and thus B<any> newly created number will be subject
3114 to the global rounding B<immediately>. This means that in the examples above, the
3115 C<3> as argument to C<bdiv()> will also get an accuracy of B<5>.
3117 It is less confusing to either calculate the result fully, and afterwards
3118 round it explicitly, or use the additional parameters to the math
3122 $x = Math::BigFloat->new(2);
3123 $y = $x->copy()->bdiv(3);
3124 print $y->bround(5),"\n"; # will give 0.66667
3129 $x = Math::BigFloat->new(2);
3130 $y = $x->copy()->bdiv(3,5); # will give 0.66667
3137 =item ffround ( +$scale )
3139 Rounds to the $scale'th place left from the '.', counting from the dot.
3140 The first digit is numbered 1.
3142 =item ffround ( -$scale )
3144 Rounds to the $scale'th place right from the '.', counting from the dot.
3148 Rounds to an integer.
3150 =item fround ( +$scale )
3152 Preserves accuracy to $scale digits from the left (aka significant digits)
3153 and pads the rest with zeros. If the number is between 1 and -1, the
3154 significant digits count from the first non-zero after the '.'
3156 =item fround ( -$scale ) and fround ( 0 )
3158 These are effectively no-ops.
3162 All rounding functions take as a second parameter a rounding mode from one of
3163 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
3165 The default rounding mode is 'even'. By using
3166 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
3167 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
3168 no longer supported.
3169 The second parameter to the round functions then overrides the default
3172 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
3173 'trunc' as rounding mode to make it equivalent to:
3178 You can override this by passing the desired rounding mode as parameter to
3181 $x = Math::BigFloat->new(2.5);
3182 $y = $x->as_number('odd'); # $y = 3
3188 $x->accuracy(5); # local for $x
3189 CLASS->accuracy(5); # global for all members of CLASS
3190 # Note: This also applies to new()!
3192 $A = $x->accuracy(); # read out accuracy that affects $x
3193 $A = CLASS->accuracy(); # read out global accuracy
3195 Set or get the global or local accuracy, aka how many significant digits the
3196 results have. If you set a global accuracy, then this also applies to new()!
3198 Warning! The accuracy I<sticks>, e.g. once you created a number under the
3199 influence of C<< CLASS->accuracy($A) >>, all results from math operations with
3200 that number will also be rounded.
3202 In most cases, you should probably round the results explicitly using one of
3203 L<round()>, L<bround()> or L<bfround()> or by passing the desired accuracy
3204 to the math operation as additional parameter:
3206 my $x = Math::BigInt->new(30000);
3207 my $y = Math::BigInt->new(7);
3208 print scalar $x->copy()->bdiv($y, 2); # print 4300
3209 print scalar $x->copy()->bdiv($y)->bround(2); # print 4300
3213 $x->precision(-2); # local for $x, round at the second digit right of the dot
3214 $x->precision(2); # ditto, round at the second digit left of the dot
3216 CLASS->precision(5); # Global for all members of CLASS
3217 # This also applies to new()!
3218 CLASS->precision(-5); # ditto
3220 $P = CLASS->precision(); # read out global precision
3221 $P = $x->precision(); # read out precision that affects $x
3223 Note: You probably want to use L<accuracy()> instead. With L<accuracy> you
3224 set the number of digits each result should have, with L<precision> you
3225 set the place where to round!
3227 =head1 Autocreating constants
3229 After C<use Math::BigFloat ':constant'> all the floating point constants
3230 in the given scope are converted to C<Math::BigFloat>. This conversion
3231 happens at compile time.
3235 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
3237 prints the value of C<2E-100>. Note that without conversion of
3238 constants the expression 2E-100 will be calculated as normal floating point
3241 Please note that ':constant' does not affect integer constants, nor binary
3242 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
3247 Math with the numbers is done (by default) by a module called
3248 Math::BigInt::Calc. This is equivalent to saying:
3250 use Math::BigFloat lib => 'Calc';
3252 You can change this by using:
3254 use Math::BigFloat lib => 'BitVect';
3256 The following would first try to find Math::BigInt::Foo, then
3257 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
3259 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
3261 Calc.pm uses as internal format an array of elements of some decimal base
3262 (usually 1e7, but this might be different for some systems) with the least
3263 significant digit first, while BitVect.pm uses a bit vector of base 2, most
3264 significant bit first. Other modules might use even different means of
3265 representing the numbers. See the respective module documentation for further
3268 Please note that Math::BigFloat does B<not> use the denoted library itself,
3269 but it merely passes the lib argument to Math::BigInt. So, instead of the need
3272 use Math::BigInt lib => 'GMP';
3275 you can roll it all into one line:
3277 use Math::BigFloat lib => 'GMP';
3279 It is also possible to just require Math::BigFloat:
3281 require Math::BigFloat;
3283 This will load the necessary things (like BigInt) when they are needed, and
3286 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
3287 you ever wanted to know about loading a different library.
3289 =head2 Using Math::BigInt::Lite
3291 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
3294 use Math::BigFloat with => 'Math::BigInt::Lite';
3296 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
3297 can combine these if you want. For instance, you may want to use
3298 Math::BigInt objects in your main script, too.
3302 use Math::BigFloat with => 'Math::BigInt::Lite';
3304 Of course, you can combine this with the C<lib> parameter.
3307 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
3309 There is no need for a "use Math::BigInt;" statement, even if you want to
3310 use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
3311 always loads it. But if you add it, add it B<before>:
3315 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
3317 Notice that the module with the last C<lib> will "win" and thus
3318 it's lib will be used if the lib is available:
3321 use Math::BigInt lib => 'Bar,Baz';
3322 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
3324 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
3325 words, Math::BigFloat will try to retain previously loaded libs when you
3326 don't specify it onem but if you specify one, it will try to load them.
3328 Actually, the lib loading order would be "Bar,Baz,Calc", and then
3329 "Foo,Bar,Baz,Calc", but independent of which lib exists, the result is the
3330 same as trying the latter load alone, except for the fact that one of Bar or
3331 Baz might be loaded needlessly in an intermidiate step (and thus hang around
3332 and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
3333 will still be tried to be loaded, but this is not as time/memory consuming as
3334 actually loading one of them. Still, this type of usage is not recommended due
3337 The old way (loading the lib only in BigInt) still works though:
3340 use Math::BigInt lib => 'Bar,Baz';
3343 You can even load Math::BigInt afterwards:
3347 use Math::BigInt lib => 'Bar,Baz';
3349 But this has the same problems like #5, it will first load Calc
3350 (Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
3351 Baz, depending on which of them works and is usable/loadable. Since this
3352 loads Calc unnec., it is not recommended.
3354 Since it also possible to just require Math::BigFloat, this poses the question
3355 about what libary this will use:
3357 require Math::BigFloat;
3358 my $x = Math::BigFloat->new(123); $x += 123;
3360 It will use Calc. Please note that the call to import() is still done, but
3361 only when you use for the first time some Math::BigFloat math (it is triggered
3362 via any constructor, so the first time you create a Math::BigFloat, the load
3363 will happen in the background). This means:
3365 require Math::BigFloat;
3366 Math::BigFloat->import ( lib => 'Foo,Bar' );
3368 would be the same as:
3370 use Math::BigFloat lib => 'Foo, Bar';
3372 But don't try to be clever to insert some operations in between:
3374 require Math::BigFloat;
3375 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
3376 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
3377 $x = Math::BigFloat->bone()+4; # now use Pari
3379 While this works, it loads Calc needlessly. But maybe you just wanted that?
3381 B<Examples #3 is highly recommended> for daily usage.
3385 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
3391 =item stringify, bstr()
3393 Both stringify and bstr() now drop the leading '+'. The old code would return
3394 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
3395 reasoning and details.
3399 The following will probably not print what you expect:
3401 print $c->bdiv(123.456),"\n";
3403 It prints both quotient and reminder since print works in list context. Also,
3404 bdiv() will modify $c, so be careful. You probably want to use
3406 print $c / 123.456,"\n";
3407 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
3413 The following will probably not print what you expect:
3415 my $c = Math::BigFloat->new('3.14159');
3416 print $c->brsft(3,10),"\n"; # prints 0.00314153.1415
3418 It prints both quotient and remainder, since print calls C<brsft()> in list
3419 context. Also, C<< $c->brsft() >> will modify $c, so be careful.
3420 You probably want to use
3422 print scalar $c->copy()->brsft(3,10),"\n";
3423 # or if you really want to modify $c
3424 print scalar $c->brsft(3,10),"\n";
3428 =item Modifying and =
3432 $x = Math::BigFloat->new(5);
3435 It will not do what you think, e.g. making a copy of $x. Instead it just makes
3436 a second reference to the B<same> object and stores it in $y. Thus anything
3437 that modifies $x will modify $y (except overloaded math operators), and vice
3438 versa. See L<Math::BigInt> for details and how to avoid that.
3442 C<bpow()> now modifies the first argument, unlike the old code which left
3443 it alone and only returned the result. This is to be consistent with
3444 C<badd()> etc. The first will modify $x, the second one won't:
3446 print bpow($x,$i),"\n"; # modify $x
3447 print $x->bpow($i),"\n"; # ditto
3448 print $x ** $i,"\n"; # leave $x alone
3450 =item precision() vs. accuracy()
3452 A common pitfall is to use L<precision()> when you want to round a result to
3453 a certain number of digits:
3457 Math::BigFloat->precision(4); # does not do what you think it does
3458 my $x = Math::BigFloat->new(12345); # rounds $x to "12000"!
3459 print "$x\n"; # print "12000"
3460 my $y = Math::BigFloat->new(3); # rounds $y to "0"!
3461 print "$y\n"; # print "0"
3462 $z = $x / $y; # 12000 / 0 => NaN!
3464 print $z->precision(),"\n"; # 4
3466 Replacing L<precision> with L<accuracy> is probably not what you want, either:
3470 Math::BigFloat->accuracy(4); # enables global rounding:
3471 my $x = Math::BigFloat->new(123456); # rounded immediately to "12350"
3472 print "$x\n"; # print "123500"
3473 my $y = Math::BigFloat->new(3); # rounded to "3
3474 print "$y\n"; # print "3"
3475 print $z = $x->copy()->bdiv($y),"\n"; # 41170
3476 print $z->accuracy(),"\n"; # 4
3478 What you want to use instead is:
3482 my $x = Math::BigFloat->new(123456); # no rounding
3483 print "$x\n"; # print "123456"
3484 my $y = Math::BigFloat->new(3); # no rounding
3485 print "$y\n"; # print "3"
3486 print $z = $x->copy()->bdiv($y,4),"\n"; # 41150
3487 print $z->accuracy(),"\n"; # undef
3489 In addition to computing what you expected, the last example also does B<not>
3490 "taint" the result with an accuracy or precision setting, which would
3491 influence any further operation.
3497 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
3498 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
3500 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
3501 because they solve the autoupgrading/downgrading issue, at least partly.
3503 The package at L<http://search.cpan.org/~tels/Math-BigInt> contains
3504 more documentation including a full version history, testcases, empty
3505 subclass files and benchmarks.
3509 This program is free software; you may redistribute it and/or modify it under
3510 the same terms as Perl itself.
3514 Mark Biggar, overloaded interface by Ilya Zakharevich.
3515 Completely rewritten by Tels L<http://bloodgate.com> in 2001 - 2006, and still