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/Math::BigInt/;
23 # $_trap_inf/$_trap_nan are internal and should never be accessed from outside
24 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode
25 $upgrade $downgrade $_trap_nan $_trap_inf/;
26 my $class = "Math::BigFloat";
29 '<=>' => sub { my $rc = $_[2] ?
30 ref($_[0])->bcmp($_[1],$_[0]) :
31 ref($_[0])->bcmp($_[0],$_[1]);
32 $rc = 1 unless defined $rc;
35 # we need '>=' to get things like "1 >= NaN" right:
36 '>=' => sub { my $rc = $_[2] ?
37 ref($_[0])->bcmp($_[1],$_[0]) :
38 ref($_[0])->bcmp($_[0],$_[1]);
39 # if there was a NaN involved, return false
40 return '' unless defined $rc;
43 'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
46 ##############################################################################
47 # global constants, flags and assorted stuff
49 # the following are public, but their usage is not recommended. Use the
50 # accessor methods instead.
52 # class constants, use Class->constant_name() to access
53 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
60 # the package we are using for our private parts, defaults to:
61 # Math::BigInt->config()->{lib}
62 my $MBI = 'Math::BigInt::FastCalc';
64 # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
66 # the same for infinity
69 # constant for easier life
72 my $IMPORT = 0; # was import() called yet? used to make require work
74 # some digits of accuracy for blog(undef,10); which we use in blog() for speed
76 '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
77 my $LOG_10_A = length($LOG_10)-1;
80 '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
81 my $LOG_2_A = length($LOG_2)-1;
82 my $HALF = '0.5'; # made into an object if nec.
84 ##############################################################################
85 # the old code had $rnd_mode, so we need to support it, too
87 sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
88 sub FETCH { return $round_mode; }
89 sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
93 # when someone sets $rnd_mode, we catch this and check the value to see
94 # whether it is valid or not.
95 $rnd_mode = 'even'; tie $rnd_mode, 'Math::BigFloat';
97 # we need both of them in this package:
98 *as_int = \&as_number;
101 ##############################################################################
104 # valid method aliases for AUTOLOAD
105 my %methods = map { $_ => 1 }
106 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
107 fint facmp fcmp fzero fnan finf finc fdec ffac fneg
108 fceil ffloor frsft flsft fone flog froot fexp
110 # valid methods that can be handed up (for AUTOLOAD)
111 my %hand_ups = map { $_ => 1 }
112 qw / is_nan is_inf is_negative is_positive is_pos is_neg
113 accuracy precision div_scale round_mode fabs fnot
114 objectify upgrade downgrade
119 sub _method_alias { exists $methods{$_[0]||''}; }
120 sub _method_hand_up { exists $hand_ups{$_[0]||''}; }
123 ##############################################################################
128 # create a new BigFloat object from a string or another bigfloat object.
131 # sign => sign (+/-), or "NaN"
133 my ($class,$wanted,@r) = @_;
135 # avoid numify-calls by not using || on $wanted!
136 return $class->bzero() if !defined $wanted; # default to 0
137 return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
139 $class->import() if $IMPORT == 0; # make require work
141 my $self = {}; bless $self, $class;
142 # shortcut for bigints and its subclasses
143 if ((ref($wanted)) && UNIVERSAL::can( $wanted, "as_number"))
145 $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy
146 $self->{_e} = $MBI->_zero();
148 $self->{sign} = $wanted->sign();
149 return $self->bnorm();
151 # else: got a string or something maskerading as number (with overload)
153 # handle '+inf', '-inf' first
154 if ($wanted =~ /^[+-]?inf\z/)
156 return $downgrade->new($wanted) if $downgrade;
158 $self->{sign} = $wanted; # set a default sign for bstr()
159 return $self->binf($wanted);
162 # shortcut for simple forms like '12' that neither have trailing nor leading
164 if ($wanted =~ /^([+-]?)([1-9][0-9]*[1-9])$/)
166 $self->{_e} = $MBI->_zero();
168 $self->{sign} = $1 || '+';
169 $self->{_m} = $MBI->_new($2);
170 return $self->round(@r) if !$downgrade;
173 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
179 Carp::croak ("$wanted is not a number initialized to $class");
182 return $downgrade->bnan() if $downgrade;
184 $self->{_e} = $MBI->_zero();
186 $self->{_m} = $MBI->_zero();
187 $self->{sign} = $nan;
191 # make integer from mantissa by adjusting exp, then convert to int
192 $self->{_e} = $MBI->_new($$ev); # exponent
193 $self->{_es} = $$es || '+';
194 my $mantissa = "$$miv$$mfv"; # create mant.
195 $mantissa =~ s/^0+(\d)/$1/; # strip leading zeros
196 $self->{_m} = $MBI->_new($mantissa); # create mant.
198 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
199 if (CORE::length($$mfv) != 0)
201 my $len = $MBI->_new( CORE::length($$mfv));
202 ($self->{_e}, $self->{_es}) =
203 _e_sub ($self->{_e}, $len, $self->{_es}, '+');
205 # we can only have trailing zeros on the mantissa if $$mfv eq ''
208 # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
209 # because that is faster, especially when _m is not stored in base 10.
210 my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/;
213 my $z = $MBI->_new($zeros);
214 # turn '120e2' into '12e3'
215 $MBI->_rsft ( $self->{_m}, $z, 10);
216 ($self->{_e}, $self->{_es}) =
217 _e_add ( $self->{_e}, $z, $self->{_es}, '+');
220 $self->{sign} = $$mis;
222 # for something like 0Ey, set y to 1, and -0 => +0
223 # Check $$miv for being '0' and $$mfv eq '', because otherwise _m could not
224 # have become 0. That's faster than to call $MBI->_is_zero().
225 $self->{sign} = '+', $self->{_e} = $MBI->_one()
226 if $$miv eq '0' and $$mfv eq '';
228 return $self->round(@r) if !$downgrade;
230 # if downgrade, inf, NaN or integers go down
232 if ($downgrade && $self->{_es} eq '+')
234 if ($MBI->_is_zero( $self->{_e} ))
236 return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
238 return $downgrade->new($self->bsstr());
240 $self->bnorm()->round(@r); # first normalize, then round
245 # if two arguments, the first one is the class to "swallow" subclasses
249 sign => $_[1]->{sign},
251 _m => $MBI->_copy($_[1]->{_m}),
252 _e => $MBI->_copy($_[1]->{_e}),
255 $self->{_a} = $_[1]->{_a} if defined $_[1]->{_a};
256 $self->{_p} = $_[1]->{_p} if defined $_[1]->{_p};
261 sign => $_[0]->{sign},
263 _m => $MBI->_copy($_[0]->{_m}),
264 _e => $MBI->_copy($_[0]->{_e}),
267 $self->{_a} = $_[0]->{_a} if defined $_[0]->{_a};
268 $self->{_p} = $_[0]->{_p} if defined $_[0]->{_p};
274 # used by parent class bone() to initialize number to NaN
280 my $class = ref($self);
281 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
284 $IMPORT=1; # call our import only once
285 $self->{_m} = $MBI->_zero();
286 $self->{_e} = $MBI->_zero();
292 # used by parent class bone() to initialize number to +-inf
298 my $class = ref($self);
299 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
302 $IMPORT=1; # call our import only once
303 $self->{_m} = $MBI->_zero();
304 $self->{_e} = $MBI->_zero();
310 # used by parent class bone() to initialize number to 1
312 $IMPORT=1; # call our import only once
313 $self->{_m} = $MBI->_one();
314 $self->{_e} = $MBI->_zero();
320 # used by parent class bone() to initialize number to 0
322 $IMPORT=1; # call our import only once
323 $self->{_m} = $MBI->_zero();
324 $self->{_e} = $MBI->_one();
330 my ($self,$class) = @_;
331 return if $class =~ /^Math::BigInt/; # we aren't one of these
332 UNIVERSAL::isa($self,$class);
337 # return (later set?) configuration data as hash ref
338 my $class = shift || 'Math::BigFloat';
340 if (@_ == 1 && ref($_[0]) ne 'HASH')
342 my $cfg = $class->SUPER::config();
343 return $cfg->{$_[0]};
346 my $cfg = $class->SUPER::config(@_);
348 # now we need only to override the ones that are different from our parent
349 $cfg->{class} = $class;
354 ##############################################################################
355 # string conversation
359 # (ref to BFLOAT or num_str ) return num_str
360 # Convert number from internal format to (non-scientific) string format.
361 # internal format is always normalized (no leading zeros, "-0" => "+0")
362 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
364 if ($x->{sign} !~ /^[+-]$/)
366 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
370 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
373 my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
376 $es = $MBI->_str($x->{_m});
377 $len = CORE::length($es);
378 my $e = $MBI->_num($x->{_e});
379 $e = -$e if $x->{_es} eq '-';
383 # if _e is bigger than a scalar, the following will blow your memory
386 my $r = abs($e) - $len;
387 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
391 substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
392 $cad = -$cad if $x->{_es} eq '-';
398 $es .= '0' x $e; $len += $e; $cad = 0;
402 $es = '-'.$es if $x->{sign} eq '-';
403 # if set accuracy or precision, pad with zeros on the right side
404 if ((defined $x->{_a}) && ($not_zero))
406 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
407 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
408 $zeros = $x->{_a} - $len if $cad != $len;
409 $es .= $dot.'0' x $zeros if $zeros > 0;
411 elsif ((($x->{_p} || 0) < 0))
413 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
414 my $zeros = -$x->{_p} + $cad;
415 $es .= $dot.'0' x $zeros if $zeros > 0;
422 # (ref to BFLOAT or num_str ) return num_str
423 # Convert number from internal format to scientific string format.
424 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
425 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
427 if ($x->{sign} !~ /^[+-]$/)
429 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
432 my $sep = 'e'.$x->{_es};
433 my $sign = $x->{sign}; $sign = '' if $sign eq '+';
434 $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
439 # Make a number from a BigFloat object
440 # simple return a string and let Perl's atoi()/atof() handle the rest
441 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
445 ##############################################################################
446 # public stuff (usually prefixed with "b")
450 # (BINT or num_str) return BINT
451 # negate number or make a negated number from string
452 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
454 return $x if $x->modify('bneg');
456 # for +0 dont negate (to have always normalized +0). Does nothing for 'NaN'
457 $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
462 # XXX TODO this must be overwritten and return NaN for non-integer values
463 # band(), bior(), bxor(), too
466 # $class->SUPER::bnot($class,@_);
471 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
474 my ($self,$x,$y) = (ref($_[0]),@_);
475 # objectify is costly, so avoid it
476 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
478 ($self,$x,$y) = objectify(2,@_);
481 return $upgrade->bcmp($x,$y) if defined $upgrade &&
482 ((!$x->isa($self)) || (!$y->isa($self)));
484 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
486 # handle +-inf and NaN
487 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
488 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
489 return +1 if $x->{sign} eq '+inf';
490 return -1 if $x->{sign} eq '-inf';
491 return -1 if $y->{sign} eq '+inf';
495 # check sign for speed first
496 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
497 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
500 my $xz = $x->is_zero();
501 my $yz = $y->is_zero();
502 return 0 if $xz && $yz; # 0 <=> 0
503 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
504 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
506 # adjust so that exponents are equal
507 my $lxm = $MBI->_len($x->{_m});
508 my $lym = $MBI->_len($y->{_m});
509 # the numify somewhat limits our length, but makes it much faster
510 my ($xes,$yes) = (1,1);
511 $xes = -1 if $x->{_es} ne '+';
512 $yes = -1 if $y->{_es} ne '+';
513 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
514 my $ly = $lym + $yes * $MBI->_num($y->{_e});
515 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
516 return $l <=> 0 if $l != 0;
518 # lengths (corrected by exponent) are equal
519 # so make mantissa equal length by padding with zero (shift left)
520 my $diff = $lxm - $lym;
521 my $xm = $x->{_m}; # not yet copy it
525 $ym = $MBI->_copy($y->{_m});
526 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
530 $xm = $MBI->_copy($x->{_m});
531 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
533 my $rc = $MBI->_acmp($xm,$ym);
534 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
540 # Compares 2 values, ignoring their signs.
541 # Returns one of undef, <0, =0, >0. (suitable for sort)
544 my ($self,$x,$y) = (ref($_[0]),@_);
545 # objectify is costly, so avoid it
546 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
548 ($self,$x,$y) = objectify(2,@_);
551 return $upgrade->bacmp($x,$y) if defined $upgrade &&
552 ((!$x->isa($self)) || (!$y->isa($self)));
554 # handle +-inf and NaN's
555 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
557 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
558 return 0 if ($x->is_inf() && $y->is_inf());
559 return 1 if ($x->is_inf() && !$y->is_inf());
564 my $xz = $x->is_zero();
565 my $yz = $y->is_zero();
566 return 0 if $xz && $yz; # 0 <=> 0
567 return -1 if $xz && !$yz; # 0 <=> +y
568 return 1 if $yz && !$xz; # +x <=> 0
570 # adjust so that exponents are equal
571 my $lxm = $MBI->_len($x->{_m});
572 my $lym = $MBI->_len($y->{_m});
573 my ($xes,$yes) = (1,1);
574 $xes = -1 if $x->{_es} ne '+';
575 $yes = -1 if $y->{_es} ne '+';
576 # the numify somewhat limits our length, but makes it much faster
577 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
578 my $ly = $lym + $yes * $MBI->_num($y->{_e});
580 return $l <=> 0 if $l != 0;
582 # lengths (corrected by exponent) are equal
583 # so make mantissa equal-length by padding with zero (shift left)
584 my $diff = $lxm - $lym;
585 my $xm = $x->{_m}; # not yet copy it
589 $ym = $MBI->_copy($y->{_m});
590 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
594 $xm = $MBI->_copy($x->{_m});
595 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
597 $MBI->_acmp($xm,$ym);
602 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
603 # return result as BFLOAT
606 my ($self,$x,$y,@r) = (ref($_[0]),@_);
607 # objectify is costly, so avoid it
608 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
610 ($self,$x,$y,@r) = objectify(2,@_);
613 return $x if $x->modify('badd');
615 # inf and NaN handling
616 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
619 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
621 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
623 # +inf++inf or -inf+-inf => same, rest is NaN
624 return $x if $x->{sign} eq $y->{sign};
627 # +-inf + something => +inf; something +-inf => +-inf
628 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
632 return $upgrade->badd($x,$y,@r) if defined $upgrade &&
633 ((!$x->isa($self)) || (!$y->isa($self)));
635 $r[3] = $y; # no push!
637 # speed: no add for 0+y or x+0
638 return $x->bround(@r) if $y->is_zero(); # x+0
639 if ($x->is_zero()) # 0+y
641 # make copy, clobbering up x (modify in place!)
642 $x->{_e} = $MBI->_copy($y->{_e});
643 $x->{_es} = $y->{_es};
644 $x->{_m} = $MBI->_copy($y->{_m});
645 $x->{sign} = $y->{sign} || $nan;
646 return $x->round(@r);
649 # take lower of the two e's and adapt m1 to it to match m2
651 $e = $MBI->_zero() if !defined $e; # if no BFLOAT?
652 $e = $MBI->_copy($e); # make copy (didn't do it yet)
656 ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
658 my $add = $MBI->_copy($y->{_m});
660 if ($es eq '-') # < 0
662 $MBI->_lsft( $x->{_m}, $e, 10);
663 ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
665 elsif (!$MBI->_is_zero($e)) # > 0
667 $MBI->_lsft($add, $e, 10);
669 # else: both e are the same, so just leave them
671 if ($x->{sign} eq $y->{sign})
674 $x->{_m} = $MBI->_add($x->{_m}, $add);
678 ($x->{_m}, $x->{sign}) =
679 _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
682 # delete trailing zeros, then round
683 $x->bnorm()->round(@r);
686 # sub bsub is inherited from Math::BigInt!
690 # increment arg by one
691 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
693 return $x if $x->modify('binc');
695 if ($x->{_es} eq '-')
697 return $x->badd($self->bone(),@r); # digits after dot
700 if (!$MBI->_is_zero($x->{_e})) # _e == 0 for NaN, inf, -inf
702 # 1e2 => 100, so after the shift below _m has a '0' as last digit
703 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
704 $x->{_e} = $MBI->_zero(); # normalize
706 # we know that the last digit of $x will be '1' or '9', depending on the
710 if ($x->{sign} eq '+')
712 $MBI->_inc($x->{_m});
713 return $x->bnorm()->bround(@r);
715 elsif ($x->{sign} eq '-')
717 $MBI->_dec($x->{_m});
718 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
719 return $x->bnorm()->bround(@r);
721 # inf, nan handling etc
722 $x->badd($self->bone(),@r); # badd() does round
727 # decrement arg by one
728 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
730 return $x if $x->modify('bdec');
732 if ($x->{_es} eq '-')
734 return $x->badd($self->bone('-'),@r); # digits after dot
737 if (!$MBI->_is_zero($x->{_e}))
739 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
740 $x->{_e} = $MBI->_zero(); # normalize
744 my $zero = $x->is_zero();
746 if (($x->{sign} eq '-') || $zero)
748 $MBI->_inc($x->{_m});
749 $x->{sign} = '-' if $zero; # 0 => 1 => -1
750 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
751 return $x->bnorm()->round(@r);
754 elsif ($x->{sign} eq '+')
756 $MBI->_dec($x->{_m});
757 return $x->bnorm()->round(@r);
759 # inf, nan handling etc
760 $x->badd($self->bone('-'),@r); # does round
767 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
769 return $x if $x->modify('blog');
771 # $base > 0, $base != 1; if $base == undef default to $base == e
774 # we need to limit the accuracy to protect against overflow
777 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
779 # also takes care of the "error in _find_round_parameters?" case
780 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
782 # no rounding at all, so must use fallback
783 if (scalar @params == 0)
785 # simulate old behaviour
786 $params[0] = $self->div_scale(); # and round to it as accuracy
787 $params[1] = undef; # P = undef
788 $scale = $params[0]+4; # at least four more for proper round
789 $params[2] = $r; # round mode by caller or undef
790 $fallback = 1; # to clear a/p afterwards
794 # the 4 below is empirical, and there might be cases where it is not
796 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
799 return $x->bzero(@params) if $x->is_one();
800 # base not defined => base == Euler's number e
803 # make object, since we don't feed it through objectify() to still get the
804 # case of $base == undef
805 $base = $self->new($base) unless ref($base);
806 # $base > 0; $base != 1
807 return $x->bnan() if $base->is_zero() || $base->is_one() ||
808 $base->{sign} ne '+';
809 # if $x == $base, we know the result must be 1.0
810 if ($x->bcmp($base) == 0)
812 $x->bone('+',@params);
815 # clear a/p after round, since user did not request it
816 delete $x->{_a}; delete $x->{_p};
822 # when user set globals, they would interfere with our calculation, so
823 # disable them and later re-enable them
825 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
826 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
827 # we also need to disable any set A or P on $x (_find_round_parameters took
828 # them already into account), since these would interfere, too
829 delete $x->{_a}; delete $x->{_p};
830 # need to disable $upgrade in BigInt, to avoid deep recursion
831 local $Math::BigInt::upgrade = undef;
832 local $Math::BigFloat::downgrade = undef;
834 # upgrade $x if $x is not a BigFloat (handle BigInt input)
836 if (!$x->isa('Math::BigFloat'))
838 $x = Math::BigFloat->new($x);
844 # If the base is defined and an integer, try to calculate integer result
845 # first. This is very fast, and in case the real result was found, we can
847 if (defined $base && $base->is_int() && $x->is_int())
849 my $i = $MBI->_copy( $x->{_m} );
850 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
851 my $int = Math::BigInt->bzero();
853 $int->blog($base->as_number());
855 if ($base->as_number()->bpow($int) == $x)
857 # found result, return it
858 $x->{_m} = $int->{value};
859 $x->{_e} = $MBI->_zero();
868 # base is undef, so base should be e (Euler's number), so first calculate the
869 # log to base e (using reduction by 10 (and probably 2)):
870 $self->_log_10($x,$scale);
872 # and if a different base was requested, convert it
875 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
876 # not ln, but some other base (don't modify $base)
877 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
881 # shortcut to not run through _find_round_parameters again
882 if (defined $params[0])
884 $x->bround($params[0],$params[2]); # then round accordingly
888 $x->bfround($params[1],$params[2]); # then round accordingly
892 # clear a/p after round, since user did not request it
893 delete $x->{_a}; delete $x->{_p};
896 $$abr = $ab; $$pbr = $pb;
903 # Given D (digits in decimal), compute N so that N! (N factorial) is
904 # at least D digits long. D should be at least 50.
907 # two constants for the Ramanujan estimate of ln(N!)
908 my $lg2 = log(2 * 3.14159265) / 2;
911 # D = 50 => N => 42, so L = 40 and R = 50
912 my $l = 40; my $r = $d;
914 # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
915 $l = $l->numify if ref($l);
916 $r = $r->numify if ref($r);
917 $lg2 = $lg2->numify if ref($lg2);
918 $lg10 = $lg10->numify if ref($lg10);
920 # binary search for the right value (could this be written as the reverse of lg(n!)?)
923 my $n = int(($r - $l) / 2) + $l;
925 int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
926 $ramanujan > $d ? $r = $n : $l = $n;
933 # Calculate n over k (binomial coefficient or "choose" function) as integer.
935 my ($self,$x,$y,@r) = (ref($_[0]),@_);
937 # objectify is costly, so avoid it
938 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
940 ($self,$x,$y,@r) = objectify(2,@_);
943 return $x if $x->modify('bnok');
945 return $x->bnan() if $x->is_nan() || $y->is_nan();
946 return $x->binf() if $x->is_inf();
948 my $u = $x->as_int();
949 $u->bnok($y->as_int());
951 $x->{_m} = $u->{value};
952 $x->{_e} = $MBI->_zero();
960 # Calculate e ** X (Euler's number to the power of X)
961 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
963 return $x if $x->modify('bexp');
965 return $x->binf() if $x->{sign} eq '+inf';
966 return $x->bzero() if $x->{sign} eq '-inf';
968 # we need to limit the accuracy to protect against overflow
971 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
973 # also takes care of the "error in _find_round_parameters?" case
974 return $x if $x->{sign} eq 'NaN';
976 # no rounding at all, so must use fallback
977 if (scalar @params == 0)
979 # simulate old behaviour
980 $params[0] = $self->div_scale(); # and round to it as accuracy
981 $params[1] = undef; # P = undef
982 $scale = $params[0]+4; # at least four more for proper round
983 $params[2] = $r; # round mode by caller or undef
984 $fallback = 1; # to clear a/p afterwards
988 # the 4 below is empirical, and there might be cases where it's not enough...
989 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
992 return $x->bone(@params) if $x->is_zero();
994 if (!$x->isa('Math::BigFloat'))
996 $x = Math::BigFloat->new($x);
1000 # when user set globals, they would interfere with our calculation, so
1001 # disable them and later re-enable them
1003 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1004 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1005 # we also need to disable any set A or P on $x (_find_round_parameters took
1006 # them already into account), since these would interfere, too
1007 delete $x->{_a}; delete $x->{_p};
1008 # need to disable $upgrade in BigInt, to avoid deep recursion
1009 local $Math::BigInt::upgrade = undef;
1010 local $Math::BigFloat::downgrade = undef;
1012 my $x_org = $x->copy();
1014 # We use the following Taylor series:
1017 # e = 1 + --- + --- + --- + --- ...
1020 # The difference for each term is X and N, which would result in:
1021 # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term
1023 # But it is faster to compute exp(1) and then raising it to the
1024 # given power, esp. if $x is really big and an integer because:
1026 # * The numerator is always 1, making the computation faster
1027 # * the series converges faster in the case of x == 1
1028 # * We can also easily check when we have reached our limit: when the
1029 # term to be added is smaller than "1E$scale", we can stop - f.i.
1030 # scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
1031 # * we can compute the *exact* result by simulating bigrat math:
1033 # 1 1 gcd(3,4) = 1 1*24 + 1*6 5
1034 # - + - = ---------- = --
1037 # We do not compute the gcd() here, but simple do:
1039 # - + - = --------- = --
1043 # a c a*d + c*b and note that c is always 1 and d = (b*f)
1047 # This leads to: which can be reduced by b to:
1048 # a 1 a*b*f + b a*f + 1
1049 # - + - = --------- = -------
1052 # The first terms in the series are:
1054 # 1 1 1 1 1 1 1 1 13700
1055 # -- + -- + -- + -- + -- + --- + --- + ---- = -----
1056 # 1 1 2 6 24 120 720 5040 5040
1058 # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!
1062 # set $x directly from a cached string form
1063 $x->{_m} = $MBI->_new(
1064 "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
1067 $x->{_e} = $MBI->_new(79);
1071 # compute A and B so that e = A / B.
1073 # After some terms we end up with this, so we use it as a starting point:
1074 my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
1075 my $F = $MBI->_new(42); my $step = 42;
1077 # Compute how many steps we need to take to get $A and $B sufficiently big
1078 my $steps = _len_to_steps($scale - 4);
1079 # print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
1080 while ($step++ <= $steps)
1082 # calculate $a * $f + 1
1083 $A = $MBI->_mul($A, $F);
1084 $A = $MBI->_inc($A);
1086 $F = $MBI->_inc($F);
1088 # compute $B as factorial of $steps (this is faster than doing it manually)
1089 my $B = $MBI->_fac($MBI->_new($steps));
1091 # print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
1093 # compute A/B with $scale digits in the result (truncate, not round)
1094 $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
1095 $A = $MBI->_div( $A, $B );
1100 $x->{_e} = $MBI->_new($scale);
1103 # $x contains now an estimate of e, with some surplus digits, so we can round
1104 if (!$x_org->is_one())
1106 # raise $x to the wanted power and round it in one step:
1107 $x->bpow($x_org, @params);
1111 # else just round the already computed result
1112 delete $x->{_a}; delete $x->{_p};
1113 # shortcut to not run through _find_round_parameters again
1114 if (defined $params[0])
1116 $x->bround($params[0],$params[2]); # then round accordingly
1120 $x->bfround($params[1],$params[2]); # then round accordingly
1125 # clear a/p after round, since user did not request it
1126 delete $x->{_a}; delete $x->{_p};
1129 $$abr = $ab; $$pbr = $pb;
1131 $x; # return modified $x
1136 # internal log function to calculate ln() based on Taylor series.
1137 # Modifies $x in place.
1138 my ($self,$x,$scale) = @_;
1140 # in case of $x == 1, result is 0
1141 return $x->bzero() if $x->is_one();
1143 # XXX TODO: rewrite this in a similiar manner to bexp()
1145 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
1149 # Taylor: | u 1 u^3 1 u^5 |
1150 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
1151 # |_ v 3 v^3 5 v^5 _|
1153 # This takes much more steps to calculate the result and is thus not used
1156 # Taylor: | u 1 u^2 1 u^3 |
1157 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
1158 # |_ x 2 x^2 3 x^3 _|
1160 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
1162 $v = $x->copy(); $v->binc(); # v = x+1
1163 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
1164 $x->bdiv($v,$scale); # first term: u/v
1165 $below = $v->copy();
1167 $u *= $u; $v *= $v; # u^2, v^2
1168 $below->bmul($v); # u^3, v^3
1170 $factor = $self->new(3); $f = $self->new(2);
1172 my $steps = 0 if DEBUG;
1173 $limit = $self->new("1E-". ($scale-1));
1176 # we calculate the next term, and add it to the last
1177 # when the next term is below our limit, it won't affect the outcome
1178 # anymore, so we stop
1180 # calculating the next term simple from over/below will result in quite
1181 # a time hog if the input has many digits, since over and below will
1182 # accumulate more and more digits, and the result will also have many
1183 # digits, but in the end it is rounded to $scale digits anyway. So if we
1184 # round $over and $below first, we save a lot of time for the division
1185 # (not with log(1.2345), but try log (123**123) to see what I mean. This
1186 # can introduce a rounding error if the division result would be f.i.
1187 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
1188 # if we truncated $over and $below we might get 0.12345. Does this matter
1189 # for the end result? So we give $over and $below 4 more digits to be
1190 # on the safe side (unscientific error handling as usual... :+D
1192 $next = $over->copy->bround($scale+4)->bdiv(
1193 $below->copy->bmul($factor)->bround($scale+4),
1197 ## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
1199 last if $next->bacmp($limit) <= 0;
1201 delete $next->{_a}; delete $next->{_p};
1203 # calculate things for the next term
1204 $over *= $u; $below *= $v; $factor->badd($f);
1207 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
1210 print "took $steps steps\n" if DEBUG;
1211 $x->bmul($f); # $x *= 2
1216 # Internal log function based on reducing input to the range of 0.1 .. 9.99
1217 # and then "correcting" the result to the proper one. Modifies $x in place.
1218 my ($self,$x,$scale) = @_;
1220 # Taking blog() from numbers greater than 10 takes a *very long* time, so we
1221 # break the computation down into parts based on the observation that:
1222 # blog(X*Y) = blog(X) + blog(Y)
1223 # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
1224 # $x is the faster it gets. Since 2*$x takes about 10 times as
1225 # long, we make it faster by about a factor of 100 by dividing $x by 10.
1227 # The same observation is valid for numbers smaller than 0.1, e.g. computing
1228 # log(1) is fastest, and the further away we get from 1, the longer it takes.
1229 # So we also 'break' this down by multiplying $x with 10 and subtract the
1230 # log(10) afterwards to get the correct result.
1232 # To get $x even closer to 1, we also divide by 2 and then use log(2) to
1233 # correct for this. For instance if $x is 2.4, we use the formula:
1234 # blog(2.4 * 2) == blog (1.2) + blog(2)
1235 # and thus calculate only blog(1.2) and blog(2), which is faster in total
1236 # than calculating blog(2.4).
1238 # In addition, the values for blog(2) and blog(10) are cached.
1240 # Calculate nr of digits before dot:
1241 my $dbd = $MBI->_num($x->{_e});
1242 $dbd = -$dbd if $x->{_es} eq '-';
1243 $dbd += $MBI->_len($x->{_m});
1245 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
1246 # infinite recursion
1248 my $calc = 1; # do some calculation?
1250 # disable the shortcut for 10, since we need log(10) and this would recurse
1252 if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
1254 $dbd = 0; # disable shortcut
1255 # we can use the cached value in these cases
1256 if ($scale <= $LOG_10_A)
1258 $x->bzero(); $x->badd($LOG_10); # modify $x in place
1259 $calc = 0; # no need to calc, but round
1261 # if we can't use the shortcut, we continue normally
1265 # disable the shortcut for 2, since we maybe have it cached
1266 if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
1268 $dbd = 0; # disable shortcut
1269 # we can use the cached value in these cases
1270 if ($scale <= $LOG_2_A)
1272 $x->bzero(); $x->badd($LOG_2); # modify $x in place
1273 $calc = 0; # no need to calc, but round
1275 # if we can't use the shortcut, we continue normally
1279 # if $x = 0.1, we know the result must be 0-log(10)
1280 if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
1281 $MBI->_is_one($x->{_m}))
1283 $dbd = 0; # disable shortcut
1284 # we can use the cached value in these cases
1285 if ($scale <= $LOG_10_A)
1287 $x->bzero(); $x->bsub($LOG_10);
1288 $calc = 0; # no need to calc, but round
1292 return if $calc == 0; # already have the result
1294 # default: these correction factors are undef and thus not used
1295 my $l_10; # value of ln(10) to A of $scale
1296 my $l_2; # value of ln(2) to A of $scale
1298 my $two = $self->new(2);
1300 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1301 # so don't do this shortcut for 1 or 0
1302 if (($dbd > 1) || ($dbd < 0))
1304 # convert our cached value to an object if not already (avoid doing this
1305 # at import() time, since not everybody needs this)
1306 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1308 #print "x = $x, dbd = $dbd, calc = $calc\n";
1309 # got more than one digit before the dot, or more than one zero after the
1311 # log(123) == log(1.23) + log(10) * 2
1312 # log(0.0123) == log(1.23) - log(10) * 2
1314 if ($scale <= $LOG_10_A)
1317 $l_10 = $LOG_10->copy(); # copy for mul
1321 # else: slower, compute and cache result
1322 # also disable downgrade for this code path
1323 local $Math::BigFloat::downgrade = undef;
1325 # shorten the time to calculate log(10) based on the following:
1326 # log(1.25 * 8) = log(1.25) + log(8)
1327 # = log(1.25) + log(2) + log(2) + log(2)
1329 # first get $l_2 (and possible compute and cache log(2))
1330 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1331 if ($scale <= $LOG_2_A)
1334 $l_2 = $LOG_2->copy(); # copy() for the mul below
1338 # else: slower, compute and cache result
1339 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1340 $LOG_2 = $l_2->copy(); # cache the result for later
1341 # the copy() is for mul below
1345 # now calculate log(1.25):
1346 $l_10 = $self->new('1.25'); $self->_log($l_10, $scale); # scale+4, actually
1348 # log(1.25) + log(2) + log(2) + log(2):
1352 $LOG_10 = $l_10->copy(); # cache the result for later
1353 # the copy() is for mul below
1356 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
1357 $l_10->bmul( $self->new($dbd)); # log(10) * (digits_before_dot-1)
1364 ($x->{_e}, $x->{_es}) =
1365 _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
1369 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1371 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1372 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1374 $HALF = $self->new($HALF) unless ref($HALF);
1376 my $twos = 0; # default: none (0 times)
1377 while ($x->bacmp($HALF) <= 0) # X <= 0.5
1379 $twos--; $x->bmul($two);
1381 while ($x->bacmp($two) >= 0) # X >= 2
1383 $twos++; $x->bdiv($two,$scale+4); # keep all digits
1385 # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
1386 # So calculate correction factor based on ln(2):
1389 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1390 if ($scale <= $LOG_2_A)
1393 $l_2 = $LOG_2->copy(); # copy() for the mul below
1397 # else: slower, compute and cache result
1398 # also disable downgrade for this code path
1399 local $Math::BigFloat::downgrade = undef;
1400 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1401 $LOG_2 = $l_2->copy(); # cache the result for later
1402 # the copy() is for mul below
1405 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
1408 $self->_log($x,$scale); # need to do the "normal" way
1409 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
1410 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
1412 # all done, $x contains now the result
1418 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1419 # does not modify arguments, but returns new object
1420 # Lowest Common Multiplicator
1422 my ($self,@arg) = objectify(0,@_);
1423 my $x = $self->new(shift @arg);
1424 while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); }
1430 # (BINT or num_str, BINT or num_str) return BINT
1431 # does not modify arguments, but returns new object
1434 $y = __PACKAGE__->new($y) if !ref($y);
1436 my $x = $y->copy()->babs(); # keep arguments
1438 return $x->bnan() if $x->{sign} !~ /^[+-]$/ # x NaN?
1439 || !$x->is_int(); # only for integers now
1443 my $t = shift; $t = $self->new($t) if !ref($t);
1444 $y = $t->copy()->babs();
1446 return $x->bnan() if $y->{sign} !~ /^[+-]$/ # y NaN?
1447 || !$y->is_int(); # only for integers now
1449 # greatest common divisor
1450 while (! $y->is_zero())
1452 ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1455 last if $x->is_one();
1460 ##############################################################################
1464 # Internal helper sub to take two positive integers and their signs and
1465 # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1466 # output ($CALC,('+'|'-'))
1467 my ($x,$y,$xs,$ys) = @_;
1469 # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1472 $x = $MBI->_add ($x, $y ); # a+b
1473 # the sign follows $xs
1477 my $a = $MBI->_acmp($x,$y);
1480 $x = $MBI->_sub ($x , $y); # abs sub
1484 $x = $MBI->_zero(); # result is 0
1489 $x = $MBI->_sub ( $y, $x, 1 ); # abs sub
1497 # Internal helper sub to take two positive integers and their signs and
1498 # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1499 # output ($CALC,('+'|'-'))
1500 my ($x,$y,$xs,$ys) = @_;
1504 _e_add($x,$y,$xs,$ys); # call add (does subtract now)
1507 ###############################################################################
1508 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1512 # return true if arg (BFLOAT or num_str) is an integer
1513 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1515 (($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1516 ($x->{_es} eq '+')) ? 1 : 0; # 1e-1 => no integer
1521 # return true if arg (BFLOAT or num_str) is zero
1522 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1524 ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m})) ? 1 : 0;
1529 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1530 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1532 $sign = '+' if !defined $sign || $sign ne '-';
1534 ($x->{sign} eq $sign &&
1535 $MBI->_is_zero($x->{_e}) &&
1536 $MBI->_is_one($x->{_m}) ) ? 1 : 0;
1541 # return true if arg (BFLOAT or num_str) is odd or false if even
1542 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1544 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1545 ($MBI->_is_zero($x->{_e})) &&
1546 ($MBI->_is_odd($x->{_m}))) ? 1 : 0;
1551 # return true if arg (BINT or num_str) is even or false if odd
1552 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1554 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1555 ($x->{_es} eq '+') && # 123.45 isn't
1556 ($MBI->_is_even($x->{_m}))) ? 1 : 0; # but 1200 is
1561 # multiply two numbers
1564 my ($self,$x,$y,@r) = (ref($_[0]),@_);
1565 # objectify is costly, so avoid it
1566 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1568 ($self,$x,$y,@r) = objectify(2,@_);
1571 return $x if $x->modify('bmul');
1573 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1576 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1578 return $x->bnan() if $x->is_zero() || $y->is_zero();
1579 # result will always be +-inf:
1580 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1581 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1582 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1583 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1584 return $x->binf('-');
1587 return $upgrade->bmul($x,$y,@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});
1594 $r[3] = $y; # no push!
1597 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1598 $x->bnorm->round(@r);
1603 # multiply two numbers and add the third to the result
1606 my ($self,$x,$y,$z,@r) = (ref($_[0]),@_);
1607 # objectify is costly, so avoid it
1608 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1610 ($self,$x,$y,$z,@r) = objectify(3,@_);
1613 return $x if $x->modify('bmuladd');
1615 return $x->bnan() if (($x->{sign} eq $nan) ||
1616 ($y->{sign} eq $nan) ||
1617 ($z->{sign} eq $nan));
1620 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1622 return $x->bnan() if $x->is_zero() || $y->is_zero();
1623 # result will always be +-inf:
1624 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1625 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1626 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1627 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1628 return $x->binf('-');
1631 return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1632 ((!$x->isa($self)) || (!$y->isa($self)));
1634 # aEb * cEd = (a*c)E(b+d)
1635 $MBI->_mul($x->{_m},$y->{_m});
1636 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1638 $r[3] = $y; # no push!
1641 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1643 # z=inf handling (z=NaN handled above)
1644 $x->{sign} = $z->{sign}, return $x if $z->{sign} =~ /^[+-]inf$/;
1646 # take lower of the two e's and adapt m1 to it to match m2
1648 $e = $MBI->_zero() if !defined $e; # if no BFLOAT?
1649 $e = $MBI->_copy($e); # make copy (didn't do it yet)
1653 ($e,$es) = _e_sub($e, $x->{_e}, $z->{_es} || '+', $x->{_es});
1655 my $add = $MBI->_copy($z->{_m});
1657 if ($es eq '-') # < 0
1659 $MBI->_lsft( $x->{_m}, $e, 10);
1660 ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
1662 elsif (!$MBI->_is_zero($e)) # > 0
1664 $MBI->_lsft($add, $e, 10);
1666 # else: both e are the same, so just leave them
1668 if ($x->{sign} eq $z->{sign})
1671 $x->{_m} = $MBI->_add($x->{_m}, $add);
1675 ($x->{_m}, $x->{sign}) =
1676 _e_add($x->{_m}, $add, $x->{sign}, $z->{sign});
1679 # delete trailing zeros, then round
1680 $x->bnorm()->round(@r);
1685 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1686 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1689 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1690 # objectify is costly, so avoid it
1691 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1693 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1696 return $x if $x->modify('bdiv');
1698 return $self->_div_inf($x,$y)
1699 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1701 # x== 0 # also: or y == 1 or y == -1
1702 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1705 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1707 # we need to limit the accuracy to protect against overflow
1709 my (@params,$scale);
1710 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1712 return $x if $x->is_nan(); # error in _find_round_parameters?
1714 # no rounding at all, so must use fallback
1715 if (scalar @params == 0)
1717 # simulate old behaviour
1718 $params[0] = $self->div_scale(); # and round to it as accuracy
1719 $scale = $params[0]+4; # at least four more for proper round
1720 $params[2] = $r; # round mode by caller or undef
1721 $fallback = 1; # to clear a/p afterwards
1725 # the 4 below is empirical, and there might be cases where it is not
1727 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1730 my $rem; $rem = $self->bzero() if wantarray;
1732 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1734 my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
1735 $scale = $lx if $lx > $scale;
1736 $scale = $ly if $ly > $scale;
1737 my $diff = $ly - $lx;
1738 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
1740 # already handled inf/NaN/-inf above:
1742 # check that $y is not 1 nor -1 and cache the result:
1743 my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));
1745 # flipping the sign of $y will also flip the sign of $x for the special
1746 # case of $x->bsub($x); so we can catch it below:
1747 my $xsign = $x->{sign};
1748 $y->{sign} =~ tr/+-/-+/;
1750 if ($xsign ne $x->{sign})
1752 # special case of $x /= $x results in 1
1753 $x->bone(); # "fixes" also sign of $y, since $x is $y
1757 # correct $y's sign again
1758 $y->{sign} =~ tr/+-/-+/;
1759 # continue with normal div code:
1761 # make copy of $x in case of list context for later reminder calculation
1762 if (wantarray && $y_not_one)
1767 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1769 # check for / +-1 ( +/- 1E0)
1772 # promote BigInts and it's subclasses (except when already a BigFloat)
1773 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1775 # calculate the result to $scale digits and then round it
1776 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1777 $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1778 $MBI->_div ($x->{_m},$y->{_m}); # a/c
1780 # correct exponent of $x
1781 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1782 # correct for 10**scale
1783 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1784 $x->bnorm(); # remove trailing 0's
1786 } # ende else $x != $y
1788 # shortcut to not run through _find_round_parameters again
1789 if (defined $params[0])
1791 delete $x->{_a}; # clear before round
1792 $x->bround($params[0],$params[2]); # then round accordingly
1796 delete $x->{_p}; # clear before round
1797 $x->bfround($params[1],$params[2]); # then round accordingly
1801 # clear a/p after round, since user did not request it
1802 delete $x->{_a}; delete $x->{_p};
1809 $rem->bmod($y,@params); # copy already done
1813 # clear a/p after round, since user did not request it
1814 delete $rem->{_a}; delete $rem->{_p};
1823 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1826 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1827 # objectify is costly, so avoid it
1828 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1830 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1833 return $x if $x->modify('bmod');
1835 # handle NaN, inf, -inf
1836 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1838 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1839 $x->{sign} = $re->{sign};
1840 $x->{_e} = $re->{_e};
1841 $x->{_m} = $re->{_m};
1842 return $x->round($a,$p,$r,$y);
1846 return $x->bnan() if $x->is_zero();
1850 return $x->bzero() if $x->is_zero()
1852 # check that $y == +1 or $y == -1:
1853 ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));
1855 my $cmp = $x->bacmp($y); # equal or $x < $y?
1856 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1858 # only $y of the operands negative?
1859 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1861 $x->{sign} = $y->{sign}; # calc sign first
1862 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1864 my $ym = $MBI->_copy($y->{_m});
1867 $MBI->_lsft( $ym, $y->{_e}, 10)
1868 if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1870 # if $y has digits after dot
1871 my $shifty = 0; # correct _e of $x by this
1872 if ($y->{_es} eq '-') # has digits after dot
1874 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1875 $shifty = $MBI->_num($y->{_e}); # no more digits after dot
1876 $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1878 # $ym is now mantissa of $y based on exponent 0
1880 my $shiftx = 0; # correct _e of $x by this
1881 if ($x->{_es} eq '-') # has digits after dot
1883 # 123.4 % 20 => 1234 % 200
1884 $shiftx = $MBI->_num($x->{_e}); # no more digits after dot
1885 $MBI->_lsft($ym, $x->{_e}, 10); # 123 => 1230
1887 # 123e1 % 20 => 1230 % 20
1888 if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1890 $MBI->_lsft( $x->{_m}, $x->{_e},10); # es => '+' here
1893 $x->{_e} = $MBI->_new($shiftx);
1895 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1896 $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1898 # now mantissas are equalized, exponent of $x is adjusted, so calc result
1900 $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1902 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1905 if ($neg != 0) # one of them negative => correct in place
1908 $x->{_m} = $r->{_m};
1909 $x->{_e} = $r->{_e};
1910 $x->{_es} = $r->{_es};
1911 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1915 $x->round($a,$p,$r,$y); # round and return
1920 # calculate $y'th root of $x
1923 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1924 # objectify is costly, so avoid it
1925 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1927 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1930 return $x if $x->modify('broot');
1932 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1933 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1934 $y->{sign} !~ /^\+$/;
1936 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1938 # we need to limit the accuracy to protect against overflow
1940 my (@params,$scale);
1941 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1943 return $x if $x->is_nan(); # error in _find_round_parameters?
1945 # no rounding at all, so must use fallback
1946 if (scalar @params == 0)
1948 # simulate old behaviour
1949 $params[0] = $self->div_scale(); # and round to it as accuracy
1950 $scale = $params[0]+4; # at least four more for proper round
1951 $params[2] = $r; # iound mode by caller or undef
1952 $fallback = 1; # to clear a/p afterwards
1956 # the 4 below is empirical, and there might be cases where it is not
1958 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1961 # when user set globals, they would interfere with our calculation, so
1962 # disable them and later re-enable them
1964 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1965 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1966 # we also need to disable any set A or P on $x (_find_round_parameters took
1967 # them already into account), since these would interfere, too
1968 delete $x->{_a}; delete $x->{_p};
1969 # need to disable $upgrade in BigInt, to avoid deep recursion
1970 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1972 # remember sign and make $x positive, since -4 ** (1/2) => -2
1973 my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
1976 if ($y->isa('Math::BigFloat'))
1978 $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
1982 $is_two = ($y == 2);
1985 # normal square root if $y == 2:
1988 $x->bsqrt($scale+4);
1990 elsif ($y->is_one('-'))
1993 my $u = $self->bone()->bdiv($x,$scale);
1994 # copy private parts over
1995 $x->{_m} = $u->{_m};
1996 $x->{_e} = $u->{_e};
1997 $x->{_es} = $u->{_es};
2001 # calculate the broot() as integer result first, and if it fits, return
2002 # it rightaway (but only if $x and $y are integer):
2004 my $done = 0; # not yet
2005 if ($y->is_int() && $x->is_int())
2007 my $i = $MBI->_copy( $x->{_m} );
2008 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2009 my $int = Math::BigInt->bzero();
2011 $int->broot($y->as_number());
2013 if ($int->copy()->bpow($y) == $x)
2015 # found result, return it
2016 $x->{_m} = $int->{value};
2017 $x->{_e} = $MBI->_zero();
2025 my $u = $self->bone()->bdiv($y,$scale+4);
2026 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
2027 $x->bpow($u,$scale+4); # el cheapo
2030 $x->bneg() if $sign == 1;
2032 # shortcut to not run through _find_round_parameters again
2033 if (defined $params[0])
2035 $x->bround($params[0],$params[2]); # then round accordingly
2039 $x->bfround($params[1],$params[2]); # then round accordingly
2043 # clear a/p after round, since user did not request it
2044 delete $x->{_a}; delete $x->{_p};
2047 $$abr = $ab; $$pbr = $pb;
2053 # calculate square root
2054 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2056 return $x if $x->modify('bsqrt');
2058 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
2059 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
2060 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
2062 # we need to limit the accuracy to protect against overflow
2064 my (@params,$scale);
2065 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2067 return $x if $x->is_nan(); # error in _find_round_parameters?
2069 # no rounding at all, so must use fallback
2070 if (scalar @params == 0)
2072 # simulate old behaviour
2073 $params[0] = $self->div_scale(); # and round to it as accuracy
2074 $scale = $params[0]+4; # at least four more for proper round
2075 $params[2] = $r; # round mode by caller or undef
2076 $fallback = 1; # to clear a/p afterwards
2080 # the 4 below is empirical, and there might be cases where it is not
2082 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2085 # when user set globals, they would interfere with our calculation, so
2086 # disable them and later re-enable them
2088 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2089 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2090 # we also need to disable any set A or P on $x (_find_round_parameters took
2091 # them already into account), since these would interfere, too
2092 delete $x->{_a}; delete $x->{_p};
2093 # need to disable $upgrade in BigInt, to avoid deep recursion
2094 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2096 my $i = $MBI->_copy( $x->{_m} );
2097 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2098 my $xas = Math::BigInt->bzero();
2101 my $gs = $xas->copy()->bsqrt(); # some guess
2103 if (($x->{_es} ne '-') # guess can't be accurate if there are
2104 # digits after the dot
2105 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
2107 # exact result, copy result over to keep $x
2108 $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
2110 # shortcut to not run through _find_round_parameters again
2111 if (defined $params[0])
2113 $x->bround($params[0],$params[2]); # then round accordingly
2117 $x->bfround($params[1],$params[2]); # then round accordingly
2121 # clear a/p after round, since user did not request it
2122 delete $x->{_a}; delete $x->{_p};
2124 # re-enable A and P, upgrade is taken care of by "local"
2125 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
2129 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
2130 # of the result by multipyling the input by 100 and then divide the integer
2131 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
2133 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
2134 my $y1 = $MBI->_copy($x->{_m});
2136 my $length = $MBI->_len($y1);
2138 # Now calculate how many digits the result of sqrt(y1) would have
2139 my $digits = int($length / 2);
2141 # But we need at least $scale digits, so calculate how many are missing
2142 my $shift = $scale - $digits;
2144 # That should never happen (we take care of integer guesses above)
2145 # $shift = 0 if $shift < 0;
2147 # Multiply in steps of 100, by shifting left two times the "missing" digits
2148 my $s2 = $shift * 2;
2150 # We now make sure that $y1 has the same odd or even number of digits than
2151 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
2152 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
2153 # steps of 10. The length of $x does not count, since an even or odd number
2154 # of digits before the dot is not changed by adding an even number of digits
2155 # after the dot (the result is still odd or even digits long).
2156 $s2++ if $MBI->_is_odd($x->{_e});
2158 $MBI->_lsft( $y1, $MBI->_new($s2), 10);
2160 # now take the square root and truncate to integer
2161 $y1 = $MBI->_sqrt($y1);
2163 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
2164 # result, which is than later rounded to the desired scale.
2166 # calculate how many zeros $x had after the '.' (or before it, depending
2167 # on sign of $dat, the result should have half as many:
2168 my $dat = $MBI->_num($x->{_e});
2169 $dat = -$dat if $x->{_es} eq '-';
2174 # no zeros after the dot (e.g. 1.23, 0.49 etc)
2175 # preserve half as many digits before the dot than the input had
2176 # (but round this "up")
2177 $dat = int(($dat+1)/2);
2181 $dat = int(($dat)/2);
2183 $dat -= $MBI->_len($y1);
2187 $x->{_e} = $MBI->_new( $dat );
2192 $x->{_e} = $MBI->_new( $dat );
2198 # shortcut to not run through _find_round_parameters again
2199 if (defined $params[0])
2201 $x->bround($params[0],$params[2]); # then round accordingly
2205 $x->bfround($params[1],$params[2]); # then round accordingly
2209 # clear a/p after round, since user did not request it
2210 delete $x->{_a}; delete $x->{_p};
2213 $$abr = $ab; $$pbr = $pb;
2219 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2220 # compute factorial number, modifies first argument
2223 my ($self,$x,@r) = (ref($_[0]),@_);
2224 # objectify is costly, so avoid it
2225 ($self,$x,@r) = objectify(1,@_) if !ref($x);
2228 return $x if $x->modify('bfac') || $x->{sign} eq '+inf';
2231 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
2232 ($x->{_es} ne '+')); # digits after dot?
2234 # use BigInt's bfac() for faster calc
2235 if (! $MBI->_is_zero($x->{_e}))
2237 $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
2238 $x->{_e} = $MBI->_zero(); # normalize
2241 $MBI->_fac($x->{_m}); # calculate factorial
2242 $x->bnorm()->round(@r); # norm again and round result
2247 # Calculate a power where $y is a non-integer, like 2 ** 0.3
2251 # if $y == 0.5, it is sqrt($x)
2252 $HALF = $self->new($HALF) unless ref($HALF);
2253 return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0;
2256 # a ** x == e ** (x * ln a)
2260 # Taylor: | u u^2 u^3 |
2261 # x ** y = 1 + | --- + --- + ----- + ... |
2264 # we need to limit the accuracy to protect against overflow
2266 my ($scale,@params);
2267 ($x,@params) = $x->_find_round_parameters(@r);
2269 return $x if $x->is_nan(); # error in _find_round_parameters?
2271 # no rounding at all, so must use fallback
2272 if (scalar @params == 0)
2274 # simulate old behaviour
2275 $params[0] = $self->div_scale(); # and round to it as accuracy
2276 $params[1] = undef; # disable P
2277 $scale = $params[0]+4; # at least four more for proper round
2278 $params[2] = $r[2]; # round mode by caller or undef
2279 $fallback = 1; # to clear a/p afterwards
2283 # the 4 below is empirical, and there might be cases where it is not
2285 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2288 # when user set globals, they would interfere with our calculation, so
2289 # disable them and later re-enable them
2291 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2292 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2293 # we also need to disable any set A or P on $x (_find_round_parameters took
2294 # them already into account), since these would interfere, too
2295 delete $x->{_a}; delete $x->{_p};
2296 # need to disable $upgrade in BigInt, to avoid deep recursion
2297 local $Math::BigInt::upgrade = undef;
2299 my ($limit,$v,$u,$below,$factor,$next,$over);
2301 $u = $x->copy()->blog(undef,$scale)->bmul($y);
2302 $v = $self->bone(); # 1
2303 $factor = $self->new(2); # 2
2304 $x->bone(); # first term: 1
2306 $below = $v->copy();
2309 $limit = $self->new("1E-". ($scale-1));
2313 # we calculate the next term, and add it to the last
2314 # when the next term is below our limit, it won't affect the outcome
2315 # anymore, so we stop:
2316 $next = $over->copy()->bdiv($below,$scale);
2317 last if $next->bacmp($limit) <= 0;
2319 # calculate things for the next term
2320 $over *= $u; $below *= $factor; $factor->binc();
2322 last if $x->{sign} !~ /^[-+]$/;
2327 # shortcut to not run through _find_round_parameters again
2328 if (defined $params[0])
2330 $x->bround($params[0],$params[2]); # then round accordingly
2334 $x->bfround($params[1],$params[2]); # then round accordingly
2338 # clear a/p after round, since user did not request it
2339 delete $x->{_a}; delete $x->{_p};
2342 $$abr = $ab; $$pbr = $pb;
2348 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2349 # compute power of two numbers, second arg is used as integer
2350 # modifies first argument
2353 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2354 # objectify is costly, so avoid it
2355 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2357 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2360 return $x if $x->modify('bpow');
2362 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
2363 return $x if $x->{sign} =~ /^[+-]inf$/;
2365 # cache the result of is_zero
2366 my $y_is_zero = $y->is_zero();
2367 return $x->bone() if $y_is_zero;
2368 return $x if $x->is_one() || $y->is_one();
2370 my $x_is_zero = $x->is_zero();
2371 return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int(); # non-integer power
2373 my $y1 = $y->as_number()->{value}; # make MBI part
2376 if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
2378 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
2379 return $MBI->_is_odd($y1) ? $x : $x->babs(1);
2383 #return $x->bone() if $y_is_zero;
2384 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
2385 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
2390 $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
2392 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
2393 $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
2394 $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
2396 $x->{sign} = $new_sign;
2398 if ($y->{sign} eq '-')
2400 # modify $x in place!
2401 my $z = $x->copy(); $x->bone();
2402 return scalar $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
2404 $x->round($a,$p,$r,$y);
2409 # takes a very large number to a very large exponent in a given very
2410 # large modulus, quickly, thanks to binary exponentation. Supports
2411 # negative exponents.
2412 my ($self,$num,$exp,$mod,@r) = objectify(3,@_);
2414 return $num if $num->modify('bmodpow');
2416 # check modulus for valid values
2417 return $num->bnan() if ($mod->{sign} ne '+' # NaN, - , -inf, +inf
2418 || $mod->is_zero());
2420 # check exponent for valid values
2421 if ($exp->{sign} =~ /\w/)
2423 # i.e., if it's NaN, +inf, or -inf...
2424 return $num->bnan();
2427 $num->bmodinv ($mod) if ($exp->{sign} eq '-');
2429 # check num for valid values (also NaN if there was no inverse but $exp < 0)
2430 return $num->bnan() if $num->{sign} !~ /^[+-]$/;
2432 # $mod is positive, sign on $exp is ignored, result also positive
2434 # XXX TODO: speed it up when all three numbers are integers
2435 $num->bpow($exp)->bmod($mod);
2438 ###############################################################################
2439 # trigonometric functions
2441 # helper function for bpi()
2445 my ($a, $s, $b) = @_;
2449 # $a and $b are negativ: -> add
2454 my $c = $MBI->_acmp($a,$b);
2455 # $a positiv, $b negativ
2463 $a = $MBI->_sub( $MBI->_copy($b), $a);
2473 my ($a, $s, $b) = @_;
2477 # $a and $b are positiv: -> add
2482 my $c = $MBI->_acmp($a,$b);
2483 # $a positiv, $b negativ
2491 $a = $MBI->_sub( $MBI->_copy($b), $a);
2501 # Calculate PI to N digits (including the 3 before the dot).
2503 # The basic algorithm is the one implemented in:
2505 # The Computer Language Shootout
2506 # http://shootout.alioth.debian.org/
2508 # contributed by Robert Bradshaw
2509 # modified by Ruud H.G.van Tol
2510 # modified by Emanuele Zeppieri
2512 # We re-implement it here by using the low-level library directly. Also,
2513 # the functions consume() and extract_digit() were inlined and some
2514 # rendundand operations ( like *= 1 ) were removed.
2523 # called like Math::BigFloat::bpi(10);
2524 $n = $self; $self = $class;
2525 # called like Math::BigFloat->bpi();
2526 $n = undef if $n eq 'Math::BigFloat';
2528 $self = ref($self) if ref($self);
2529 $n = 40 if !defined $n || $n < 1;
2531 my $z0 = $MBI->_one();
2532 my $z1 = $MBI->_zero();
2533 my $z2 = $MBI->_one();
2534 my $ten = $MBI->_ten();
2535 my $three = $MBI->_new(3);
2536 my ($s, $d, $e, $r); my $k = 0; my $z1_sign = 0;
2543 if ($MBI->_acmp($z0,$z2) != 1)
2545 # my $o = $z0 * 3 + $z1;
2546 my $o = $MBI->_mul( $MBI->_copy($z0), $three);
2547 $z1_sign == 0 ? $MBI->_sub( $o, $z1) : $MBI->_add( $o, $z1);
2548 ($d,$r) = $MBI->_div( $MBI->_copy($o), $z2 );
2549 $d = $MBI->_num($d);
2550 $e = $MBI->_num( scalar $MBI->_div( $MBI->_add($o, $z0), $z2 ) );
2554 my $k2 = $MBI->_new( 2*$k+1 );
2555 # mul works regardless of the sign of $z1 since $k is always positive
2556 $MBI->_mul( $z1, $k2 );
2557 ($z1, $z1_sign) = _signed_add($z1, $z1_sign, $MBI->_mul( $MBI->_new(4*$k+2), $z0 ) );
2558 $MBI->_mul( $z0, $MBI->_new($k) );
2559 $MBI->_mul( $z2, $k2 );
2561 $MBI->_mul( $z1, $ten );
2562 ($z1, $z1_sign) = _signed_sub($z1, $z1_sign, $MBI->_mul( $MBI->_new(10*$d), $z2 ) );
2563 $MBI->_mul( $z0, $ten );
2567 my $x = $self->new(0);
2569 $x->{_e} = $MBI->_new(length($s)-1);
2570 $x->{_m} = $MBI->_new($s);
2577 # Calculate a cosinus of x.
2578 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2580 # Taylor: x^2 x^4 x^6 x^8
2581 # cos = 1 - --- + --- - --- + --- ...
2584 # we need to limit the accuracy to protect against overflow
2586 my ($scale,@params);
2587 ($x,@params) = $x->_find_round_parameters(@r);
2589 # constant object or error in _find_round_parameters?
2590 return $x if $x->modify('bcos') || $x->is_nan();
2592 return $x->bone(@r) if $x->is_zero();
2594 # no rounding at all, so must use fallback
2595 if (scalar @params == 0)
2597 # simulate old behaviour
2598 $params[0] = $self->div_scale(); # and round to it as accuracy
2599 $params[1] = undef; # disable P
2600 $scale = $params[0]+4; # at least four more for proper round
2601 $params[2] = $r[2]; # round mode by caller or undef
2602 $fallback = 1; # to clear a/p afterwards
2606 # the 4 below is empirical, and there might be cases where it is not
2608 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2611 # when user set globals, they would interfere with our calculation, so
2612 # disable them and later re-enable them
2614 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2615 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2616 # we also need to disable any set A or P on $x (_find_round_parameters took
2617 # them already into account), since these would interfere, too
2618 delete $x->{_a}; delete $x->{_p};
2619 # need to disable $upgrade in BigInt, to avoid deep recursion
2620 local $Math::BigInt::upgrade = undef;
2623 my $over = $x * $x; # X ^ 2
2624 my $x2 = $over->copy(); # X ^ 2; difference between terms
2625 my $sign = 1; # start with -=
2626 my $below = $self->new(2); my $factorial = $self->new(3);
2627 $x->bone(); delete $x->{a}; delete $x->{p};
2629 my $limit = $self->new("1E-". ($scale-1));
2633 # we calculate the next term, and add it to the last
2634 # when the next term is below our limit, it won't affect the outcome
2635 # anymore, so we stop:
2636 my $next = $over->copy()->bdiv($below,$scale);
2637 last if $next->bacmp($limit) <= 0;
2647 $sign = 1-$sign; # alternate
2648 # calculate things for the next term
2649 $over->bmul($x2); # $x*$x
2650 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2651 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2654 # shortcut to not run through _find_round_parameters again
2655 if (defined $params[0])
2657 $x->bround($params[0],$params[2]); # then round accordingly
2661 $x->bfround($params[1],$params[2]); # then round accordingly
2665 # clear a/p after round, since user did not request it
2666 delete $x->{_a}; delete $x->{_p};
2669 $$abr = $ab; $$pbr = $pb;
2675 # Calculate a sinus of x.
2676 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2678 # taylor: x^3 x^5 x^7 x^9
2679 # sin = x - --- + --- - --- + --- ...
2682 # we need to limit the accuracy to protect against overflow
2684 my ($scale,@params);
2685 ($x,@params) = $x->_find_round_parameters(@r);
2687 # constant object or error in _find_round_parameters?
2688 return $x if $x->modify('bsin') || $x->is_nan();
2690 return $x->bzero(@r) if $x->is_zero();
2692 # no rounding at all, so must use fallback
2693 if (scalar @params == 0)
2695 # simulate old behaviour
2696 $params[0] = $self->div_scale(); # and round to it as accuracy
2697 $params[1] = undef; # disable P
2698 $scale = $params[0]+4; # at least four more for proper round
2699 $params[2] = $r[2]; # round mode by caller or undef
2700 $fallback = 1; # to clear a/p afterwards
2704 # the 4 below is empirical, and there might be cases where it is not
2706 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2709 # when user set globals, they would interfere with our calculation, so
2710 # disable them and later re-enable them
2712 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2713 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2714 # we also need to disable any set A or P on $x (_find_round_parameters took
2715 # them already into account), since these would interfere, too
2716 delete $x->{_a}; delete $x->{_p};
2717 # need to disable $upgrade in BigInt, to avoid deep recursion
2718 local $Math::BigInt::upgrade = undef;
2721 my $over = $x * $x; # X ^ 2
2722 my $x2 = $over->copy(); # X ^ 2; difference between terms
2723 $over->bmul($x); # X ^ 3 as starting value
2724 my $sign = 1; # start with -=
2725 my $below = $self->new(6); my $factorial = $self->new(4);
2726 delete $x->{a}; delete $x->{p};
2728 my $limit = $self->new("1E-". ($scale-1));
2732 # we calculate the next term, and add it to the last
2733 # when the next term is below our limit, it won't affect the outcome
2734 # anymore, so we stop:
2735 my $next = $over->copy()->bdiv($below,$scale);
2736 last if $next->bacmp($limit) <= 0;
2746 $sign = 1-$sign; # alternate
2747 # calculate things for the next term
2748 $over->bmul($x2); # $x*$x
2749 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2750 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2753 # shortcut to not run through _find_round_parameters again
2754 if (defined $params[0])
2756 $x->bround($params[0],$params[2]); # then round accordingly
2760 $x->bfround($params[1],$params[2]); # then round accordingly
2764 # clear a/p after round, since user did not request it
2765 delete $x->{_a}; delete $x->{_p};
2768 $$abr = $ab; $$pbr = $pb;
2774 # Calculate a arcus tangens of x.
2778 # taylor: x^3 x^5 x^7 x^9
2779 # atan = x - --- + --- - --- + --- ...
2783 # This series is only valid if -1 < x < 1, so for other x we need to
2784 # find a different way.
2786 if ($x < -1 || $x > 1)
2788 die("$x is out of range for batan()!");
2791 # we need to limit the accuracy to protect against overflow
2793 my ($scale,@params);
2794 ($x,@params) = $x->_find_round_parameters(@r);
2796 # constant object or error in _find_round_parameters?
2797 return $x if $x->modify('batan') || $x->is_nan();
2799 # no rounding at all, so must use fallback
2800 if (scalar @params == 0)
2802 # simulate old behaviour
2803 $params[0] = $self->div_scale(); # and round to it as accuracy
2804 $params[1] = undef; # disable P
2805 $scale = $params[0]+4; # at least four more for proper round
2806 $params[2] = $r[2]; # round mode by caller or undef
2807 $fallback = 1; # to clear a/p afterwards
2811 # the 4 below is empirical, and there might be cases where it is not
2813 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2816 # when user set globals, they would interfere with our calculation, so
2817 # disable them and later re-enable them
2819 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2820 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2821 # we also need to disable any set A or P on $x (_find_round_parameters took
2822 # them already into account), since these would interfere, too
2823 delete $x->{_a}; delete $x->{_p};
2824 # need to disable $upgrade in BigInt, to avoid deep recursion
2825 local $Math::BigInt::upgrade = undef;
2828 my $over = $x * $x; # X ^ 2
2829 my $x2 = $over->copy(); # X ^ 2; difference between terms
2830 $over->bmul($x); # X ^ 3 as starting value
2831 my $sign = 1; # start with -=
2832 my $below = $self->new(3);
2833 my $two = $self->new(2);
2834 $x->bone(); delete $x->{a}; delete $x->{p};
2836 my $limit = $self->new("1E-". ($scale-1));
2840 # we calculate the next term, and add it to the last
2841 # when the next term is below our limit, it won't affect the outcome
2842 # anymore, so we stop:
2843 my $next = $over->copy()->bdiv($below,$scale);
2844 last if $next->bacmp($limit) <= 0;
2854 $sign = 1-$sign; # alternate
2855 # calculate things for the next term
2856 $over->bmul($x2); # $x*$x
2857 $below->badd($two); # n += 2
2860 # shortcut to not run through _find_round_parameters again
2861 if (defined $params[0])
2863 $x->bround($params[0],$params[2]); # then round accordingly
2867 $x->bfround($params[1],$params[2]); # then round accordingly
2871 # clear a/p after round, since user did not request it
2872 delete $x->{_a}; delete $x->{_p};
2875 $$abr = $ab; $$pbr = $pb;
2879 ###############################################################################
2880 # rounding functions
2884 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
2885 # $n == 0 means round to integer
2886 # expects and returns normalized numbers!
2887 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2889 my ($scale,$mode) = $x->_scale_p(@_);
2890 return $x if !defined $scale || $x->modify('bfround'); # no-op
2892 # never round a 0, +-inf, NaN
2895 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
2898 return $x if $x->{sign} !~ /^[+-]$/;
2900 # don't round if x already has lower precision
2901 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
2903 $x->{_p} = $scale; # remember round in any case
2904 delete $x->{_a}; # and clear A
2907 # round right from the '.'
2909 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round
2911 $scale = -$scale; # positive for simplicity
2912 my $len = $MBI->_len($x->{_m}); # length of mantissa
2914 # the following poses a restriction on _e, but if _e is bigger than a
2915 # scalar, you got other problems (memory etc) anyway
2916 my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e}))); # digits after dot
2917 my $zad = 0; # zeros after dot
2918 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
2920 # p rint "scale $scale dad $dad zad $zad len $len\n";
2921 # number bsstr len zad dad
2922 # 0.123 123e-3 3 0 3
2923 # 0.0123 123e-4 3 1 4
2926 # 1.2345 12345e-4 5 0 4
2928 # do not round after/right of the $dad
2929 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
2931 # round to zero if rounding inside the $zad, but not for last zero like:
2932 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
2933 return $x->bzero() if $scale < $zad;
2934 if ($scale == $zad) # for 0.006, scale -3 and trunc
2940 # adjust round-point to be inside mantissa
2943 $scale = $scale-$zad;
2947 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
2948 $scale = $dbd+$scale;
2954 # round left from the '.'
2956 # 123 => 100 means length(123) = 3 - $scale (2) => 1
2958 my $dbt = $MBI->_len($x->{_m});
2960 my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
2961 # should be the same, so treat it as this
2962 $scale = 1 if $scale == 0;
2963 # shortcut if already integer
2964 return $x if $scale == 1 && $dbt <= $dbd;
2965 # maximum digits before dot
2970 # not enough digits before dot, so round to zero
2973 elsif ( $scale == $dbd )
2980 $scale = $dbd - $scale;
2983 # pass sign to bround for rounding modes '+inf' and '-inf'
2984 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2985 $m->bround($scale,$mode);
2986 $x->{_m} = $m->{value}; # get our mantissa back
2992 # accuracy: preserve $N digits, and overwrite the rest with 0's
2993 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2995 if (($_[0] || 0) < 0)
2997 require Carp; Carp::croak ('bround() needs positive accuracy');
3000 my ($scale,$mode) = $x->_scale_a(@_);
3001 return $x if !defined $scale || $x->modify('bround'); # no-op
3003 # scale is now either $x->{_a}, $accuracy, or the user parameter
3004 # test whether $x already has lower accuracy, do nothing in this case
3005 # but do round if the accuracy is the same, since a math operation might
3006 # want to round a number with A=5 to 5 digits afterwards again
3007 return $x if defined $x->{_a} && $x->{_a} < $scale;
3009 # scale < 0 makes no sense
3010 # scale == 0 => keep all digits
3011 # never round a +-inf, NaN
3012 return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
3014 # 1: never round a 0
3015 # 2: if we should keep more digits than the mantissa has, do nothing
3016 if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
3018 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
3022 # pass sign to bround for '+inf' and '-inf' rounding modes
3023 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
3025 $m->bround($scale,$mode); # round mantissa
3026 $x->{_m} = $m->{value}; # get our mantissa back
3027 $x->{_a} = $scale; # remember rounding
3028 delete $x->{_p}; # and clear P
3029 $x->bnorm(); # del trailing zeros gen. by bround()
3034 # return integer less or equal then $x
3035 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3037 return $x if $x->modify('bfloor');
3039 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3041 # if $x has digits after dot
3042 if ($x->{_es} eq '-')
3044 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3045 $x->{_e} = $MBI->_zero(); # trunc/norm
3046 $x->{_es} = '+'; # abs e
3047 $MBI->_inc($x->{_m}) if $x->{sign} eq '-'; # increment if negative
3049 $x->round($a,$p,$r);
3054 # return integer greater or equal then $x
3055 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3057 return $x if $x->modify('bceil');
3058 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3060 # if $x has digits after dot
3061 if ($x->{_es} eq '-')
3063 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3064 $x->{_e} = $MBI->_zero(); # trunc/norm
3065 $x->{_es} = '+'; # abs e
3066 $MBI->_inc($x->{_m}) if $x->{sign} eq '+'; # increment if positive
3068 $x->round($a,$p,$r);
3073 # shift right by $y (divide by power of $n)
3076 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3077 # objectify is costly, so avoid it
3078 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3080 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3083 return $x if $x->modify('brsft');
3084 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3086 $n = 2 if !defined $n; $n = $self->new($n);
3089 return $x->blsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3091 # the following call to bdiv() will return either quo or (quo,reminder):
3092 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
3097 # shift left by $y (multiply by power of $n)
3100 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3101 # objectify is costly, so avoid it
3102 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3104 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3107 return $x if $x->modify('blsft');
3108 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3110 $n = 2 if !defined $n; $n = $self->new($n);
3113 return $x->brsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3115 $x->bmul($n->bpow($y),$a,$p,$r,$y);
3118 ###############################################################################
3122 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
3127 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
3128 # or falling back to MBI::bxxx()
3129 my $name = $AUTOLOAD;
3131 $name =~ s/(.*):://; # split package
3132 my $c = $1 || $class;
3134 $c->import() if $IMPORT == 0;
3135 if (!_method_alias($name))
3139 # delayed load of Carp and avoid recursion
3141 Carp::croak ("$c: Can't call a method without name");
3143 if (!_method_hand_up($name))
3145 # delayed load of Carp and avoid recursion
3147 Carp::croak ("Can't call $c\-\>$name, not a valid method");
3149 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
3151 return &{"Math::BigInt"."::$name"}(@_);
3153 my $bname = $name; $bname =~ s/^f/b/;
3161 # return a copy of the exponent
3162 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3164 if ($x->{sign} !~ /^[+-]$/)
3166 my $s = $x->{sign}; $s =~ s/^[+-]//;
3167 return Math::BigInt->new($s); # -inf, +inf => +inf
3169 Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
3174 # return a copy of the mantissa
3175 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3177 if ($x->{sign} !~ /^[+-]$/)
3179 my $s = $x->{sign}; $s =~ s/^[+]//;
3180 return Math::BigInt->new($s); # -inf, +inf => +inf
3182 my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
3183 $m->bneg() if $x->{sign} eq '-';
3190 # return a copy of both the exponent and the mantissa
3191 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3193 if ($x->{sign} !~ /^[+-]$/)
3195 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
3196 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
3198 my $m = Math::BigInt->bzero();
3199 $m->{value} = $MBI->_copy($x->{_m});
3200 $m->bneg() if $x->{sign} eq '-';
3201 ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
3204 ##############################################################################
3205 # private stuff (internal use only)
3211 my $lib = ''; my @a;
3212 my $lib_kind = 'try';
3214 for ( my $i = 0; $i < $l ; $i++)
3216 if ( $_[$i] eq ':constant' )
3218 # This causes overlord er load to step in. 'binary' and 'integer'
3219 # are handled by BigInt.
3220 overload::constant float => sub { $self->new(shift); };
3222 elsif ($_[$i] eq 'upgrade')
3224 # this causes upgrading
3225 $upgrade = $_[$i+1]; # or undef to disable
3228 elsif ($_[$i] eq 'downgrade')
3230 # this causes downgrading
3231 $downgrade = $_[$i+1]; # or undef to disable
3234 elsif ($_[$i] =~ /^(lib|try|only)\z/)
3236 # alternative library
3237 $lib = $_[$i+1] || ''; # default Calc
3238 $lib_kind = $1; # lib, try or only
3241 elsif ($_[$i] eq 'with')
3243 # alternative class for our private parts()
3244 # XXX: no longer supported
3245 # $MBI = $_[$i+1] || 'Math::BigInt';
3254 $lib =~ tr/a-zA-Z0-9,://cd; # restrict to sane characters
3255 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
3256 my $mbilib = eval { Math::BigInt->config()->{lib} };
3257 if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
3259 # MBI already loaded
3260 Math::BigInt->import( $lib_kind, "$lib,$mbilib", 'objectify');
3264 # MBI not loaded, or with ne "Math::BigInt::Calc"
3265 $lib .= ",$mbilib" if defined $mbilib;
3266 $lib =~ s/^,//; # don't leave empty
3268 # replacement library can handle lib statement, but also could ignore it
3270 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
3271 # used in the same script, or eval inside import(). So we require MBI:
3272 require Math::BigInt;
3273 Math::BigInt->import( $lib_kind => $lib, 'objectify' );
3277 require Carp; Carp::croak ("Couldn't load $lib: $! $@");
3279 # find out which one was actually loaded
3280 $MBI = Math::BigInt->config()->{lib};
3282 # register us with MBI to get notified of future lib changes
3283 Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
3285 $self->export_to_level(1,$self,@a); # export wanted functions
3290 # adjust m and e so that m is smallest possible
3291 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
3293 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
3295 my $zeros = $MBI->_zeros($x->{_m}); # correct for trailing zeros
3298 my $z = $MBI->_new($zeros);
3299 $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
3300 if ($x->{_es} eq '-')
3302 if ($MBI->_acmp($x->{_e},$z) >= 0)
3304 $x->{_e} = $MBI->_sub ($x->{_e}, $z);
3305 $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
3309 $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
3315 $x->{_e} = $MBI->_add ($x->{_e}, $z);
3320 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
3321 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
3322 $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
3323 if $MBI->_is_zero($x->{_m});
3326 $x; # MBI bnorm is no-op, so dont call it
3329 ##############################################################################
3333 # return number as hexadecimal string (only for integers defined)
3334 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3336 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
3337 return '0x0' if $x->is_zero();
3339 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
3341 my $z = $MBI->_copy($x->{_m});
3342 if (! $MBI->_is_zero($x->{_e})) # > 0
3344 $MBI->_lsft( $z, $x->{_e},10);
3346 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3352 # return number as binary digit string (only for integers defined)
3353 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3355 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
3356 return '0b0' if $x->is_zero();
3358 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
3360 my $z = $MBI->_copy($x->{_m});
3361 if (! $MBI->_is_zero($x->{_e})) # > 0
3363 $MBI->_lsft( $z, $x->{_e},10);
3365 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3371 # return number as octal digit string (only for integers defined)
3372 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3374 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
3375 return '0' if $x->is_zero();
3377 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
3379 my $z = $MBI->_copy($x->{_m});
3380 if (! $MBI->_is_zero($x->{_e})) # > 0
3382 $MBI->_lsft( $z, $x->{_e},10);
3384 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3390 # return copy as a bigint representation of this BigFloat number
3391 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3393 return $x if $x->modify('as_number');
3395 my $z = $MBI->_copy($x->{_m});
3396 if ($x->{_es} eq '-') # < 0
3398 $MBI->_rsft( $z, $x->{_e},10);
3400 elsif (! $MBI->_is_zero($x->{_e})) # > 0
3402 $MBI->_lsft( $z, $x->{_e},10);
3404 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3411 my $class = ref($x) || $x;
3412 $x = $class->new(shift) unless ref($x);
3414 return 1 if $MBI->_is_zero($x->{_m});
3416 my $len = $MBI->_len($x->{_m});
3417 $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
3421 $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
3432 Math::BigFloat - Arbitrary size floating point math package
3439 my $x = Math::BigFloat->new($str); # defaults to 0
3440 my $y = $x->copy(); # make a true copy
3441 my $nan = Math::BigFloat->bnan(); # create a NotANumber
3442 my $zero = Math::BigFloat->bzero(); # create a +0
3443 my $inf = Math::BigFloat->binf(); # create a +inf
3444 my $inf = Math::BigFloat->binf('-'); # create a -inf
3445 my $one = Math::BigFloat->bone(); # create a +1
3446 my $mone = Math::BigFloat->bone('-'); # create a -1
3448 my $pi = Math::BigFloat->bpi(100); # PI to 100 digits
3450 # the following examples compute their result to 100 digits accuracy:
3451 my $cos = Math::BigFloat->new(1)->bcos(100); # cosinus(1)
3452 my $sin = Math::BigFloat->new(1)->bsin(100); # sinus(1)
3453 my $atan = Math::BigFloat->new(1)->batan(100); # arcus tangens(1)
3456 $x->is_zero(); # true if arg is +0
3457 $x->is_nan(); # true if arg is NaN
3458 $x->is_one(); # true if arg is +1
3459 $x->is_one('-'); # true if arg is -1
3460 $x->is_odd(); # true if odd, false for even
3461 $x->is_even(); # true if even, false for odd
3462 $x->is_pos(); # true if >= 0
3463 $x->is_neg(); # true if < 0
3464 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
3466 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
3467 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
3468 $x->sign(); # return the sign, either +,- or NaN
3469 $x->digit($n); # return the nth digit, counting from right
3470 $x->digit(-$n); # return the nth digit, counting from left
3472 # The following all modify their first argument. If you want to preserve
3473 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
3474 # necessary when mixing $a = $b assignments with non-overloaded math.
3477 $x->bzero(); # set $i to 0
3478 $x->bnan(); # set $i to NaN
3479 $x->bone(); # set $x to +1
3480 $x->bone('-'); # set $x to -1
3481 $x->binf(); # set $x to inf
3482 $x->binf('-'); # set $x to -inf
3484 $x->bneg(); # negation
3485 $x->babs(); # absolute value
3486 $x->bnorm(); # normalize (no-op)
3487 $x->bnot(); # two's complement (bit wise not)
3488 $x->binc(); # increment x by 1
3489 $x->bdec(); # decrement x by 1
3491 $x->badd($y); # addition (add $y to $x)
3492 $x->bsub($y); # subtraction (subtract $y from $x)
3493 $x->bmul($y); # multiplication (multiply $x by $y)
3494 $x->bdiv($y); # divide, set $x to quotient
3495 # return (quo,rem) or quo if scalar
3497 $x->bmod($y); # modulus ($x % $y)
3498 $x->bpow($y); # power of arguments ($x ** $y)
3499 $x->bmodpow($exp,$mod); # modular exponentation (($num**$exp) % $mod))
3500 $x->blsft($y, $n); # left shift by $y places in base $n
3501 $x->brsft($y, $n); # right shift by $y places in base $n
3502 # returns (quo,rem) or quo if in scalar context
3504 $x->blog(); # logarithm of $x to base e (Euler's number)
3505 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
3506 $x->bexp(); # calculate e ** $x where e is Euler's number
3508 $x->band($y); # bit-wise and
3509 $x->bior($y); # bit-wise inclusive or
3510 $x->bxor($y); # bit-wise exclusive or
3511 $x->bnot(); # bit-wise not (two's complement)
3513 $x->bsqrt(); # calculate square-root
3514 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
3515 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
3517 $x->bround($N); # accuracy: preserve $N digits
3518 $x->bfround($N); # precision: round to the $Nth digit
3520 $x->bfloor(); # return integer less or equal than $x
3521 $x->bceil(); # return integer greater or equal than $x
3523 # The following do not modify their arguments:
3525 bgcd(@values); # greatest common divisor
3526 blcm(@values); # lowest common multiplicator
3528 $x->bstr(); # return string
3529 $x->bsstr(); # return string in scientific notation
3531 $x->as_int(); # return $x as BigInt
3532 $x->exponent(); # return exponent as BigInt
3533 $x->mantissa(); # return mantissa as BigInt
3534 $x->parts(); # return (mantissa,exponent) as BigInt
3536 $x->length(); # number of digits (w/o sign and '.')
3537 ($l,$f) = $x->length(); # number of digits, and length of fraction
3539 $x->precision(); # return P of $x (or global, if P of $x undef)
3540 $x->precision($n); # set P of $x to $n
3541 $x->accuracy(); # return A of $x (or global, if A of $x undef)
3542 $x->accuracy($n); # set A $x to $n
3544 # these get/set the appropriate global value for all BigFloat objects
3545 Math::BigFloat->precision(); # Precision
3546 Math::BigFloat->accuracy(); # Accuracy
3547 Math::BigFloat->round_mode(); # rounding mode
3551 All operators (including basic math operations) are overloaded if you
3552 declare your big floating point numbers as
3554 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
3556 Operations with overloaded operators preserve the arguments, which is
3557 exactly what you expect.
3559 =head2 Canonical notation
3561 Input to these routines are either BigFloat objects, or strings of the
3562 following four forms:
3576 C</^[+-]\d+E[+-]?\d+$/>
3580 C</^[+-]\d*\.\d+E[+-]?\d+$/>
3584 all with optional leading and trailing zeros and/or spaces. Additionally,
3585 numbers are allowed to have an underscore between any two digits.
3587 Empty strings as well as other illegal numbers results in 'NaN'.
3589 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
3590 are always stored in normalized form. On a string, it creates a BigFloat
3595 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
3597 The string output will always have leading and trailing zeros stripped and drop
3598 a plus sign. C<bstr()> will give you always the form with a decimal point,
3599 while C<bsstr()> (s for scientific) gives you the scientific notation.
3601 Input bstr() bsstr()
3603 ' -123 123 123' '-123123123' '-123123123E0'
3604 '00.0123' '0.0123' '123E-4'
3605 '123.45E-2' '1.2345' '12345E-4'
3606 '10E+3' '10000' '1E4'
3608 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
3609 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
3610 return either undef, <0, 0 or >0 and are suited for sort.
3612 Actual math is done by using the class defined with C<with => Class;> (which
3613 defaults to BigInts) to represent the mantissa and exponent.
3615 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
3616 represent the result when input arguments are not numbers, as well as
3617 the result of dividing by zero.
3619 =head2 C<mantissa()>, C<exponent()> and C<parts()>
3621 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
3622 as BigInts such that:
3624 $m = $x->mantissa();
3625 $e = $x->exponent();
3626 $y = $m * ( 10 ** $e );
3627 print "ok\n" if $x == $y;
3629 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
3631 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
3633 Currently the mantissa is reduced as much as possible, favouring higher
3634 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
3635 This might change in the future, so do not depend on it.
3637 =head2 Accuracy vs. Precision
3639 See also: L<Rounding|Rounding>.
3641 Math::BigFloat supports both precision (rounding to a certain place before or
3642 after the dot) and accuracy (rounding to a certain number of digits). For a
3643 full documentation, examples and tips on these topics please see the large
3644 section about rounding in L<Math::BigInt>.
3646 Since things like C<sqrt(2)> or C<1 / 3> must presented with a limited
3647 accuracy lest a operation consumes all resources, each operation produces
3648 no more than the requested number of digits.
3650 If there is no gloabl precision or accuracy set, B<and> the operation in
3651 question was not called with a requested precision or accuracy, B<and> the
3652 input $x has no accuracy or precision set, then a fallback parameter will
3653 be used. For historical reasons, it is called C<div_scale> and can be accessed
3656 $d = Math::BigFloat->div_scale(); # query
3657 Math::BigFloat->div_scale($n); # set to $n digits
3659 The default value for C<div_scale> is 40.
3661 In case the result of one operation has more digits than specified,
3662 it is rounded. The rounding mode taken is either the default mode, or the one
3663 supplied to the operation after the I<scale>:
3665 $x = Math::BigFloat->new(2);
3666 Math::BigFloat->accuracy(5); # 5 digits max
3667 $y = $x->copy()->bdiv(3); # will give 0.66667
3668 $y = $x->copy()->bdiv(3,6); # will give 0.666667
3669 $y = $x->copy()->bdiv(3,6,undef,'odd'); # will give 0.666667
3670 Math::BigFloat->round_mode('zero');
3671 $y = $x->copy()->bdiv(3,6); # will also give 0.666667
3673 Note that C<< Math::BigFloat->accuracy() >> and C<< Math::BigFloat->precision() >>
3674 set the global variables, and thus B<any> newly created number will be subject
3675 to the global rounding B<immediately>. This means that in the examples above, the
3676 C<3> as argument to C<bdiv()> will also get an accuracy of B<5>.
3678 It is less confusing to either calculate the result fully, and afterwards
3679 round it explicitly, or use the additional parameters to the math
3683 $x = Math::BigFloat->new(2);
3684 $y = $x->copy()->bdiv(3);
3685 print $y->bround(5),"\n"; # will give 0.66667
3690 $x = Math::BigFloat->new(2);
3691 $y = $x->copy()->bdiv(3,5); # will give 0.66667
3698 =item ffround ( +$scale )
3700 Rounds to the $scale'th place left from the '.', counting from the dot.
3701 The first digit is numbered 1.
3703 =item ffround ( -$scale )
3705 Rounds to the $scale'th place right from the '.', counting from the dot.
3709 Rounds to an integer.
3711 =item fround ( +$scale )
3713 Preserves accuracy to $scale digits from the left (aka significant digits)
3714 and pads the rest with zeros. If the number is between 1 and -1, the
3715 significant digits count from the first non-zero after the '.'
3717 =item fround ( -$scale ) and fround ( 0 )
3719 These are effectively no-ops.
3723 All rounding functions take as a second parameter a rounding mode from one of
3724 the following: 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'.
3726 The default rounding mode is 'even'. By using
3727 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
3728 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
3729 no longer supported.
3730 The second parameter to the round functions then overrides the default
3733 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
3734 'trunc' as rounding mode to make it equivalent to:
3739 You can override this by passing the desired rounding mode as parameter to
3742 $x = Math::BigFloat->new(2.5);
3743 $y = $x->as_number('odd'); # $y = 3
3747 Math::BigFloat supports all methods that Math::BigInt supports, except it
3748 calculates non-integer results when possible. Please see L<Math::BigInt>
3749 for a full description of each method. Below are just the most important
3754 $x->accuracy(5); # local for $x
3755 CLASS->accuracy(5); # global for all members of CLASS
3756 # Note: This also applies to new()!
3758 $A = $x->accuracy(); # read out accuracy that affects $x
3759 $A = CLASS->accuracy(); # read out global accuracy
3761 Set or get the global or local accuracy, aka how many significant digits the
3762 results have. If you set a global accuracy, then this also applies to new()!
3764 Warning! The accuracy I<sticks>, e.g. once you created a number under the
3765 influence of C<< CLASS->accuracy($A) >>, all results from math operations with
3766 that number will also be rounded.
3768 In most cases, you should probably round the results explicitly using one of
3769 L<round()>, L<bround()> or L<bfround()> or by passing the desired accuracy
3770 to the math operation as additional parameter:
3772 my $x = Math::BigInt->new(30000);
3773 my $y = Math::BigInt->new(7);
3774 print scalar $x->copy()->bdiv($y, 2); # print 4300
3775 print scalar $x->copy()->bdiv($y)->bround(2); # print 4300
3779 $x->precision(-2); # local for $x, round at the second digit right of the dot
3780 $x->precision(2); # ditto, round at the second digit left of the dot
3782 CLASS->precision(5); # Global for all members of CLASS
3783 # This also applies to new()!
3784 CLASS->precision(-5); # ditto
3786 $P = CLASS->precision(); # read out global precision
3787 $P = $x->precision(); # read out precision that affects $x
3789 Note: You probably want to use L<accuracy()> instead. With L<accuracy> you
3790 set the number of digits each result should have, with L<precision> you
3791 set the place where to round!
3795 $x->bexp($accuracy); # calculate e ** X
3797 Calculates the expression C<e ** $x> where C<e> is Euler's number.
3799 This method was added in v1.82 of Math::BigInt (April 2007).
3803 $x->bnok($y); # x over y (binomial coefficient n over k)
3805 Calculates the binomial coefficient n over k, also called the "choose"
3806 function. The result is equivalent to:
3812 This method was added in v1.84 of Math::BigInt (April 2007).
3816 print Math::BigFloat->bpi(100), "\n";
3818 Calculate PI to N digits (including the 3 before the dot).
3820 This method was added in v1.87 of Math::BigInt (June 2007).
3824 my $x = Math::BigFloat->new(1);
3825 print $x->bcos(100), "\n";
3827 Calculate the cosinus of $x, modifying $x in place.
3829 This method was added in v1.87 of Math::BigInt (June 2007).
3833 my $x = Math::BigFloat->new(1);
3834 print $x->bsin(100), "\n";
3836 Calculate the sinus of $x, modifying $x in place.
3838 This method was added in v1.87 of Math::BigInt (June 2007).
3842 my $x = Math::BigFloat->new(1);
3843 print $x->batan(100), "\n";
3845 Calculate the arcus tanges of $x, modifying $x in place.
3847 This method was added in v1.87 of Math::BigInt (June 2007).
3853 Multiply $x by $y, and then add $z to the result.
3855 This method was added in v1.87 of Math::BigInt (June 2007).
3857 =head1 Autocreating constants
3859 After C<use Math::BigFloat ':constant'> all the floating point constants
3860 in the given scope are converted to C<Math::BigFloat>. This conversion
3861 happens at compile time.
3865 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
3867 prints the value of C<2E-100>. Note that without conversion of
3868 constants the expression 2E-100 will be calculated as normal floating point
3871 Please note that ':constant' does not affect integer constants, nor binary
3872 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
3877 Math with the numbers is done (by default) by a module called
3878 Math::BigInt::Calc. This is equivalent to saying:
3880 use Math::BigFloat lib => 'Calc';
3882 You can change this by using:
3884 use Math::BigFloat lib => 'GMP';
3886 Note: The keyword 'lib' will warn when the requested library could not be
3887 loaded. To suppress the warning use 'try' instead:
3889 use Math::BigFloat try => 'GMP';
3891 To turn the warning into a die(), use 'only' instead:
3893 use Math::BigFloat only => 'GMP';
3895 The following would first try to find Math::BigInt::Foo, then
3896 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
3898 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
3900 See the respective low-level library documentation for further details.
3902 Please note that Math::BigFloat does B<not> use the denoted library itself,
3903 but it merely passes the lib argument to Math::BigInt. So, instead of the need
3906 use Math::BigInt lib => 'GMP';
3909 you can roll it all into one line:
3911 use Math::BigFloat lib => 'GMP';
3913 It is also possible to just require Math::BigFloat:
3915 require Math::BigFloat;
3917 This will load the necessary things (like BigInt) when they are needed, and
3920 See L<Math::BigInt> for more details than you ever wanted to know about using
3921 a different low-level library.
3923 =head2 Using Math::BigInt::Lite
3925 For backwards compatibility reasons it is still possible to
3926 request a different storage class for use with Math::BigFloat:
3928 use Math::BigFloat with => 'Math::BigInt::Lite';
3930 However, this request is ignored, as the current code now uses the low-level
3931 math libary for directly storing the number parts.
3935 C<Math::BigFloat> exports nothing by default, but can export the C<bpi()> method:
3937 use Math::BigFloat qw/bpi/;
3939 print bpi(10), "\n";
3943 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
3947 Do not try to be clever to insert some operations in between switching
3950 require Math::BigFloat;
3951 my $matter = Math::BigFloat->bone() + 4; # load BigInt and Calc
3952 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
3953 my $anti_matter = Math::BigFloat->bone()+4; # now use Pari
3955 This will create objects with numbers stored in two different backend libraries,
3956 and B<VERY BAD THINGS> will happen when you use these together:
3958 my $flash_and_bang = $matter + $anti_matter; # Don't do this!
3962 =item stringify, bstr()
3964 Both stringify and bstr() now drop the leading '+'. The old code would return
3965 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
3966 reasoning and details.
3970 The following will probably not print what you expect:
3972 print $c->bdiv(123.456),"\n";
3974 It prints both quotient and reminder since print works in list context. Also,
3975 bdiv() will modify $c, so be careful. You probably want to use
3977 print $c / 123.456,"\n";
3978 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
3984 The following will probably not print what you expect:
3986 my $c = Math::BigFloat->new('3.14159');
3987 print $c->brsft(3,10),"\n"; # prints 0.00314153.1415
3989 It prints both quotient and remainder, since print calls C<brsft()> in list
3990 context. Also, C<< $c->brsft() >> will modify $c, so be careful.
3991 You probably want to use
3993 print scalar $c->copy()->brsft(3,10),"\n";
3994 # or if you really want to modify $c
3995 print scalar $c->brsft(3,10),"\n";
3999 =item Modifying and =
4003 $x = Math::BigFloat->new(5);
4006 It will not do what you think, e.g. making a copy of $x. Instead it just makes
4007 a second reference to the B<same> object and stores it in $y. Thus anything
4008 that modifies $x will modify $y (except overloaded math operators), and vice
4009 versa. See L<Math::BigInt> for details and how to avoid that.
4013 C<bpow()> now modifies the first argument, unlike the old code which left
4014 it alone and only returned the result. This is to be consistent with
4015 C<badd()> etc. The first will modify $x, the second one won't:
4017 print bpow($x,$i),"\n"; # modify $x
4018 print $x->bpow($i),"\n"; # ditto
4019 print $x ** $i,"\n"; # leave $x alone
4021 =item precision() vs. accuracy()
4023 A common pitfall is to use L<precision()> when you want to round a result to
4024 a certain number of digits:
4028 Math::BigFloat->precision(4); # does not do what you think it does
4029 my $x = Math::BigFloat->new(12345); # rounds $x to "12000"!
4030 print "$x\n"; # print "12000"
4031 my $y = Math::BigFloat->new(3); # rounds $y to "0"!
4032 print "$y\n"; # print "0"
4033 $z = $x / $y; # 12000 / 0 => NaN!
4035 print $z->precision(),"\n"; # 4
4037 Replacing L<precision> with L<accuracy> is probably not what you want, either:
4041 Math::BigFloat->accuracy(4); # enables global rounding:
4042 my $x = Math::BigFloat->new(123456); # rounded immediately to "12350"
4043 print "$x\n"; # print "123500"
4044 my $y = Math::BigFloat->new(3); # rounded to "3
4045 print "$y\n"; # print "3"
4046 print $z = $x->copy()->bdiv($y),"\n"; # 41170
4047 print $z->accuracy(),"\n"; # 4
4049 What you want to use instead is:
4053 my $x = Math::BigFloat->new(123456); # no rounding
4054 print "$x\n"; # print "123456"
4055 my $y = Math::BigFloat->new(3); # no rounding
4056 print "$y\n"; # print "3"
4057 print $z = $x->copy()->bdiv($y,4),"\n"; # 41150
4058 print $z->accuracy(),"\n"; # undef
4060 In addition to computing what you expected, the last example also does B<not>
4061 "taint" the result with an accuracy or precision setting, which would
4062 influence any further operation.
4068 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
4069 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
4071 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
4072 because they solve the autoupgrading/downgrading issue, at least partly.
4074 The package at L<http://search.cpan.org/~tels/Math-BigInt> contains
4075 more documentation including a full version history, testcases, empty
4076 subclass files and benchmarks.
4080 This program is free software; you may redistribute it and/or modify it under
4081 the same terms as Perl itself.
4085 Mark Biggar, overloaded interface by Ilya Zakharevich.
4086 Completely rewritten by Tels L<http://bloodgate.com> in 2001 - 2006, and still