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 # one of 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'
61 # the package we are using for our private parts, defaults to:
62 # Math::BigInt->config()->{lib}
63 my $MBI = 'Math::BigInt::FastCalc';
65 # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
67 # the same for infinity
70 # constant for easier life
73 my $IMPORT = 0; # was import() called yet? used to make require work
75 # some digits of accuracy for blog(undef,10); which we use in blog() for speed
77 '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
78 my $LOG_10_A = length($LOG_10)-1;
81 '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
82 my $LOG_2_A = length($LOG_2)-1;
83 my $HALF = '0.5'; # made into an object if nec.
85 ##############################################################################
86 # the old code had $rnd_mode, so we need to support it, too
88 sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
89 sub FETCH { return $round_mode; }
90 sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
94 # when someone sets $rnd_mode, we catch this and check the value to see
95 # whether it is valid or not.
96 $rnd_mode = 'even'; tie $rnd_mode, 'Math::BigFloat';
98 # we need both of them in this package:
99 *as_int = \&as_number;
102 ##############################################################################
105 # valid method aliases for AUTOLOAD
106 my %methods = map { $_ => 1 }
107 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
108 fint facmp fcmp fzero fnan finf finc fdec ffac fneg
109 fceil ffloor frsft flsft fone flog froot fexp
111 # valid methods that can be handed up (for AUTOLOAD)
112 my %hand_ups = map { $_ => 1 }
113 qw / is_nan is_inf is_negative is_positive is_pos is_neg
114 accuracy precision div_scale round_mode fabs fnot
115 objectify upgrade downgrade
120 sub _method_alias { exists $methods{$_[0]||''}; }
121 sub _method_hand_up { exists $hand_ups{$_[0]||''}; }
124 ##############################################################################
129 # create a new BigFloat object from a string or another bigfloat object.
132 # sign => sign (+/-), or "NaN"
134 my ($class,$wanted,@r) = @_;
136 # avoid numify-calls by not using || on $wanted!
137 return $class->bzero() if !defined $wanted; # default to 0
138 return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
140 $class->import() if $IMPORT == 0; # make require work
142 my $self = {}; bless $self, $class;
143 # shortcut for bigints and its subclasses
144 if ((ref($wanted)) && UNIVERSAL::can( $wanted, "as_number"))
146 $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy
147 $self->{_e} = $MBI->_zero();
149 $self->{sign} = $wanted->sign();
150 return $self->bnorm();
152 # else: got a string or something masquerading as number (with overload)
154 # handle '+inf', '-inf' first
155 if ($wanted =~ /^[+-]?inf\z/)
157 return $downgrade->new($wanted) if $downgrade;
159 $self->{sign} = $wanted; # set a default sign for bstr()
160 return $self->binf($wanted);
163 # shortcut for simple forms like '12' that neither have trailing nor leading
165 if ($wanted =~ /^([+-]?)([1-9][0-9]*[1-9])$/)
167 $self->{_e} = $MBI->_zero();
169 $self->{sign} = $1 || '+';
170 $self->{_m} = $MBI->_new($2);
171 return $self->round(@r) if !$downgrade;
174 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
180 Carp::croak ("$wanted is not a number initialized to $class");
183 return $downgrade->bnan() if $downgrade;
185 $self->{_e} = $MBI->_zero();
187 $self->{_m} = $MBI->_zero();
188 $self->{sign} = $nan;
192 # make integer from mantissa by adjusting exp, then convert to int
193 $self->{_e} = $MBI->_new($$ev); # exponent
194 $self->{_es} = $$es || '+';
195 my $mantissa = "$$miv$$mfv"; # create mant.
196 $mantissa =~ s/^0+(\d)/$1/; # strip leading zeros
197 $self->{_m} = $MBI->_new($mantissa); # create mant.
199 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
200 if (CORE::length($$mfv) != 0)
202 my $len = $MBI->_new( CORE::length($$mfv));
203 ($self->{_e}, $self->{_es}) =
204 _e_sub ($self->{_e}, $len, $self->{_es}, '+');
206 # we can only have trailing zeros on the mantissa if $$mfv eq ''
209 # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
210 # because that is faster, especially when _m is not stored in base 10.
211 my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/;
214 my $z = $MBI->_new($zeros);
215 # turn '120e2' into '12e3'
216 $MBI->_rsft ( $self->{_m}, $z, 10);
217 ($self->{_e}, $self->{_es}) =
218 _e_add ( $self->{_e}, $z, $self->{_es}, '+');
221 $self->{sign} = $$mis;
223 # for something like 0Ey, set y to 1, and -0 => +0
224 # Check $$miv for being '0' and $$mfv eq '', because otherwise _m could not
225 # have become 0. That's faster than to call $MBI->_is_zero().
226 $self->{sign} = '+', $self->{_e} = $MBI->_one()
227 if $$miv eq '0' and $$mfv eq '';
229 return $self->round(@r) if !$downgrade;
231 # if downgrade, inf, NaN or integers go down
233 if ($downgrade && $self->{_es} eq '+')
235 if ($MBI->_is_zero( $self->{_e} ))
237 return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
239 return $downgrade->new($self->bsstr());
241 $self->bnorm()->round(@r); # first normalize, then round
246 # if two arguments, the first one is the class to "swallow" subclasses
250 sign => $_[1]->{sign},
252 _m => $MBI->_copy($_[1]->{_m}),
253 _e => $MBI->_copy($_[1]->{_e}),
256 $self->{_a} = $_[1]->{_a} if defined $_[1]->{_a};
257 $self->{_p} = $_[1]->{_p} if defined $_[1]->{_p};
262 sign => $_[0]->{sign},
264 _m => $MBI->_copy($_[0]->{_m}),
265 _e => $MBI->_copy($_[0]->{_e}),
268 $self->{_a} = $_[0]->{_a} if defined $_[0]->{_a};
269 $self->{_p} = $_[0]->{_p} if defined $_[0]->{_p};
275 # used by parent class bone() to initialize number to NaN
281 my $class = ref($self);
282 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
285 $IMPORT=1; # call our import only once
286 $self->{_m} = $MBI->_zero();
287 $self->{_e} = $MBI->_zero();
293 # used by parent class bone() to initialize number to +-inf
299 my $class = ref($self);
300 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
303 $IMPORT=1; # call our import only once
304 $self->{_m} = $MBI->_zero();
305 $self->{_e} = $MBI->_zero();
311 # used by parent class bone() to initialize number to 1
313 $IMPORT=1; # call our import only once
314 $self->{_m} = $MBI->_one();
315 $self->{_e} = $MBI->_zero();
321 # used by parent class bone() to initialize number to 0
323 $IMPORT=1; # call our import only once
324 $self->{_m} = $MBI->_zero();
325 $self->{_e} = $MBI->_one();
331 my ($self,$class) = @_;
332 return if $class =~ /^Math::BigInt/; # we aren't one of these
333 UNIVERSAL::isa($self,$class);
338 # return (later set?) configuration data as hash ref
339 my $class = shift || 'Math::BigFloat';
341 if (@_ == 1 && ref($_[0]) ne 'HASH')
343 my $cfg = $class->SUPER::config();
344 return $cfg->{$_[0]};
347 my $cfg = $class->SUPER::config(@_);
349 # now we need only to override the ones that are different from our parent
350 $cfg->{class} = $class;
355 ##############################################################################
360 # (ref to BFLOAT or num_str ) return num_str
361 # Convert number from internal format to (non-scientific) string format.
362 # internal format is always normalized (no leading zeros, "-0" => "+0")
363 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
365 if ($x->{sign} !~ /^[+-]$/)
367 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
371 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
374 my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
377 $es = $MBI->_str($x->{_m});
378 $len = CORE::length($es);
379 my $e = $MBI->_num($x->{_e});
380 $e = -$e if $x->{_es} eq '-';
384 # if _e is bigger than a scalar, the following will blow your memory
387 my $r = abs($e) - $len;
388 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
392 substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
393 $cad = -$cad if $x->{_es} eq '-';
399 $es .= '0' x $e; $len += $e; $cad = 0;
403 $es = '-'.$es if $x->{sign} eq '-';
404 # if set accuracy or precision, pad with zeros on the right side
405 if ((defined $x->{_a}) && ($not_zero))
407 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
408 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
409 $zeros = $x->{_a} - $len if $cad != $len;
410 $es .= $dot.'0' x $zeros if $zeros > 0;
412 elsif ((($x->{_p} || 0) < 0))
414 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
415 my $zeros = -$x->{_p} + $cad;
416 $es .= $dot.'0' x $zeros if $zeros > 0;
423 # (ref to BFLOAT or num_str ) return num_str
424 # Convert number from internal format to scientific string format.
425 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
426 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
428 if ($x->{sign} !~ /^[+-]$/)
430 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
433 my $sep = 'e'.$x->{_es};
434 my $sign = $x->{sign}; $sign = '' if $sign eq '+';
435 $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
440 # Make a number from a BigFloat object
441 # simple return a string and let Perl's atoi()/atof() handle the rest
442 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
446 ##############################################################################
447 # public stuff (usually prefixed with "b")
451 # (BINT or num_str) return BINT
452 # negate number or make a negated number from string
453 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
455 return $x if $x->modify('bneg');
457 # for +0 dont negate (to have always normalized +0). Does nothing for 'NaN'
458 $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
463 # XXX TODO this must be overwritten and return NaN for non-integer values
464 # band(), bior(), bxor(), too
467 # $class->SUPER::bnot($class,@_);
472 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
475 my ($self,$x,$y) = (ref($_[0]),@_);
476 # objectify is costly, so avoid it
477 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
479 ($self,$x,$y) = objectify(2,@_);
482 return $upgrade->bcmp($x,$y) if defined $upgrade &&
483 ((!$x->isa($self)) || (!$y->isa($self)));
485 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
487 # handle +-inf and NaN
488 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
489 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
490 return +1 if $x->{sign} eq '+inf';
491 return -1 if $x->{sign} eq '-inf';
492 return -1 if $y->{sign} eq '+inf';
496 # check sign for speed first
497 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
498 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
501 my $xz = $x->is_zero();
502 my $yz = $y->is_zero();
503 return 0 if $xz && $yz; # 0 <=> 0
504 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
505 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
507 # adjust so that exponents are equal
508 my $lxm = $MBI->_len($x->{_m});
509 my $lym = $MBI->_len($y->{_m});
510 # the numify somewhat limits our length, but makes it much faster
511 my ($xes,$yes) = (1,1);
512 $xes = -1 if $x->{_es} ne '+';
513 $yes = -1 if $y->{_es} ne '+';
514 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
515 my $ly = $lym + $yes * $MBI->_num($y->{_e});
516 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
517 return $l <=> 0 if $l != 0;
519 # lengths (corrected by exponent) are equal
520 # so make mantissa equal length by padding with zero (shift left)
521 my $diff = $lxm - $lym;
522 my $xm = $x->{_m}; # not yet copy it
526 $ym = $MBI->_copy($y->{_m});
527 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
531 $xm = $MBI->_copy($x->{_m});
532 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
534 my $rc = $MBI->_acmp($xm,$ym);
535 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
541 # Compares 2 values, ignoring their signs.
542 # Returns one of undef, <0, =0, >0. (suitable for sort)
545 my ($self,$x,$y) = (ref($_[0]),@_);
546 # objectify is costly, so avoid it
547 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
549 ($self,$x,$y) = objectify(2,@_);
552 return $upgrade->bacmp($x,$y) if defined $upgrade &&
553 ((!$x->isa($self)) || (!$y->isa($self)));
555 # handle +-inf and NaN's
556 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
558 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
559 return 0 if ($x->is_inf() && $y->is_inf());
560 return 1 if ($x->is_inf() && !$y->is_inf());
565 my $xz = $x->is_zero();
566 my $yz = $y->is_zero();
567 return 0 if $xz && $yz; # 0 <=> 0
568 return -1 if $xz && !$yz; # 0 <=> +y
569 return 1 if $yz && !$xz; # +x <=> 0
571 # adjust so that exponents are equal
572 my $lxm = $MBI->_len($x->{_m});
573 my $lym = $MBI->_len($y->{_m});
574 my ($xes,$yes) = (1,1);
575 $xes = -1 if $x->{_es} ne '+';
576 $yes = -1 if $y->{_es} ne '+';
577 # the numify somewhat limits our length, but makes it much faster
578 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
579 my $ly = $lym + $yes * $MBI->_num($y->{_e});
581 return $l <=> 0 if $l != 0;
583 # lengths (corrected by exponent) are equal
584 # so make mantissa equal-length by padding with zero (shift left)
585 my $diff = $lxm - $lym;
586 my $xm = $x->{_m}; # not yet copy it
590 $ym = $MBI->_copy($y->{_m});
591 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
595 $xm = $MBI->_copy($x->{_m});
596 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
598 $MBI->_acmp($xm,$ym);
603 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
604 # return result as BFLOAT
607 my ($self,$x,$y,@r) = (ref($_[0]),@_);
608 # objectify is costly, so avoid it
609 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
611 ($self,$x,$y,@r) = objectify(2,@_);
614 return $x if $x->modify('badd');
616 # inf and NaN handling
617 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
620 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
622 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
624 # +inf++inf or -inf+-inf => same, rest is NaN
625 return $x if $x->{sign} eq $y->{sign};
628 # +-inf + something => +inf; something +-inf => +-inf
629 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
633 return $upgrade->badd($x,$y,@r) if defined $upgrade &&
634 ((!$x->isa($self)) || (!$y->isa($self)));
636 $r[3] = $y; # no push!
638 # speed: no add for 0+y or x+0
639 return $x->bround(@r) if $y->is_zero(); # x+0
640 if ($x->is_zero()) # 0+y
642 # make copy, clobbering up x (modify in place!)
643 $x->{_e} = $MBI->_copy($y->{_e});
644 $x->{_es} = $y->{_es};
645 $x->{_m} = $MBI->_copy($y->{_m});
646 $x->{sign} = $y->{sign} || $nan;
647 return $x->round(@r);
650 # take lower of the two e's and adapt m1 to it to match m2
652 $e = $MBI->_zero() if !defined $e; # if no BFLOAT?
653 $e = $MBI->_copy($e); # make copy (didn't do it yet)
657 ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
659 my $add = $MBI->_copy($y->{_m});
661 if ($es eq '-') # < 0
663 $MBI->_lsft( $x->{_m}, $e, 10);
664 ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
666 elsif (!$MBI->_is_zero($e)) # > 0
668 $MBI->_lsft($add, $e, 10);
670 # else: both e are the same, so just leave them
672 if ($x->{sign} eq $y->{sign})
675 $x->{_m} = $MBI->_add($x->{_m}, $add);
679 ($x->{_m}, $x->{sign}) =
680 _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
683 # delete trailing zeros, then round
684 $x->bnorm()->round(@r);
687 # sub bsub is inherited from Math::BigInt!
691 # increment arg by one
692 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
694 return $x if $x->modify('binc');
696 if ($x->{_es} eq '-')
698 return $x->badd($self->bone(),@r); # digits after dot
701 if (!$MBI->_is_zero($x->{_e})) # _e == 0 for NaN, inf, -inf
703 # 1e2 => 100, so after the shift below _m has a '0' as last digit
704 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
705 $x->{_e} = $MBI->_zero(); # normalize
707 # we know that the last digit of $x will be '1' or '9', depending on the
711 if ($x->{sign} eq '+')
713 $MBI->_inc($x->{_m});
714 return $x->bnorm()->bround(@r);
716 elsif ($x->{sign} eq '-')
718 $MBI->_dec($x->{_m});
719 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
720 return $x->bnorm()->bround(@r);
722 # inf, nan handling etc
723 $x->badd($self->bone(),@r); # badd() does round
728 # decrement arg by one
729 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
731 return $x if $x->modify('bdec');
733 if ($x->{_es} eq '-')
735 return $x->badd($self->bone('-'),@r); # digits after dot
738 if (!$MBI->_is_zero($x->{_e}))
740 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
741 $x->{_e} = $MBI->_zero(); # normalize
745 my $zero = $x->is_zero();
747 if (($x->{sign} eq '-') || $zero)
749 $MBI->_inc($x->{_m});
750 $x->{sign} = '-' if $zero; # 0 => 1 => -1
751 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
752 return $x->bnorm()->round(@r);
755 elsif ($x->{sign} eq '+')
757 $MBI->_dec($x->{_m});
758 return $x->bnorm()->round(@r);
760 # inf, nan handling etc
761 $x->badd($self->bone('-'),@r); # does round
768 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
770 return $x if $x->modify('blog');
772 # $base > 0, $base != 1; if $base == undef default to $base == e
775 # we need to limit the accuracy to protect against overflow
778 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
780 # also takes care of the "error in _find_round_parameters?" case
781 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
783 # no rounding at all, so must use fallback
784 if (scalar @params == 0)
786 # simulate old behaviour
787 $params[0] = $self->div_scale(); # and round to it as accuracy
788 $params[1] = undef; # P = undef
789 $scale = $params[0]+4; # at least four more for proper round
790 $params[2] = $r; # round mode by caller or undef
791 $fallback = 1; # to clear a/p afterwards
795 # the 4 below is empirical, and there might be cases where it is not
797 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
800 return $x->bzero(@params) if $x->is_one();
801 # base not defined => base == Euler's number e
804 # make object, since we don't feed it through objectify() to still get the
805 # case of $base == undef
806 $base = $self->new($base) unless ref($base);
807 # $base > 0; $base != 1
808 return $x->bnan() if $base->is_zero() || $base->is_one() ||
809 $base->{sign} ne '+';
810 # if $x == $base, we know the result must be 1.0
811 if ($x->bcmp($base) == 0)
813 $x->bone('+',@params);
816 # clear a/p after round, since user did not request it
817 delete $x->{_a}; delete $x->{_p};
823 # when user set globals, they would interfere with our calculation, so
824 # disable them and later re-enable them
826 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
827 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
828 # we also need to disable any set A or P on $x (_find_round_parameters took
829 # them already into account), since these would interfere, too
830 delete $x->{_a}; delete $x->{_p};
831 # need to disable $upgrade in BigInt, to avoid deep recursion
832 local $Math::BigInt::upgrade = undef;
833 local $Math::BigFloat::downgrade = undef;
835 # upgrade $x if $x is not a BigFloat (handle BigInt input)
837 if (!$x->isa('Math::BigFloat'))
839 $x = Math::BigFloat->new($x);
845 # If the base is defined and an integer, try to calculate integer result
846 # first. This is very fast, and in case the real result was found, we can
848 if (defined $base && $base->is_int() && $x->is_int())
850 my $i = $MBI->_copy( $x->{_m} );
851 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
852 my $int = Math::BigInt->bzero();
854 $int->blog($base->as_number());
856 if ($base->as_number()->bpow($int) == $x)
858 # found result, return it
859 $x->{_m} = $int->{value};
860 $x->{_e} = $MBI->_zero();
869 # base is undef, so base should be e (Euler's number), so first calculate the
870 # log to base e (using reduction by 10 (and probably 2)):
871 $self->_log_10($x,$scale);
873 # and if a different base was requested, convert it
876 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
877 # not ln, but some other base (don't modify $base)
878 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
882 # shortcut to not run through _find_round_parameters again
883 if (defined $params[0])
885 $x->bround($params[0],$params[2]); # then round accordingly
889 $x->bfround($params[1],$params[2]); # then round accordingly
893 # clear a/p after round, since user did not request it
894 delete $x->{_a}; delete $x->{_p};
897 $$abr = $ab; $$pbr = $pb;
904 # Given D (digits in decimal), compute N so that N! (N factorial) is
905 # at least D digits long. D should be at least 50.
908 # two constants for the Ramanujan estimate of ln(N!)
909 my $lg2 = log(2 * 3.14159265) / 2;
912 # D = 50 => N => 42, so L = 40 and R = 50
913 my $l = 40; my $r = $d;
915 # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
916 $l = $l->numify if ref($l);
917 $r = $r->numify if ref($r);
918 $lg2 = $lg2->numify if ref($lg2);
919 $lg10 = $lg10->numify if ref($lg10);
921 # binary search for the right value (could this be written as the reverse of lg(n!)?)
924 my $n = int(($r - $l) / 2) + $l;
926 int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
927 $ramanujan > $d ? $r = $n : $l = $n;
934 # Calculate n over k (binomial coefficient or "choose" function) as integer.
936 my ($self,$x,$y,@r) = (ref($_[0]),@_);
938 # objectify is costly, so avoid it
939 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
941 ($self,$x,$y,@r) = objectify(2,@_);
944 return $x if $x->modify('bnok');
946 return $x->bnan() if $x->is_nan() || $y->is_nan();
947 return $x->binf() if $x->is_inf();
949 my $u = $x->as_int();
950 $u->bnok($y->as_int());
952 $x->{_m} = $u->{value};
953 $x->{_e} = $MBI->_zero();
961 # Calculate e ** X (Euler's number to the power of X)
962 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
964 return $x if $x->modify('bexp');
966 return $x->binf() if $x->{sign} eq '+inf';
967 return $x->bzero() if $x->{sign} eq '-inf';
969 # we need to limit the accuracy to protect against overflow
972 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
974 # also takes care of the "error in _find_round_parameters?" case
975 return $x if $x->{sign} eq 'NaN';
977 # no rounding at all, so must use fallback
978 if (scalar @params == 0)
980 # simulate old behaviour
981 $params[0] = $self->div_scale(); # and round to it as accuracy
982 $params[1] = undef; # P = undef
983 $scale = $params[0]+4; # at least four more for proper round
984 $params[2] = $r; # round mode by caller or undef
985 $fallback = 1; # to clear a/p afterwards
989 # the 4 below is empirical, and there might be cases where it's not enough...
990 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
993 return $x->bone(@params) if $x->is_zero();
995 if (!$x->isa('Math::BigFloat'))
997 $x = Math::BigFloat->new($x);
1001 # when user set globals, they would interfere with our calculation, so
1002 # disable them and later re-enable them
1004 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1005 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1006 # we also need to disable any set A or P on $x (_find_round_parameters took
1007 # them already into account), since these would interfere, too
1008 delete $x->{_a}; delete $x->{_p};
1009 # need to disable $upgrade in BigInt, to avoid deep recursion
1010 local $Math::BigInt::upgrade = undef;
1011 local $Math::BigFloat::downgrade = undef;
1013 my $x_org = $x->copy();
1015 # We use the following Taylor series:
1018 # e = 1 + --- + --- + --- + --- ...
1021 # The difference for each term is X and N, which would result in:
1022 # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term
1024 # But it is faster to compute exp(1) and then raising it to the
1025 # given power, esp. if $x is really big and an integer because:
1027 # * The numerator is always 1, making the computation faster
1028 # * the series converges faster in the case of x == 1
1029 # * We can also easily check when we have reached our limit: when the
1030 # term to be added is smaller than "1E$scale", we can stop - f.i.
1031 # scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
1032 # * we can compute the *exact* result by simulating bigrat math:
1034 # 1 1 gcd(3,4) = 1 1*24 + 1*6 5
1035 # - + - = ---------- = --
1038 # We do not compute the gcd() here, but simple do:
1040 # - + - = --------- = --
1044 # a c a*d + c*b and note that c is always 1 and d = (b*f)
1048 # This leads to: which can be reduced by b to:
1049 # a 1 a*b*f + b a*f + 1
1050 # - + - = --------- = -------
1053 # The first terms in the series are:
1055 # 1 1 1 1 1 1 1 1 13700
1056 # -- + -- + -- + -- + -- + --- + --- + ---- = -----
1057 # 1 1 2 6 24 120 720 5040 5040
1059 # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!
1063 # set $x directly from a cached string form
1064 $x->{_m} = $MBI->_new(
1065 "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
1068 $x->{_e} = $MBI->_new(79);
1072 # compute A and B so that e = A / B.
1074 # After some terms we end up with this, so we use it as a starting point:
1075 my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
1076 my $F = $MBI->_new(42); my $step = 42;
1078 # Compute how many steps we need to take to get $A and $B sufficiently big
1079 my $steps = _len_to_steps($scale - 4);
1080 # print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
1081 while ($step++ <= $steps)
1083 # calculate $a * $f + 1
1084 $A = $MBI->_mul($A, $F);
1085 $A = $MBI->_inc($A);
1087 $F = $MBI->_inc($F);
1089 # compute $B as factorial of $steps (this is faster than doing it manually)
1090 my $B = $MBI->_fac($MBI->_new($steps));
1092 # print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
1094 # compute A/B with $scale digits in the result (truncate, not round)
1095 $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
1096 $A = $MBI->_div( $A, $B );
1101 $x->{_e} = $MBI->_new($scale);
1104 # $x contains now an estimate of e, with some surplus digits, so we can round
1105 if (!$x_org->is_one())
1107 # raise $x to the wanted power and round it in one step:
1108 $x->bpow($x_org, @params);
1112 # else just round the already computed result
1113 delete $x->{_a}; delete $x->{_p};
1114 # shortcut to not run through _find_round_parameters again
1115 if (defined $params[0])
1117 $x->bround($params[0],$params[2]); # then round accordingly
1121 $x->bfround($params[1],$params[2]); # then round accordingly
1126 # clear a/p after round, since user did not request it
1127 delete $x->{_a}; delete $x->{_p};
1130 $$abr = $ab; $$pbr = $pb;
1132 $x; # return modified $x
1137 # internal log function to calculate ln() based on Taylor series.
1138 # Modifies $x in place.
1139 my ($self,$x,$scale) = @_;
1141 # in case of $x == 1, result is 0
1142 return $x->bzero() if $x->is_one();
1144 # XXX TODO: rewrite this in a similar manner to bexp()
1146 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
1150 # Taylor: | u 1 u^3 1 u^5 |
1151 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
1152 # |_ v 3 v^3 5 v^5 _|
1154 # This takes much more steps to calculate the result and is thus not used
1157 # Taylor: | u 1 u^2 1 u^3 |
1158 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
1159 # |_ x 2 x^2 3 x^3 _|
1161 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
1163 $v = $x->copy(); $v->binc(); # v = x+1
1164 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
1165 $x->bdiv($v,$scale); # first term: u/v
1166 $below = $v->copy();
1168 $u *= $u; $v *= $v; # u^2, v^2
1169 $below->bmul($v); # u^3, v^3
1171 $factor = $self->new(3); $f = $self->new(2);
1173 my $steps = 0 if DEBUG;
1174 $limit = $self->new("1E-". ($scale-1));
1177 # we calculate the next term, and add it to the last
1178 # when the next term is below our limit, it won't affect the outcome
1179 # anymore, so we stop
1181 # calculating the next term simple from over/below will result in quite
1182 # a time hog if the input has many digits, since over and below will
1183 # accumulate more and more digits, and the result will also have many
1184 # digits, but in the end it is rounded to $scale digits anyway. So if we
1185 # round $over and $below first, we save a lot of time for the division
1186 # (not with log(1.2345), but try log (123**123) to see what I mean. This
1187 # can introduce a rounding error if the division result would be f.i.
1188 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
1189 # if we truncated $over and $below we might get 0.12345. Does this matter
1190 # for the end result? So we give $over and $below 4 more digits to be
1191 # on the safe side (unscientific error handling as usual... :+D
1193 $next = $over->copy->bround($scale+4)->bdiv(
1194 $below->copy->bmul($factor)->bround($scale+4),
1198 ## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
1200 last if $next->bacmp($limit) <= 0;
1202 delete $next->{_a}; delete $next->{_p};
1204 # calculate things for the next term
1205 $over *= $u; $below *= $v; $factor->badd($f);
1208 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
1211 print "took $steps steps\n" if DEBUG;
1212 $x->bmul($f); # $x *= 2
1217 # Internal log function based on reducing input to the range of 0.1 .. 9.99
1218 # and then "correcting" the result to the proper one. Modifies $x in place.
1219 my ($self,$x,$scale) = @_;
1221 # Taking blog() from numbers greater than 10 takes a *very long* time, so we
1222 # break the computation down into parts based on the observation that:
1223 # blog(X*Y) = blog(X) + blog(Y)
1224 # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
1225 # $x is the faster it gets. Since 2*$x takes about 10 times as
1226 # long, we make it faster by about a factor of 100 by dividing $x by 10.
1228 # The same observation is valid for numbers smaller than 0.1, e.g. computing
1229 # log(1) is fastest, and the further away we get from 1, the longer it takes.
1230 # So we also 'break' this down by multiplying $x with 10 and subtract the
1231 # log(10) afterwards to get the correct result.
1233 # To get $x even closer to 1, we also divide by 2 and then use log(2) to
1234 # correct for this. For instance if $x is 2.4, we use the formula:
1235 # blog(2.4 * 2) == blog (1.2) + blog(2)
1236 # and thus calculate only blog(1.2) and blog(2), which is faster in total
1237 # than calculating blog(2.4).
1239 # In addition, the values for blog(2) and blog(10) are cached.
1241 # Calculate nr of digits before dot:
1242 my $dbd = $MBI->_num($x->{_e});
1243 $dbd = -$dbd if $x->{_es} eq '-';
1244 $dbd += $MBI->_len($x->{_m});
1246 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
1247 # infinite recursion
1249 my $calc = 1; # do some calculation?
1251 # disable the shortcut for 10, since we need log(10) and this would recurse
1253 if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
1255 $dbd = 0; # disable shortcut
1256 # we can use the cached value in these cases
1257 if ($scale <= $LOG_10_A)
1259 $x->bzero(); $x->badd($LOG_10); # modify $x in place
1260 $calc = 0; # no need to calc, but round
1262 # if we can't use the shortcut, we continue normally
1266 # disable the shortcut for 2, since we maybe have it cached
1267 if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
1269 $dbd = 0; # disable shortcut
1270 # we can use the cached value in these cases
1271 if ($scale <= $LOG_2_A)
1273 $x->bzero(); $x->badd($LOG_2); # modify $x in place
1274 $calc = 0; # no need to calc, but round
1276 # if we can't use the shortcut, we continue normally
1280 # if $x = 0.1, we know the result must be 0-log(10)
1281 if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
1282 $MBI->_is_one($x->{_m}))
1284 $dbd = 0; # disable shortcut
1285 # we can use the cached value in these cases
1286 if ($scale <= $LOG_10_A)
1288 $x->bzero(); $x->bsub($LOG_10);
1289 $calc = 0; # no need to calc, but round
1293 return if $calc == 0; # already have the result
1295 # default: these correction factors are undef and thus not used
1296 my $l_10; # value of ln(10) to A of $scale
1297 my $l_2; # value of ln(2) to A of $scale
1299 my $two = $self->new(2);
1301 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1302 # so don't do this shortcut for 1 or 0
1303 if (($dbd > 1) || ($dbd < 0))
1305 # convert our cached value to an object if not already (avoid doing this
1306 # at import() time, since not everybody needs this)
1307 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1309 #print "x = $x, dbd = $dbd, calc = $calc\n";
1310 # got more than one digit before the dot, or more than one zero after the
1312 # log(123) == log(1.23) + log(10) * 2
1313 # log(0.0123) == log(1.23) - log(10) * 2
1315 if ($scale <= $LOG_10_A)
1318 $l_10 = $LOG_10->copy(); # copy for mul
1322 # else: slower, compute and cache result
1323 # also disable downgrade for this code path
1324 local $Math::BigFloat::downgrade = undef;
1326 # shorten the time to calculate log(10) based on the following:
1327 # log(1.25 * 8) = log(1.25) + log(8)
1328 # = log(1.25) + log(2) + log(2) + log(2)
1330 # first get $l_2 (and possible compute and cache log(2))
1331 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1332 if ($scale <= $LOG_2_A)
1335 $l_2 = $LOG_2->copy(); # copy() for the mul below
1339 # else: slower, compute and cache result
1340 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1341 $LOG_2 = $l_2->copy(); # cache the result for later
1342 # the copy() is for mul below
1346 # now calculate log(1.25):
1347 $l_10 = $self->new('1.25'); $self->_log($l_10, $scale); # scale+4, actually
1349 # log(1.25) + log(2) + log(2) + log(2):
1353 $LOG_10 = $l_10->copy(); # cache the result for later
1354 # the copy() is for mul below
1357 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
1358 $l_10->bmul( $self->new($dbd)); # log(10) * (digits_before_dot-1)
1365 ($x->{_e}, $x->{_es}) =
1366 _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
1370 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1372 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1373 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1375 $HALF = $self->new($HALF) unless ref($HALF);
1377 my $twos = 0; # default: none (0 times)
1378 while ($x->bacmp($HALF) <= 0) # X <= 0.5
1380 $twos--; $x->bmul($two);
1382 while ($x->bacmp($two) >= 0) # X >= 2
1384 $twos++; $x->bdiv($two,$scale+4); # keep all digits
1386 # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
1387 # So calculate correction factor based on ln(2):
1390 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1391 if ($scale <= $LOG_2_A)
1394 $l_2 = $LOG_2->copy(); # copy() for the mul below
1398 # else: slower, compute and cache result
1399 # also disable downgrade for this code path
1400 local $Math::BigFloat::downgrade = undef;
1401 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1402 $LOG_2 = $l_2->copy(); # cache the result for later
1403 # the copy() is for mul below
1406 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
1409 $self->_log($x,$scale); # need to do the "normal" way
1410 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
1411 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
1413 # all done, $x contains now the result
1419 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1420 # does not modify arguments, but returns new object
1421 # Lowest Common Multiplicator
1423 my ($self,@arg) = objectify(0,@_);
1424 my $x = $self->new(shift @arg);
1425 while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); }
1431 # (BINT or num_str, BINT or num_str) return BINT
1432 # does not modify arguments, but returns new object
1435 $y = __PACKAGE__->new($y) if !ref($y);
1437 my $x = $y->copy()->babs(); # keep arguments
1439 return $x->bnan() if $x->{sign} !~ /^[+-]$/ # x NaN?
1440 || !$x->is_int(); # only for integers now
1444 my $t = shift; $t = $self->new($t) if !ref($t);
1445 $y = $t->copy()->babs();
1447 return $x->bnan() if $y->{sign} !~ /^[+-]$/ # y NaN?
1448 || !$y->is_int(); # only for integers now
1450 # greatest common divisor
1451 while (! $y->is_zero())
1453 ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1456 last if $x->is_one();
1461 ##############################################################################
1465 # Internal helper sub to take two positive integers and their signs and
1466 # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1467 # output ($CALC,('+'|'-'))
1468 my ($x,$y,$xs,$ys) = @_;
1470 # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1473 $x = $MBI->_add ($x, $y ); # a+b
1474 # the sign follows $xs
1478 my $a = $MBI->_acmp($x,$y);
1481 $x = $MBI->_sub ($x , $y); # abs sub
1485 $x = $MBI->_zero(); # result is 0
1490 $x = $MBI->_sub ( $y, $x, 1 ); # abs sub
1498 # Internal helper sub to take two positive integers and their signs and
1499 # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1500 # output ($CALC,('+'|'-'))
1501 my ($x,$y,$xs,$ys) = @_;
1505 _e_add($x,$y,$xs,$ys); # call add (does subtract now)
1508 ###############################################################################
1509 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1513 # return true if arg (BFLOAT or num_str) is an integer
1514 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1516 (($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1517 ($x->{_es} eq '+')) ? 1 : 0; # 1e-1 => no integer
1522 # return true if arg (BFLOAT or num_str) is zero
1523 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1525 ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m})) ? 1 : 0;
1530 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1531 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1533 $sign = '+' if !defined $sign || $sign ne '-';
1535 ($x->{sign} eq $sign &&
1536 $MBI->_is_zero($x->{_e}) &&
1537 $MBI->_is_one($x->{_m}) ) ? 1 : 0;
1542 # return true if arg (BFLOAT or num_str) is odd or false if even
1543 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1545 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1546 ($MBI->_is_zero($x->{_e})) &&
1547 ($MBI->_is_odd($x->{_m}))) ? 1 : 0;
1552 # return true if arg (BINT or num_str) is even or false if odd
1553 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1555 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1556 ($x->{_es} eq '+') && # 123.45 isn't
1557 ($MBI->_is_even($x->{_m}))) ? 1 : 0; # but 1200 is
1562 # multiply two numbers
1565 my ($self,$x,$y,@r) = (ref($_[0]),@_);
1566 # objectify is costly, so avoid it
1567 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1569 ($self,$x,$y,@r) = objectify(2,@_);
1572 return $x if $x->modify('bmul');
1574 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1577 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1579 return $x->bnan() if $x->is_zero() || $y->is_zero();
1580 # result will always be +-inf:
1581 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1582 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1583 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1584 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1585 return $x->binf('-');
1588 return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1589 ((!$x->isa($self)) || (!$y->isa($self)));
1591 # aEb * cEd = (a*c)E(b+d)
1592 $MBI->_mul($x->{_m},$y->{_m});
1593 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1595 $r[3] = $y; # no push!
1598 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1599 $x->bnorm->round(@r);
1604 # multiply two numbers and add the third to the result
1607 my ($self,$x,$y,$z,@r) = (ref($_[0]),@_);
1608 # objectify is costly, so avoid it
1609 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1611 ($self,$x,$y,$z,@r) = objectify(3,@_);
1614 return $x if $x->modify('bmuladd');
1616 return $x->bnan() if (($x->{sign} eq $nan) ||
1617 ($y->{sign} eq $nan) ||
1618 ($z->{sign} eq $nan));
1621 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1623 return $x->bnan() if $x->is_zero() || $y->is_zero();
1624 # result will always be +-inf:
1625 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1626 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1627 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1628 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1629 return $x->binf('-');
1632 return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1633 ((!$x->isa($self)) || (!$y->isa($self)));
1635 # aEb * cEd = (a*c)E(b+d)
1636 $MBI->_mul($x->{_m},$y->{_m});
1637 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1639 $r[3] = $y; # no push!
1642 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1644 # z=inf handling (z=NaN handled above)
1645 $x->{sign} = $z->{sign}, return $x if $z->{sign} =~ /^[+-]inf$/;
1647 # take lower of the two e's and adapt m1 to it to match m2
1649 $e = $MBI->_zero() if !defined $e; # if no BFLOAT?
1650 $e = $MBI->_copy($e); # make copy (didn't do it yet)
1654 ($e,$es) = _e_sub($e, $x->{_e}, $z->{_es} || '+', $x->{_es});
1656 my $add = $MBI->_copy($z->{_m});
1658 if ($es eq '-') # < 0
1660 $MBI->_lsft( $x->{_m}, $e, 10);
1661 ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
1663 elsif (!$MBI->_is_zero($e)) # > 0
1665 $MBI->_lsft($add, $e, 10);
1667 # else: both e are the same, so just leave them
1669 if ($x->{sign} eq $z->{sign})
1672 $x->{_m} = $MBI->_add($x->{_m}, $add);
1676 ($x->{_m}, $x->{sign}) =
1677 _e_add($x->{_m}, $add, $x->{sign}, $z->{sign});
1680 # delete trailing zeros, then round
1681 $x->bnorm()->round(@r);
1686 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1687 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1690 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1691 # objectify is costly, so avoid it
1692 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1694 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1697 return $x if $x->modify('bdiv');
1699 return $self->_div_inf($x,$y)
1700 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1702 # x== 0 # also: or y == 1 or y == -1
1703 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1706 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1708 # we need to limit the accuracy to protect against overflow
1710 my (@params,$scale);
1711 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1713 return $x if $x->is_nan(); # error in _find_round_parameters?
1715 # no rounding at all, so must use fallback
1716 if (scalar @params == 0)
1718 # simulate old behaviour
1719 $params[0] = $self->div_scale(); # and round to it as accuracy
1720 $scale = $params[0]+4; # at least four more for proper round
1721 $params[2] = $r; # round mode by caller or undef
1722 $fallback = 1; # to clear a/p afterwards
1726 # the 4 below is empirical, and there might be cases where it is not
1728 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1731 my $rem; $rem = $self->bzero() if wantarray;
1733 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1735 my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
1736 $scale = $lx if $lx > $scale;
1737 $scale = $ly if $ly > $scale;
1738 my $diff = $ly - $lx;
1739 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
1741 # already handled inf/NaN/-inf above:
1743 # check that $y is not 1 nor -1 and cache the result:
1744 my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));
1746 # flipping the sign of $y will also flip the sign of $x for the special
1747 # case of $x->bsub($x); so we can catch it below:
1748 my $xsign = $x->{sign};
1749 $y->{sign} =~ tr/+-/-+/;
1751 if ($xsign ne $x->{sign})
1753 # special case of $x /= $x results in 1
1754 $x->bone(); # "fixes" also sign of $y, since $x is $y
1758 # correct $y's sign again
1759 $y->{sign} =~ tr/+-/-+/;
1760 # continue with normal div code:
1762 # make copy of $x in case of list context for later reminder calculation
1763 if (wantarray && $y_not_one)
1768 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1770 # check for / +-1 ( +/- 1E0)
1773 # promote BigInts and it's subclasses (except when already a BigFloat)
1774 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1776 # calculate the result to $scale digits and then round it
1777 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1778 $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1779 $MBI->_div ($x->{_m},$y->{_m}); # a/c
1781 # correct exponent of $x
1782 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1783 # correct for 10**scale
1784 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1785 $x->bnorm(); # remove trailing 0's
1787 } # ende else $x != $y
1789 # shortcut to not run through _find_round_parameters again
1790 if (defined $params[0])
1792 delete $x->{_a}; # clear before round
1793 $x->bround($params[0],$params[2]); # then round accordingly
1797 delete $x->{_p}; # clear before round
1798 $x->bfround($params[1],$params[2]); # then round accordingly
1802 # clear a/p after round, since user did not request it
1803 delete $x->{_a}; delete $x->{_p};
1810 $rem->bmod($y,@params); # copy already done
1814 # clear a/p after round, since user did not request it
1815 delete $rem->{_a}; delete $rem->{_p};
1824 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1827 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1828 # objectify is costly, so avoid it
1829 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1831 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1834 return $x if $x->modify('bmod');
1836 # handle NaN, inf, -inf
1837 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1839 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1840 $x->{sign} = $re->{sign};
1841 $x->{_e} = $re->{_e};
1842 $x->{_m} = $re->{_m};
1843 return $x->round($a,$p,$r,$y);
1847 return $x->bnan() if $x->is_zero();
1851 return $x->bzero() if $x->is_zero()
1853 # check that $y == +1 or $y == -1:
1854 ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));
1856 my $cmp = $x->bacmp($y); # equal or $x < $y?
1857 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1859 # only $y of the operands negative?
1860 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1862 $x->{sign} = $y->{sign}; # calc sign first
1863 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1865 my $ym = $MBI->_copy($y->{_m});
1868 $MBI->_lsft( $ym, $y->{_e}, 10)
1869 if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1871 # if $y has digits after dot
1872 my $shifty = 0; # correct _e of $x by this
1873 if ($y->{_es} eq '-') # has digits after dot
1875 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1876 $shifty = $MBI->_num($y->{_e}); # no more digits after dot
1877 $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1879 # $ym is now mantissa of $y based on exponent 0
1881 my $shiftx = 0; # correct _e of $x by this
1882 if ($x->{_es} eq '-') # has digits after dot
1884 # 123.4 % 20 => 1234 % 200
1885 $shiftx = $MBI->_num($x->{_e}); # no more digits after dot
1886 $MBI->_lsft($ym, $x->{_e}, 10); # 123 => 1230
1888 # 123e1 % 20 => 1230 % 20
1889 if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1891 $MBI->_lsft( $x->{_m}, $x->{_e},10); # es => '+' here
1894 $x->{_e} = $MBI->_new($shiftx);
1896 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1897 $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1899 # now mantissas are equalized, exponent of $x is adjusted, so calc result
1901 $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1903 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1906 if ($neg != 0) # one of them negative => correct in place
1909 $x->{_m} = $r->{_m};
1910 $x->{_e} = $r->{_e};
1911 $x->{_es} = $r->{_es};
1912 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
1916 $x->round($a,$p,$r,$y); # round and return
1921 # calculate $y'th root of $x
1924 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1925 # objectify is costly, so avoid it
1926 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1928 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1931 return $x if $x->modify('broot');
1933 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1934 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1935 $y->{sign} !~ /^\+$/;
1937 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1939 # we need to limit the accuracy to protect against overflow
1941 my (@params,$scale);
1942 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1944 return $x if $x->is_nan(); # error in _find_round_parameters?
1946 # no rounding at all, so must use fallback
1947 if (scalar @params == 0)
1949 # simulate old behaviour
1950 $params[0] = $self->div_scale(); # and round to it as accuracy
1951 $scale = $params[0]+4; # at least four more for proper round
1952 $params[2] = $r; # iound mode by caller or undef
1953 $fallback = 1; # to clear a/p afterwards
1957 # the 4 below is empirical, and there might be cases where it is not
1959 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1962 # when user set globals, they would interfere with our calculation, so
1963 # disable them and later re-enable them
1965 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1966 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1967 # we also need to disable any set A or P on $x (_find_round_parameters took
1968 # them already into account), since these would interfere, too
1969 delete $x->{_a}; delete $x->{_p};
1970 # need to disable $upgrade in BigInt, to avoid deep recursion
1971 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1973 # remember sign and make $x positive, since -4 ** (1/2) => -2
1974 my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
1977 if ($y->isa('Math::BigFloat'))
1979 $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
1983 $is_two = ($y == 2);
1986 # normal square root if $y == 2:
1989 $x->bsqrt($scale+4);
1991 elsif ($y->is_one('-'))
1994 my $u = $self->bone()->bdiv($x,$scale);
1995 # copy private parts over
1996 $x->{_m} = $u->{_m};
1997 $x->{_e} = $u->{_e};
1998 $x->{_es} = $u->{_es};
2002 # calculate the broot() as integer result first, and if it fits, return
2003 # it rightaway (but only if $x and $y are integer):
2005 my $done = 0; # not yet
2006 if ($y->is_int() && $x->is_int())
2008 my $i = $MBI->_copy( $x->{_m} );
2009 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2010 my $int = Math::BigInt->bzero();
2012 $int->broot($y->as_number());
2014 if ($int->copy()->bpow($y) == $x)
2016 # found result, return it
2017 $x->{_m} = $int->{value};
2018 $x->{_e} = $MBI->_zero();
2026 my $u = $self->bone()->bdiv($y,$scale+4);
2027 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
2028 $x->bpow($u,$scale+4); # el cheapo
2031 $x->bneg() if $sign == 1;
2033 # shortcut to not run through _find_round_parameters again
2034 if (defined $params[0])
2036 $x->bround($params[0],$params[2]); # then round accordingly
2040 $x->bfround($params[1],$params[2]); # then round accordingly
2044 # clear a/p after round, since user did not request it
2045 delete $x->{_a}; delete $x->{_p};
2048 $$abr = $ab; $$pbr = $pb;
2054 # calculate square root
2055 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2057 return $x if $x->modify('bsqrt');
2059 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
2060 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
2061 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
2063 # we need to limit the accuracy to protect against overflow
2065 my (@params,$scale);
2066 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2068 return $x if $x->is_nan(); # error in _find_round_parameters?
2070 # no rounding at all, so must use fallback
2071 if (scalar @params == 0)
2073 # simulate old behaviour
2074 $params[0] = $self->div_scale(); # and round to it as accuracy
2075 $scale = $params[0]+4; # at least four more for proper round
2076 $params[2] = $r; # round mode by caller or undef
2077 $fallback = 1; # to clear a/p afterwards
2081 # the 4 below is empirical, and there might be cases where it is not
2083 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2086 # when user set globals, they would interfere with our calculation, so
2087 # disable them and later re-enable them
2089 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2090 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2091 # we also need to disable any set A or P on $x (_find_round_parameters took
2092 # them already into account), since these would interfere, too
2093 delete $x->{_a}; delete $x->{_p};
2094 # need to disable $upgrade in BigInt, to avoid deep recursion
2095 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2097 my $i = $MBI->_copy( $x->{_m} );
2098 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2099 my $xas = Math::BigInt->bzero();
2102 my $gs = $xas->copy()->bsqrt(); # some guess
2104 if (($x->{_es} ne '-') # guess can't be accurate if there are
2105 # digits after the dot
2106 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
2108 # exact result, copy result over to keep $x
2109 $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
2111 # shortcut to not run through _find_round_parameters again
2112 if (defined $params[0])
2114 $x->bround($params[0],$params[2]); # then round accordingly
2118 $x->bfround($params[1],$params[2]); # then round accordingly
2122 # clear a/p after round, since user did not request it
2123 delete $x->{_a}; delete $x->{_p};
2125 # re-enable A and P, upgrade is taken care of by "local"
2126 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
2130 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
2131 # of the result by multiplying the input by 100 and then divide the integer
2132 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
2134 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
2135 my $y1 = $MBI->_copy($x->{_m});
2137 my $length = $MBI->_len($y1);
2139 # Now calculate how many digits the result of sqrt(y1) would have
2140 my $digits = int($length / 2);
2142 # But we need at least $scale digits, so calculate how many are missing
2143 my $shift = $scale - $digits;
2145 # This happens if the input had enough digits
2146 # (we take care of integer guesses above)
2147 $shift = 0 if $shift < 0;
2149 # Multiply in steps of 100, by shifting left two times the "missing" digits
2150 my $s2 = $shift * 2;
2152 # We now make sure that $y1 has the same odd or even number of digits than
2153 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
2154 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
2155 # steps of 10. The length of $x does not count, since an even or odd number
2156 # of digits before the dot is not changed by adding an even number of digits
2157 # after the dot (the result is still odd or even digits long).
2158 $s2++ if $MBI->_is_odd($x->{_e});
2160 $MBI->_lsft( $y1, $MBI->_new($s2), 10);
2162 # now take the square root and truncate to integer
2163 $y1 = $MBI->_sqrt($y1);
2165 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
2166 # result, which is than later rounded to the desired scale.
2168 # calculate how many zeros $x had after the '.' (or before it, depending
2169 # on sign of $dat, the result should have half as many:
2170 my $dat = $MBI->_num($x->{_e});
2171 $dat = -$dat if $x->{_es} eq '-';
2176 # no zeros after the dot (e.g. 1.23, 0.49 etc)
2177 # preserve half as many digits before the dot than the input had
2178 # (but round this "up")
2179 $dat = int(($dat+1)/2);
2183 $dat = int(($dat)/2);
2185 $dat -= $MBI->_len($y1);
2189 $x->{_e} = $MBI->_new( $dat );
2194 $x->{_e} = $MBI->_new( $dat );
2200 # shortcut to not run through _find_round_parameters again
2201 if (defined $params[0])
2203 $x->bround($params[0],$params[2]); # then round accordingly
2207 $x->bfround($params[1],$params[2]); # then round accordingly
2211 # clear a/p after round, since user did not request it
2212 delete $x->{_a}; delete $x->{_p};
2215 $$abr = $ab; $$pbr = $pb;
2221 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2222 # compute factorial number, modifies first argument
2225 my ($self,$x,@r) = (ref($_[0]),@_);
2226 # objectify is costly, so avoid it
2227 ($self,$x,@r) = objectify(1,@_) if !ref($x);
2230 return $x if $x->modify('bfac') || $x->{sign} eq '+inf';
2233 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
2234 ($x->{_es} ne '+')); # digits after dot?
2236 # use BigInt's bfac() for faster calc
2237 if (! $MBI->_is_zero($x->{_e}))
2239 $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
2240 $x->{_e} = $MBI->_zero(); # normalize
2243 $MBI->_fac($x->{_m}); # calculate factorial
2244 $x->bnorm()->round(@r); # norm again and round result
2249 # Calculate a power where $y is a non-integer, like 2 ** 0.3
2253 # if $y == 0.5, it is sqrt($x)
2254 $HALF = $self->new($HALF) unless ref($HALF);
2255 return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0;
2258 # a ** x == e ** (x * ln a)
2262 # Taylor: | u u^2 u^3 |
2263 # x ** y = 1 + | --- + --- + ----- + ... |
2266 # we need to limit the accuracy to protect against overflow
2268 my ($scale,@params);
2269 ($x,@params) = $x->_find_round_parameters(@r);
2271 return $x if $x->is_nan(); # error in _find_round_parameters?
2273 # no rounding at all, so must use fallback
2274 if (scalar @params == 0)
2276 # simulate old behaviour
2277 $params[0] = $self->div_scale(); # and round to it as accuracy
2278 $params[1] = undef; # disable P
2279 $scale = $params[0]+4; # at least four more for proper round
2280 $params[2] = $r[2]; # round mode by caller or undef
2281 $fallback = 1; # to clear a/p afterwards
2285 # the 4 below is empirical, and there might be cases where it is not
2287 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2290 # when user set globals, they would interfere with our calculation, so
2291 # disable them and later re-enable them
2293 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2294 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2295 # we also need to disable any set A or P on $x (_find_round_parameters took
2296 # them already into account), since these would interfere, too
2297 delete $x->{_a}; delete $x->{_p};
2298 # need to disable $upgrade in BigInt, to avoid deep recursion
2299 local $Math::BigInt::upgrade = undef;
2301 my ($limit,$v,$u,$below,$factor,$next,$over);
2303 $u = $x->copy()->blog(undef,$scale)->bmul($y);
2304 $v = $self->bone(); # 1
2305 $factor = $self->new(2); # 2
2306 $x->bone(); # first term: 1
2308 $below = $v->copy();
2311 $limit = $self->new("1E-". ($scale-1));
2315 # we calculate the next term, and add it to the last
2316 # when the next term is below our limit, it won't affect the outcome
2317 # anymore, so we stop:
2318 $next = $over->copy()->bdiv($below,$scale);
2319 last if $next->bacmp($limit) <= 0;
2321 # calculate things for the next term
2322 $over *= $u; $below *= $factor; $factor->binc();
2324 last if $x->{sign} !~ /^[-+]$/;
2329 # shortcut to not run through _find_round_parameters again
2330 if (defined $params[0])
2332 $x->bround($params[0],$params[2]); # then round accordingly
2336 $x->bfround($params[1],$params[2]); # then round accordingly
2340 # clear a/p after round, since user did not request it
2341 delete $x->{_a}; delete $x->{_p};
2344 $$abr = $ab; $$pbr = $pb;
2350 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2351 # compute power of two numbers, second arg is used as integer
2352 # modifies first argument
2355 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2356 # objectify is costly, so avoid it
2357 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2359 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2362 return $x if $x->modify('bpow');
2364 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
2365 return $x if $x->{sign} =~ /^[+-]inf$/;
2367 # cache the result of is_zero
2368 my $y_is_zero = $y->is_zero();
2369 return $x->bone() if $y_is_zero;
2370 return $x if $x->is_one() || $y->is_one();
2372 my $x_is_zero = $x->is_zero();
2373 return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int(); # non-integer power
2375 my $y1 = $y->as_number()->{value}; # make MBI part
2378 if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
2380 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
2381 return $MBI->_is_odd($y1) ? $x : $x->babs(1);
2385 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
2386 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
2391 $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
2393 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
2394 $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
2395 $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
2397 $x->{sign} = $new_sign;
2399 if ($y->{sign} eq '-')
2401 # modify $x in place!
2402 my $z = $x->copy(); $x->bone();
2403 return scalar $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
2405 $x->round($a,$p,$r,$y);
2410 # takes a very large number to a very large exponent in a given very
2411 # large modulus, quickly, thanks to binary exponentation. Supports
2412 # negative exponents.
2413 my ($self,$num,$exp,$mod,@r) = objectify(3,@_);
2415 return $num if $num->modify('bmodpow');
2417 # check modulus for valid values
2418 return $num->bnan() if ($mod->{sign} ne '+' # NaN, - , -inf, +inf
2419 || $mod->is_zero());
2421 # check exponent for valid values
2422 if ($exp->{sign} =~ /\w/)
2424 # i.e., if it's NaN, +inf, or -inf...
2425 return $num->bnan();
2428 $num->bmodinv ($mod) if ($exp->{sign} eq '-');
2430 # check num for valid values (also NaN if there was no inverse but $exp < 0)
2431 return $num->bnan() if $num->{sign} !~ /^[+-]$/;
2433 # $mod is positive, sign on $exp is ignored, result also positive
2435 # XXX TODO: speed it up when all three numbers are integers
2436 $num->bpow($exp)->bmod($mod);
2439 ###############################################################################
2440 # trigonometric functions
2442 # helper function for bpi() and batan2(), calculates arcus tanges (1/x)
2446 # return a/b so that a/b approximates atan(1/x) to at least limit digits
2447 my ($self, $x, $limit) = @_;
2449 # Taylor: x^3 x^5 x^7 x^9
2450 # atan = x - --- + --- - --- + --- - ...
2454 # atan 1/x = - - ------- + ------- - ------- + ...
2455 # x x^3 * 3 x^5 * 5 x^7 * 7
2458 # atan 1/x = - - --------- + ---------- - ----------- + ...
2459 # 5 3 * 125 5 * 3125 7 * 78125
2461 # Subtraction/addition of a rational:
2464 # - +- - = ----------
2469 # a 1 a * d * c +- b
2470 # ----- +- ------------------ = ----------------
2473 # since b1 = b0 * (d-2) * c
2475 # a 1 a * d +- b / c
2476 # ----- +- ------------------ = ----------------
2483 # stop if length($u) > limit
2490 my $a = $MBI->_one();
2491 my $b = $MBI->_copy($x);
2493 my $x2 = $MBI->_mul( $MBI->_copy($x), $b); # x2 = x * x
2494 my $d = $MBI->_new( 3 ); # d = 3
2495 my $c = $MBI->_mul( $MBI->_copy($x), $x2); # c = x ^ 3
2496 my $two = $MBI->_new( 2 );
2498 # run the first step unconditionally
2499 my $u = $MBI->_mul( $MBI->_copy($d), $c);
2500 $a = $MBI->_mul($a, $u);
2501 $a = $MBI->_sub($a, $b);
2502 $b = $MBI->_mul($b, $u);
2503 $d = $MBI->_add($d, $two);
2504 $c = $MBI->_mul($c, $x2);
2506 # a is now a * (d-3) * c
2507 # b is now b * (d-2) * c
2509 # run the second step unconditionally
2510 $u = $MBI->_mul( $MBI->_copy($d), $c);
2511 $a = $MBI->_mul($a, $u);
2512 $a = $MBI->_add($a, $b);
2513 $b = $MBI->_mul($b, $u);
2514 $d = $MBI->_add($d, $two);
2515 $c = $MBI->_mul($c, $x2);
2517 # a is now a * (d-3) * (d-5) * c * c
2518 # b is now b * (d-2) * (d-4) * c * c
2520 # so we can remove c * c from both a and b to shorten the numbers involved:
2521 $a = $MBI->_div($a, $x2);
2522 $b = $MBI->_div($b, $x2);
2523 $a = $MBI->_div($a, $x2);
2524 $b = $MBI->_div($b, $x2);
2527 my $sign = 0; # 0 => -, 1 => +
2531 # if (($i++ % 100) == 0)
2533 # print "a=",$MBI->_str($a),"\n";
2534 # print "b=",$MBI->_str($b),"\n";
2536 # print "d=",$MBI->_str($d),"\n";
2537 # print "x2=",$MBI->_str($x2),"\n";
2538 # print "c=",$MBI->_str($c),"\n";
2540 my $u = $MBI->_mul( $MBI->_copy($d), $c);
2541 # use _alen() for libs like GMP where _len() would be O(N^2)
2542 last if $MBI->_alen($u) > $limit;
2543 my ($bc,$r) = $MBI->_div( $MBI->_copy($b), $c);
2544 if ($MBI->_is_zero($r))
2546 # b / c is an integer, so we can remove c from all terms
2547 # this happens almost every time:
2548 $a = $MBI->_mul($a, $d);
2549 $a = $MBI->_sub($a, $bc) if $sign == 0;
2550 $a = $MBI->_add($a, $bc) if $sign == 1;
2551 $b = $MBI->_mul($b, $d);
2555 # b / c is not an integer, so we keep c in the terms
2556 # this happens very rarely, for instance for x = 5, this happens only
2557 # at the following steps:
2558 # 1, 5, 14, 32, 72, 157, 340, ...
2559 $a = $MBI->_mul($a, $u);
2560 $a = $MBI->_sub($a, $b) if $sign == 0;
2561 $a = $MBI->_add($a, $b) if $sign == 1;
2562 $b = $MBI->_mul($b, $u);
2564 $d = $MBI->_add($d, $two);
2565 $c = $MBI->_mul($c, $x2);
2570 # print "Took $step steps for ", $MBI->_str($x),"\n";
2571 # print "a=",$MBI->_str($a),"\n"; print "b=",$MBI->_str($b),"\n";
2572 # return a/b so that a/b approximates atan(1/x)
2585 # called like Math::BigFloat::bpi(10);
2586 $n = $self; $self = $class;
2587 # called like Math::BigFloat->bpi();
2588 $n = undef if $n eq 'Math::BigFloat';
2590 $self = ref($self) if ref($self);
2591 my $fallback = defined $n ? 0 : 1;
2592 $n = 40 if !defined $n || $n < 1;
2594 # after 黃見利 (Hwang Chien-Lih) (1997)
2595 # pi/4 = 183 * atan(1/239) + 32 * atan(1/1023) – 68 * atan(1/5832)
2596 # + 12 * atan(1/110443) - 12 * atan(1/4841182) - 100 * atan(1/6826318)
2598 # a few more to prevent rounding errors
2601 my ($a,$b) = $self->_atan_inv( $MBI->_new(239),$n);
2602 my ($c,$d) = $self->_atan_inv( $MBI->_new(1023),$n);
2603 my ($e,$f) = $self->_atan_inv( $MBI->_new(5832),$n);
2604 my ($g,$h) = $self->_atan_inv( $MBI->_new(110443),$n);
2605 my ($i,$j) = $self->_atan_inv( $MBI->_new(4841182),$n);
2606 my ($k,$l) = $self->_atan_inv( $MBI->_new(6826318),$n);
2608 $MBI->_mul($a, $MBI->_new(732));
2609 $MBI->_mul($c, $MBI->_new(128));
2610 $MBI->_mul($e, $MBI->_new(272));
2611 $MBI->_mul($g, $MBI->_new(48));
2612 $MBI->_mul($i, $MBI->_new(48));
2613 $MBI->_mul($k, $MBI->_new(400));
2615 my $x = $self->bone(); $x->{_m} = $a; my $x_d = $self->bone(); $x_d->{_m} = $b;
2616 my $y = $self->bone(); $y->{_m} = $c; my $y_d = $self->bone(); $y_d->{_m} = $d;
2617 my $z = $self->bone(); $z->{_m} = $e; my $z_d = $self->bone(); $z_d->{_m} = $f;
2618 my $u = $self->bone(); $u->{_m} = $g; my $u_d = $self->bone(); $u_d->{_m} = $h;
2619 my $v = $self->bone(); $v->{_m} = $i; my $v_d = $self->bone(); $v_d->{_m} = $j;
2620 my $w = $self->bone(); $w->{_m} = $k; my $w_d = $self->bone(); $w_d->{_m} = $l;
2628 delete $x->{_a}; delete $y->{_a}; delete $z->{_a};
2629 delete $u->{_a}; delete $v->{_a}; delete $w->{_a};
2630 $x->badd($y)->bsub($z)->badd($u)->bsub($v)->bsub($w);
2633 delete $x->{_a} if $fallback == 1;
2639 # Calculate a cosinus of x.
2640 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2642 # Taylor: x^2 x^4 x^6 x^8
2643 # cos = 1 - --- + --- - --- + --- ...
2646 # we need to limit the accuracy to protect against overflow
2648 my ($scale,@params);
2649 ($x,@params) = $x->_find_round_parameters(@r);
2651 # constant object or error in _find_round_parameters?
2652 return $x if $x->modify('bcos') || $x->is_nan();
2654 return $x->bone(@r) if $x->is_zero();
2656 # no rounding at all, so must use fallback
2657 if (scalar @params == 0)
2659 # simulate old behaviour
2660 $params[0] = $self->div_scale(); # and round to it as accuracy
2661 $params[1] = undef; # disable P
2662 $scale = $params[0]+4; # at least four more for proper round
2663 $params[2] = $r[2]; # round mode by caller or undef
2664 $fallback = 1; # to clear a/p afterwards
2668 # the 4 below is empirical, and there might be cases where it is not
2670 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2673 # when user set globals, they would interfere with our calculation, so
2674 # disable them and later re-enable them
2676 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2677 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2678 # we also need to disable any set A or P on $x (_find_round_parameters took
2679 # them already into account), since these would interfere, too
2680 delete $x->{_a}; delete $x->{_p};
2681 # need to disable $upgrade in BigInt, to avoid deep recursion
2682 local $Math::BigInt::upgrade = undef;
2685 my $over = $x * $x; # X ^ 2
2686 my $x2 = $over->copy(); # X ^ 2; difference between terms
2687 my $sign = 1; # start with -=
2688 my $below = $self->new(2); my $factorial = $self->new(3);
2689 $x->bone(); delete $x->{_a}; delete $x->{_p};
2691 my $limit = $self->new("1E-". ($scale-1));
2695 # we calculate the next term, and add it to the last
2696 # when the next term is below our limit, it won't affect the outcome
2697 # anymore, so we stop:
2698 my $next = $over->copy()->bdiv($below,$scale);
2699 last if $next->bacmp($limit) <= 0;
2709 $sign = 1-$sign; # alternate
2710 # calculate things for the next term
2711 $over->bmul($x2); # $x*$x
2712 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2713 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2716 # shortcut to not run through _find_round_parameters again
2717 if (defined $params[0])
2719 $x->bround($params[0],$params[2]); # then round accordingly
2723 $x->bfround($params[1],$params[2]); # then round accordingly
2727 # clear a/p after round, since user did not request it
2728 delete $x->{_a}; delete $x->{_p};
2731 $$abr = $ab; $$pbr = $pb;
2737 # Calculate a sinus of x.
2738 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2740 # taylor: x^3 x^5 x^7 x^9
2741 # sin = x - --- + --- - --- + --- ...
2744 # we need to limit the accuracy to protect against overflow
2746 my ($scale,@params);
2747 ($x,@params) = $x->_find_round_parameters(@r);
2749 # constant object or error in _find_round_parameters?
2750 return $x if $x->modify('bsin') || $x->is_nan();
2752 return $x->bzero(@r) if $x->is_zero();
2754 # no rounding at all, so must use fallback
2755 if (scalar @params == 0)
2757 # simulate old behaviour
2758 $params[0] = $self->div_scale(); # and round to it as accuracy
2759 $params[1] = undef; # disable P
2760 $scale = $params[0]+4; # at least four more for proper round
2761 $params[2] = $r[2]; # round mode by caller or undef
2762 $fallback = 1; # to clear a/p afterwards
2766 # the 4 below is empirical, and there might be cases where it is not
2768 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2771 # when user set globals, they would interfere with our calculation, so
2772 # disable them and later re-enable them
2774 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2775 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2776 # we also need to disable any set A or P on $x (_find_round_parameters took
2777 # them already into account), since these would interfere, too
2778 delete $x->{_a}; delete $x->{_p};
2779 # need to disable $upgrade in BigInt, to avoid deep recursion
2780 local $Math::BigInt::upgrade = undef;
2783 my $over = $x * $x; # X ^ 2
2784 my $x2 = $over->copy(); # X ^ 2; difference between terms
2785 $over->bmul($x); # X ^ 3 as starting value
2786 my $sign = 1; # start with -=
2787 my $below = $self->new(6); my $factorial = $self->new(4);
2788 delete $x->{_a}; delete $x->{_p};
2790 my $limit = $self->new("1E-". ($scale-1));
2794 # we calculate the next term, and add it to the last
2795 # when the next term is below our limit, it won't affect the outcome
2796 # anymore, so we stop:
2797 my $next = $over->copy()->bdiv($below,$scale);
2798 last if $next->bacmp($limit) <= 0;
2808 $sign = 1-$sign; # alternate
2809 # calculate things for the next term
2810 $over->bmul($x2); # $x*$x
2811 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2812 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2815 # shortcut to not run through _find_round_parameters again
2816 if (defined $params[0])
2818 $x->bround($params[0],$params[2]); # then round accordingly
2822 $x->bfround($params[1],$params[2]); # then round accordingly
2826 # clear a/p after round, since user did not request it
2827 delete $x->{_a}; delete $x->{_p};
2830 $$abr = $ab; $$pbr = $pb;
2836 # calculate arcus tangens of ($y/$x)
2839 my ($self,$y,$x,@r) = (ref($_[0]),@_);
2840 # objectify is costly, so avoid it
2841 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2843 ($self,$y,$x,@r) = objectify(2,@_);
2846 return $y if $y->modify('batan2');
2848 return $y->bnan() if ($y->{sign} eq $nan) || ($x->{sign} eq $nan);
2854 return $y->bzero(@r) if ($x->is_inf('+') && !$y->is_inf()) || ($y->is_zero() && $x->{sign} eq '+');
2857 # != 0 -inf result is +- pi
2858 if ($x->is_inf() || $y->is_inf())
2861 my $pi = $self->bpi(@r);
2864 # upgrade to BigRat etc.
2865 return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
2866 if ($x->{sign} eq '-inf')
2869 $MBI->_mul($pi->{_m}, $MBI->_new(3));
2870 $MBI->_div($pi->{_m}, $MBI->_new(4));
2872 elsif ($x->{sign} eq '+inf')
2875 $MBI->_div($pi->{_m}, $MBI->_new(4));
2880 $MBI->_div($pi->{_m}, $MBI->_new(2));
2882 $y->{sign} = substr($y->{sign},0,1); # keep +/-
2884 # modify $y in place
2885 $y->{_m} = $pi->{_m};
2886 $y->{_e} = $pi->{_e};
2887 $y->{_es} = $pi->{_es};
2888 # keep the sign of $y
2892 return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
2899 my $pi = $self->bpi(@r);
2900 # modify $y in place
2901 $y->{_m} = $pi->{_m};
2902 $y->{_e} = $pi->{_e};
2903 $y->{_es} = $pi->{_es};
2909 # +y 0 result is PI/2
2910 # -y 0 result is -PI/2
2914 my $pi = $self->bpi(@r);
2915 # modify $y in place
2916 $y->{_m} = $pi->{_m};
2917 $y->{_e} = $pi->{_e};
2918 $y->{_es} = $pi->{_es};
2919 # -y => -PI/2, +y => PI/2
2920 $MBI->_div($y->{_m}, $MBI->_new(2));
2924 # we need to limit the accuracy to protect against overflow
2926 my ($scale,@params);
2927 ($y,@params) = $y->_find_round_parameters(@r);
2929 # error in _find_round_parameters?
2930 return $y if $y->is_nan();
2932 # no rounding at all, so must use fallback
2933 if (scalar @params == 0)
2935 # simulate old behaviour
2936 $params[0] = $self->div_scale(); # and round to it as accuracy
2937 $params[1] = undef; # disable P
2938 $scale = $params[0]+4; # at least four more for proper round
2939 $params[2] = $r[2]; # round mode by caller or undef
2940 $fallback = 1; # to clear a/p afterwards
2944 # the 4 below is empirical, and there might be cases where it is not
2946 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2949 # inlined is_one() && is_one('-')
2950 if ($MBI->_is_one($y->{_m}) && $MBI->_is_zero($y->{_e}))
2952 # shortcut: 1 1 result is PI/4
2953 # inlined is_one() && is_one('-')
2954 if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
2957 my $pi_4 = $self->bpi( $scale - 3);
2958 # modify $y in place
2959 $y->{_m} = $pi_4->{_m};
2960 $y->{_e} = $pi_4->{_e};
2961 $y->{_es} = $pi_4->{_es};
2966 $y->{sign} = $x->{sign} eq $y->{sign} ? '+' : '-';
2967 $MBI->_div($y->{_m}, $MBI->_new(4));
2970 # shortcut: 1 int(X) result is _atan_inv(X)
2973 if ($x->{_es} eq '+')
2975 my $x1 = $MBI->_copy($x->{_m});
2976 $MBI->_lsft($x1, $x->{_e},10) unless $MBI->_is_zero($x->{_e});
2978 my ($a,$b) = $self->_atan_inv($x1, $scale);
2979 my $y_sign = $y->{sign};
2981 $y->bone(); $y->{_m} = $a; my $y_d = $self->bone(); $y_d->{_m} = $b;
2983 $y->{sign} = $y_sign;
2988 # handle all other cases
2993 # -x -y -PI/2 to -PI
2995 my $y_sign = $y->{sign};
2998 $y->bdiv($x, $scale) unless $x->is_one();
3002 $y->{sign} = $y_sign;
3009 # Calculate a arcus tangens of x.
3013 # taylor: x^3 x^5 x^7 x^9
3014 # atan = x - --- + --- - --- + --- ...
3017 # we need to limit the accuracy to protect against overflow
3019 my ($scale,@params);
3020 ($x,@params) = $x->_find_round_parameters(@r);
3022 # constant object or error in _find_round_parameters?
3023 return $x if $x->modify('batan') || $x->is_nan();
3025 if ($x->{sign} =~ /^[+-]inf\z/)
3027 # +inf result is PI/2
3028 # -inf result is -PI/2
3030 my $pi = $self->bpi(@r);
3031 # modify $x in place
3032 $x->{_m} = $pi->{_m};
3033 $x->{_e} = $pi->{_e};
3034 $x->{_es} = $pi->{_es};
3035 # -y => -PI/2, +y => PI/2
3036 $x->{sign} = substr($x->{sign},0,1); # +inf => +
3037 $MBI->_div($x->{_m}, $MBI->_new(2));
3041 return $x->bzero(@r) if $x->is_zero();
3043 # no rounding at all, so must use fallback
3044 if (scalar @params == 0)
3046 # simulate old behaviour
3047 $params[0] = $self->div_scale(); # and round to it as accuracy
3048 $params[1] = undef; # disable P
3049 $scale = $params[0]+4; # at least four more for proper round
3050 $params[2] = $r[2]; # round mode by caller or undef
3051 $fallback = 1; # to clear a/p afterwards
3055 # the 4 below is empirical, and there might be cases where it is not
3057 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3061 # inlined is_one() && is_one('-')
3062 if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3064 my $pi = $self->bpi($scale - 3);
3065 # modify $x in place
3066 $x->{_m} = $pi->{_m};
3067 $x->{_e} = $pi->{_e};
3068 $x->{_es} = $pi->{_es};
3069 # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4)
3070 $MBI->_div($x->{_m}, $MBI->_new(4));
3074 # This series is only valid if -1 < x < 1, so for other x we need to
3075 # to calculate PI/2 - atan(1/x):
3076 my $one = $MBI->_new(1);
3078 if ($x->{_es} eq '+' && ($MBI->_acmp($x->{_m},$one) >= 0))
3081 $pi = $self->bpi($scale - 3);
3082 $MBI->_div($pi->{_m}, $MBI->_new(2));
3084 my $x_copy = $x->copy();
3085 # modify $x in place
3086 $x->bone(); $x->bdiv($x_copy,$scale);
3089 # when user set globals, they would interfere with our calculation, so
3090 # disable them and later re-enable them
3092 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
3093 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
3094 # we also need to disable any set A or P on $x (_find_round_parameters took
3095 # them already into account), since these would interfere, too
3096 delete $x->{_a}; delete $x->{_p};
3097 # need to disable $upgrade in BigInt, to avoid deep recursion
3098 local $Math::BigInt::upgrade = undef;
3101 my $over = $x * $x; # X ^ 2
3102 my $x2 = $over->copy(); # X ^ 2; difference between terms
3103 $over->bmul($x); # X ^ 3 as starting value
3104 my $sign = 1; # start with -=
3105 my $below = $self->new(3);
3106 my $two = $self->new(2);
3107 delete $x->{_a}; delete $x->{_p};
3109 my $limit = $self->new("1E-". ($scale-1));
3113 # we calculate the next term, and add it to the last
3114 # when the next term is below our limit, it won't affect the outcome
3115 # anymore, so we stop:
3116 my $next = $over->copy()->bdiv($below,$scale);
3117 last if $next->bacmp($limit) <= 0;
3127 $sign = 1-$sign; # alternate
3128 # calculate things for the next term
3129 $over->bmul($x2); # $x*$x
3130 $below->badd($two); # n += 2
3135 my $x_copy = $x->copy();
3136 # modify $x in place
3137 $x->{_m} = $pi->{_m};
3138 $x->{_e} = $pi->{_e};
3139 $x->{_es} = $pi->{_es};
3144 # shortcut to not run through _find_round_parameters again
3145 if (defined $params[0])
3147 $x->bround($params[0],$params[2]); # then round accordingly
3151 $x->bfround($params[1],$params[2]); # then round accordingly
3155 # clear a/p after round, since user did not request it
3156 delete $x->{_a}; delete $x->{_p};
3159 $$abr = $ab; $$pbr = $pb;
3163 ###############################################################################
3164 # rounding functions
3168 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
3169 # $n == 0 means round to integer
3170 # expects and returns normalized numbers!
3171 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
3173 my ($scale,$mode) = $x->_scale_p(@_);
3174 return $x if !defined $scale || $x->modify('bfround'); # no-op
3176 # never round a 0, +-inf, NaN
3179 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
3182 return $x if $x->{sign} !~ /^[+-]$/;
3184 # don't round if x already has lower precision
3185 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
3187 $x->{_p} = $scale; # remember round in any case
3188 delete $x->{_a}; # and clear A
3191 # round right from the '.'
3193 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round
3195 $scale = -$scale; # positive for simplicity
3196 my $len = $MBI->_len($x->{_m}); # length of mantissa
3198 # the following poses a restriction on _e, but if _e is bigger than a
3199 # scalar, you got other problems (memory etc) anyway
3200 my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e}))); # digits after dot
3201 my $zad = 0; # zeros after dot
3202 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
3204 # p rint "scale $scale dad $dad zad $zad len $len\n";
3205 # number bsstr len zad dad
3206 # 0.123 123e-3 3 0 3
3207 # 0.0123 123e-4 3 1 4
3210 # 1.2345 12345e-4 5 0 4
3212 # do not round after/right of the $dad
3213 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
3215 # round to zero if rounding inside the $zad, but not for last zero like:
3216 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
3217 return $x->bzero() if $scale < $zad;
3218 if ($scale == $zad) # for 0.006, scale -3 and trunc
3224 # adjust round-point to be inside mantissa
3227 $scale = $scale-$zad;
3231 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
3232 $scale = $dbd+$scale;
3238 # round left from the '.'
3240 # 123 => 100 means length(123) = 3 - $scale (2) => 1
3242 my $dbt = $MBI->_len($x->{_m});
3244 my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
3245 # should be the same, so treat it as this
3246 $scale = 1 if $scale == 0;
3247 # shortcut if already integer
3248 return $x if $scale == 1 && $dbt <= $dbd;
3249 # maximum digits before dot
3254 # not enough digits before dot, so round to zero
3257 elsif ( $scale == $dbd )
3264 $scale = $dbd - $scale;
3267 # pass sign to bround for rounding modes '+inf' and '-inf'
3268 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
3269 $m->bround($scale,$mode);
3270 $x->{_m} = $m->{value}; # get our mantissa back
3276 # accuracy: preserve $N digits, and overwrite the rest with 0's
3277 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
3279 if (($_[0] || 0) < 0)
3281 require Carp; Carp::croak ('bround() needs positive accuracy');
3284 my ($scale,$mode) = $x->_scale_a(@_);
3285 return $x if !defined $scale || $x->modify('bround'); # no-op
3287 # scale is now either $x->{_a}, $accuracy, or the user parameter
3288 # test whether $x already has lower accuracy, do nothing in this case
3289 # but do round if the accuracy is the same, since a math operation might
3290 # want to round a number with A=5 to 5 digits afterwards again
3291 return $x if defined $x->{_a} && $x->{_a} < $scale;
3293 # scale < 0 makes no sense
3294 # scale == 0 => keep all digits
3295 # never round a +-inf, NaN
3296 return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
3298 # 1: never round a 0
3299 # 2: if we should keep more digits than the mantissa has, do nothing
3300 if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
3302 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
3306 # pass sign to bround for '+inf' and '-inf' rounding modes
3307 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
3309 $m->bround($scale,$mode); # round mantissa
3310 $x->{_m} = $m->{value}; # get our mantissa back
3311 $x->{_a} = $scale; # remember rounding
3312 delete $x->{_p}; # and clear P
3313 $x->bnorm(); # del trailing zeros gen. by bround()
3318 # return integer less or equal then $x
3319 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3321 return $x if $x->modify('bfloor');
3323 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3325 # if $x has digits after dot
3326 if ($x->{_es} eq '-')
3328 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3329 $x->{_e} = $MBI->_zero(); # trunc/norm
3330 $x->{_es} = '+'; # abs e
3331 $MBI->_inc($x->{_m}) if $x->{sign} eq '-'; # increment if negative
3333 $x->round($a,$p,$r);
3338 # return integer greater or equal then $x
3339 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3341 return $x if $x->modify('bceil');
3342 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3344 # if $x has digits after dot
3345 if ($x->{_es} eq '-')
3347 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3348 $x->{_e} = $MBI->_zero(); # trunc/norm
3349 $x->{_es} = '+'; # abs e
3350 $MBI->_inc($x->{_m}) if $x->{sign} eq '+'; # increment if positive
3352 $x->round($a,$p,$r);
3357 # shift right by $y (divide by power of $n)
3360 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3361 # objectify is costly, so avoid it
3362 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3364 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3367 return $x if $x->modify('brsft');
3368 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3370 $n = 2 if !defined $n; $n = $self->new($n);
3373 return $x->blsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3375 # the following call to bdiv() will return either quo or (quo,reminder):
3376 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
3381 # shift left by $y (multiply by power of $n)
3384 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3385 # objectify is costly, so avoid it
3386 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3388 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3391 return $x if $x->modify('blsft');
3392 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3394 $n = 2 if !defined $n; $n = $self->new($n);
3397 return $x->brsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3399 $x->bmul($n->bpow($y),$a,$p,$r,$y);
3402 ###############################################################################
3406 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
3411 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
3412 # or falling back to MBI::bxxx()
3413 my $name = $AUTOLOAD;
3415 $name =~ s/(.*):://; # split package
3416 my $c = $1 || $class;
3418 $c->import() if $IMPORT == 0;
3419 if (!_method_alias($name))
3423 # delayed load of Carp and avoid recursion
3425 Carp::croak ("$c: Can't call a method without name");
3427 if (!_method_hand_up($name))
3429 # delayed load of Carp and avoid recursion
3431 Carp::croak ("Can't call $c\-\>$name, not a valid method");
3433 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
3435 return &{"Math::BigInt"."::$name"}(@_);
3437 my $bname = $name; $bname =~ s/^f/b/;
3445 # return a copy of the exponent
3446 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3448 if ($x->{sign} !~ /^[+-]$/)
3450 my $s = $x->{sign}; $s =~ s/^[+-]//;
3451 return Math::BigInt->new($s); # -inf, +inf => +inf
3453 Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
3458 # return a copy of the mantissa
3459 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3461 if ($x->{sign} !~ /^[+-]$/)
3463 my $s = $x->{sign}; $s =~ s/^[+]//;
3464 return Math::BigInt->new($s); # -inf, +inf => +inf
3466 my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
3467 $m->bneg() if $x->{sign} eq '-';
3474 # return a copy of both the exponent and the mantissa
3475 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3477 if ($x->{sign} !~ /^[+-]$/)
3479 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
3480 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
3482 my $m = Math::BigInt->bzero();
3483 $m->{value} = $MBI->_copy($x->{_m});
3484 $m->bneg() if $x->{sign} eq '-';
3485 ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
3488 ##############################################################################
3489 # private stuff (internal use only)
3495 my $lib = ''; my @a;
3496 my $lib_kind = 'try';
3498 for ( my $i = 0; $i < $l ; $i++)
3500 if ( $_[$i] eq ':constant' )
3502 # This causes overlord er load to step in. 'binary' and 'integer'
3503 # are handled by BigInt.
3504 overload::constant float => sub { $self->new(shift); };
3506 elsif ($_[$i] eq 'upgrade')
3508 # this causes upgrading
3509 $upgrade = $_[$i+1]; # or undef to disable
3512 elsif ($_[$i] eq 'downgrade')
3514 # this causes downgrading
3515 $downgrade = $_[$i+1]; # or undef to disable
3518 elsif ($_[$i] =~ /^(lib|try|only)\z/)
3520 # alternative library
3521 $lib = $_[$i+1] || ''; # default Calc
3522 $lib_kind = $1; # lib, try or only
3525 elsif ($_[$i] eq 'with')
3527 # alternative class for our private parts()
3528 # XXX: no longer supported
3529 # $MBI = $_[$i+1] || 'Math::BigInt';
3538 $lib =~ tr/a-zA-Z0-9,://cd; # restrict to sane characters
3539 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
3540 my $mbilib = eval { Math::BigInt->config()->{lib} };
3541 if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
3543 # MBI already loaded
3544 Math::BigInt->import( $lib_kind, "$lib,$mbilib", 'objectify');
3548 # MBI not loaded, or with ne "Math::BigInt::Calc"
3549 $lib .= ",$mbilib" if defined $mbilib;
3550 $lib =~ s/^,//; # don't leave empty
3552 # replacement library can handle lib statement, but also could ignore it
3554 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
3555 # used in the same script, or eval inside import(). So we require MBI:
3556 require Math::BigInt;
3557 Math::BigInt->import( $lib_kind => $lib, 'objectify' );
3561 require Carp; Carp::croak ("Couldn't load $lib: $! $@");
3563 # find out which one was actually loaded
3564 $MBI = Math::BigInt->config()->{lib};
3566 # register us with MBI to get notified of future lib changes
3567 Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
3569 $self->export_to_level(1,$self,@a); # export wanted functions
3574 # adjust m and e so that m is smallest possible
3575 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
3577 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
3579 my $zeros = $MBI->_zeros($x->{_m}); # correct for trailing zeros
3582 my $z = $MBI->_new($zeros);
3583 $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
3584 if ($x->{_es} eq '-')
3586 if ($MBI->_acmp($x->{_e},$z) >= 0)
3588 $x->{_e} = $MBI->_sub ($x->{_e}, $z);
3589 $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
3593 $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
3599 $x->{_e} = $MBI->_add ($x->{_e}, $z);
3604 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
3605 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
3606 $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
3607 if $MBI->_is_zero($x->{_m});
3610 $x; # MBI bnorm is no-op, so dont call it
3613 ##############################################################################
3617 # return number as hexadecimal string (only for integers defined)
3618 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3620 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
3621 return '0x0' if $x->is_zero();
3623 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
3625 my $z = $MBI->_copy($x->{_m});
3626 if (! $MBI->_is_zero($x->{_e})) # > 0
3628 $MBI->_lsft( $z, $x->{_e},10);
3630 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3636 # return number as binary digit string (only for integers defined)
3637 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3639 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
3640 return '0b0' if $x->is_zero();
3642 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
3644 my $z = $MBI->_copy($x->{_m});
3645 if (! $MBI->_is_zero($x->{_e})) # > 0
3647 $MBI->_lsft( $z, $x->{_e},10);
3649 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3655 # return number as octal digit string (only for integers defined)
3656 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3658 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
3659 return '0' if $x->is_zero();
3661 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
3663 my $z = $MBI->_copy($x->{_m});
3664 if (! $MBI->_is_zero($x->{_e})) # > 0
3666 $MBI->_lsft( $z, $x->{_e},10);
3668 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3674 # return copy as a bigint representation of this BigFloat number
3675 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3677 return $x if $x->modify('as_number');
3679 if (!$x->isa('Math::BigFloat'))
3681 # if the object can as_number(), use it
3682 return $x->as_number() if $x->can('as_number');
3683 # otherwise, get us a float and then a number
3684 $x = $x->can('as_float') ? $x->as_float() : $self->new(0+"$x");
3687 return Math::BigInt->binf($x->sign()) if $x->is_inf();
3688 return Math::BigInt->bnan() if $x->is_nan();
3690 my $z = $MBI->_copy($x->{_m});
3691 if ($x->{_es} eq '-') # < 0
3693 $MBI->_rsft( $z, $x->{_e},10);
3695 elsif (! $MBI->_is_zero($x->{_e})) # > 0
3697 $MBI->_lsft( $z, $x->{_e},10);
3699 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3706 my $class = ref($x) || $x;
3707 $x = $class->new(shift) unless ref($x);
3709 return 1 if $MBI->_is_zero($x->{_m});
3711 my $len = $MBI->_len($x->{_m});
3712 $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
3716 $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
3727 Math::BigFloat - Arbitrary size floating point math package
3734 my $x = Math::BigFloat->new($str); # defaults to 0
3735 my $y = $x->copy(); # make a true copy
3736 my $nan = Math::BigFloat->bnan(); # create a NotANumber
3737 my $zero = Math::BigFloat->bzero(); # create a +0
3738 my $inf = Math::BigFloat->binf(); # create a +inf
3739 my $inf = Math::BigFloat->binf('-'); # create a -inf
3740 my $one = Math::BigFloat->bone(); # create a +1
3741 my $mone = Math::BigFloat->bone('-'); # create a -1
3743 my $pi = Math::BigFloat->bpi(100); # PI to 100 digits
3745 # the following examples compute their result to 100 digits accuracy:
3746 my $cos = Math::BigFloat->new(1)->bcos(100); # cosinus(1)
3747 my $sin = Math::BigFloat->new(1)->bsin(100); # sinus(1)
3748 my $atan = Math::BigFloat->new(1)->batan(100); # arcus tangens(1)
3750 my $atan2 = Math::BigFloat->new( 1 )->batan2( 1 ,100); # batan(1)
3751 my $atan2 = Math::BigFloat->new( 1 )->batan2( 8 ,100); # batan(1/8)
3752 my $atan2 = Math::BigFloat->new( -2 )->batan2( 1 ,100); # batan(-2)
3755 $x->is_zero(); # true if arg is +0
3756 $x->is_nan(); # true if arg is NaN
3757 $x->is_one(); # true if arg is +1
3758 $x->is_one('-'); # true if arg is -1
3759 $x->is_odd(); # true if odd, false for even
3760 $x->is_even(); # true if even, false for odd
3761 $x->is_pos(); # true if >= 0
3762 $x->is_neg(); # true if < 0
3763 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
3765 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
3766 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
3767 $x->sign(); # return the sign, either +,- or NaN
3768 $x->digit($n); # return the nth digit, counting from right
3769 $x->digit(-$n); # return the nth digit, counting from left
3771 # The following all modify their first argument. If you want to preserve
3772 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
3773 # necessary when mixing $a = $b assignments with non-overloaded math.
3776 $x->bzero(); # set $i to 0
3777 $x->bnan(); # set $i to NaN
3778 $x->bone(); # set $x to +1
3779 $x->bone('-'); # set $x to -1
3780 $x->binf(); # set $x to inf
3781 $x->binf('-'); # set $x to -inf
3783 $x->bneg(); # negation
3784 $x->babs(); # absolute value
3785 $x->bnorm(); # normalize (no-op)
3786 $x->bnot(); # two's complement (bit wise not)
3787 $x->binc(); # increment x by 1
3788 $x->bdec(); # decrement x by 1
3790 $x->badd($y); # addition (add $y to $x)
3791 $x->bsub($y); # subtraction (subtract $y from $x)
3792 $x->bmul($y); # multiplication (multiply $x by $y)
3793 $x->bdiv($y); # divide, set $x to quotient
3794 # return (quo,rem) or quo if scalar
3796 $x->bmod($y); # modulus ($x % $y)
3797 $x->bpow($y); # power of arguments ($x ** $y)
3798 $x->bmodpow($exp,$mod); # modular exponentation (($num**$exp) % $mod))
3799 $x->blsft($y, $n); # left shift by $y places in base $n
3800 $x->brsft($y, $n); # right shift by $y places in base $n
3801 # returns (quo,rem) or quo if in scalar context
3803 $x->blog(); # logarithm of $x to base e (Euler's number)
3804 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
3805 $x->bexp(); # calculate e ** $x where e is Euler's number
3807 $x->band($y); # bit-wise and
3808 $x->bior($y); # bit-wise inclusive or
3809 $x->bxor($y); # bit-wise exclusive or
3810 $x->bnot(); # bit-wise not (two's complement)
3812 $x->bsqrt(); # calculate square-root
3813 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
3814 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
3816 $x->bround($N); # accuracy: preserve $N digits
3817 $x->bfround($N); # precision: round to the $Nth digit
3819 $x->bfloor(); # return integer less or equal than $x
3820 $x->bceil(); # return integer greater or equal than $x
3822 # The following do not modify their arguments:
3824 bgcd(@values); # greatest common divisor
3825 blcm(@values); # lowest common multiplicator
3827 $x->bstr(); # return string
3828 $x->bsstr(); # return string in scientific notation
3830 $x->as_int(); # return $x as BigInt
3831 $x->exponent(); # return exponent as BigInt
3832 $x->mantissa(); # return mantissa as BigInt
3833 $x->parts(); # return (mantissa,exponent) as BigInt
3835 $x->length(); # number of digits (w/o sign and '.')
3836 ($l,$f) = $x->length(); # number of digits, and length of fraction
3838 $x->precision(); # return P of $x (or global, if P of $x undef)
3839 $x->precision($n); # set P of $x to $n
3840 $x->accuracy(); # return A of $x (or global, if A of $x undef)
3841 $x->accuracy($n); # set A $x to $n
3843 # these get/set the appropriate global value for all BigFloat objects
3844 Math::BigFloat->precision(); # Precision
3845 Math::BigFloat->accuracy(); # Accuracy
3846 Math::BigFloat->round_mode(); # rounding mode
3850 All operators (including basic math operations) are overloaded if you
3851 declare your big floating point numbers as
3853 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
3855 Operations with overloaded operators preserve the arguments, which is
3856 exactly what you expect.
3858 =head2 Canonical notation
3860 Input to these routines are either BigFloat objects, or strings of the
3861 following four forms:
3875 C</^[+-]\d+E[+-]?\d+$/>
3879 C</^[+-]\d*\.\d+E[+-]?\d+$/>
3883 all with optional leading and trailing zeros and/or spaces. Additionally,
3884 numbers are allowed to have an underscore between any two digits.
3886 Empty strings as well as other illegal numbers results in 'NaN'.
3888 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
3889 are always stored in normalized form. On a string, it creates a BigFloat
3894 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
3896 The string output will always have leading and trailing zeros stripped and drop
3897 a plus sign. C<bstr()> will give you always the form with a decimal point,
3898 while C<bsstr()> (s for scientific) gives you the scientific notation.
3900 Input bstr() bsstr()
3902 ' -123 123 123' '-123123123' '-123123123E0'
3903 '00.0123' '0.0123' '123E-4'
3904 '123.45E-2' '1.2345' '12345E-4'
3905 '10E+3' '10000' '1E4'
3907 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
3908 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
3909 return either undef, <0, 0 or >0 and are suited for sort.
3911 Actual math is done by using the class defined with C<< with => Class; >> (which
3912 defaults to BigInts) to represent the mantissa and exponent.
3914 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
3915 represent the result when input arguments are not numbers, as well as
3916 the result of dividing by zero.
3918 =head2 C<mantissa()>, C<exponent()> and C<parts()>
3920 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
3921 as BigInts such that:
3923 $m = $x->mantissa();
3924 $e = $x->exponent();
3925 $y = $m * ( 10 ** $e );
3926 print "ok\n" if $x == $y;
3928 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
3930 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
3932 Currently the mantissa is reduced as much as possible, favouring higher
3933 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
3934 This might change in the future, so do not depend on it.
3936 =head2 Accuracy vs. Precision
3938 See also: L<Rounding|Rounding>.
3940 Math::BigFloat supports both precision (rounding to a certain place before or
3941 after the dot) and accuracy (rounding to a certain number of digits). For a
3942 full documentation, examples and tips on these topics please see the large
3943 section about rounding in L<Math::BigInt>.
3945 Since things like C<sqrt(2)> or C<1 / 3> must presented with a limited
3946 accuracy lest a operation consumes all resources, each operation produces
3947 no more than the requested number of digits.
3949 If there is no global precision or accuracy set, B<and> the operation in
3950 question was not called with a requested precision or accuracy, B<and> the
3951 input $x has no accuracy or precision set, then a fallback parameter will
3952 be used. For historical reasons, it is called C<div_scale> and can be accessed
3955 $d = Math::BigFloat->div_scale(); # query
3956 Math::BigFloat->div_scale($n); # set to $n digits
3958 The default value for C<div_scale> is 40.
3960 In case the result of one operation has more digits than specified,
3961 it is rounded. The rounding mode taken is either the default mode, or the one
3962 supplied to the operation after the I<scale>:
3964 $x = Math::BigFloat->new(2);
3965 Math::BigFloat->accuracy(5); # 5 digits max
3966 $y = $x->copy()->bdiv(3); # will give 0.66667
3967 $y = $x->copy()->bdiv(3,6); # will give 0.666667
3968 $y = $x->copy()->bdiv(3,6,undef,'odd'); # will give 0.666667
3969 Math::BigFloat->round_mode('zero');
3970 $y = $x->copy()->bdiv(3,6); # will also give 0.666667
3972 Note that C<< Math::BigFloat->accuracy() >> and C<< Math::BigFloat->precision() >>
3973 set the global variables, and thus B<any> newly created number will be subject
3974 to the global rounding B<immediately>. This means that in the examples above, the
3975 C<3> as argument to C<bdiv()> will also get an accuracy of B<5>.
3977 It is less confusing to either calculate the result fully, and afterwards
3978 round it explicitly, or use the additional parameters to the math
3982 $x = Math::BigFloat->new(2);
3983 $y = $x->copy()->bdiv(3);
3984 print $y->bround(5),"\n"; # will give 0.66667
3989 $x = Math::BigFloat->new(2);
3990 $y = $x->copy()->bdiv(3,5); # will give 0.66667
3997 =item ffround ( +$scale )
3999 Rounds to the $scale'th place left from the '.', counting from the dot.
4000 The first digit is numbered 1.
4002 =item ffround ( -$scale )
4004 Rounds to the $scale'th place right from the '.', counting from the dot.
4008 Rounds to an integer.
4010 =item fround ( +$scale )
4012 Preserves accuracy to $scale digits from the left (aka significant digits)
4013 and pads the rest with zeros. If the number is between 1 and -1, the
4014 significant digits count from the first non-zero after the '.'
4016 =item fround ( -$scale ) and fround ( 0 )
4018 These are effectively no-ops.
4022 All rounding functions take as a second parameter a rounding mode from one of
4023 the following: 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'.
4025 The default rounding mode is 'even'. By using
4026 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
4027 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
4028 no longer supported.
4029 The second parameter to the round functions then overrides the default
4032 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
4033 'trunc' as rounding mode to make it equivalent to:
4038 You can override this by passing the desired rounding mode as parameter to
4041 $x = Math::BigFloat->new(2.5);
4042 $y = $x->as_number('odd'); # $y = 3
4046 Math::BigFloat supports all methods that Math::BigInt supports, except it
4047 calculates non-integer results when possible. Please see L<Math::BigInt>
4048 for a full description of each method. Below are just the most important
4053 $x->accuracy(5); # local for $x
4054 CLASS->accuracy(5); # global for all members of CLASS
4055 # Note: This also applies to new()!
4057 $A = $x->accuracy(); # read out accuracy that affects $x
4058 $A = CLASS->accuracy(); # read out global accuracy
4060 Set or get the global or local accuracy, aka how many significant digits the
4061 results have. If you set a global accuracy, then this also applies to new()!
4063 Warning! The accuracy I<sticks>, e.g. once you created a number under the
4064 influence of C<< CLASS->accuracy($A) >>, all results from math operations with
4065 that number will also be rounded.
4067 In most cases, you should probably round the results explicitly using one of
4068 L<round()>, L<bround()> or L<bfround()> or by passing the desired accuracy
4069 to the math operation as additional parameter:
4071 my $x = Math::BigInt->new(30000);
4072 my $y = Math::BigInt->new(7);
4073 print scalar $x->copy()->bdiv($y, 2); # print 4300
4074 print scalar $x->copy()->bdiv($y)->bround(2); # print 4300
4078 $x->precision(-2); # local for $x, round at the second digit right of the dot
4079 $x->precision(2); # ditto, round at the second digit left of the dot
4081 CLASS->precision(5); # Global for all members of CLASS
4082 # This also applies to new()!
4083 CLASS->precision(-5); # ditto
4085 $P = CLASS->precision(); # read out global precision
4086 $P = $x->precision(); # read out precision that affects $x
4088 Note: You probably want to use L<accuracy()> instead. With L<accuracy> you
4089 set the number of digits each result should have, with L<precision> you
4090 set the place where to round!
4094 $x->bexp($accuracy); # calculate e ** X
4096 Calculates the expression C<e ** $x> where C<e> is Euler's number.
4098 This method was added in v1.82 of Math::BigInt (April 2007).
4102 $x->bnok($y); # x over y (binomial coefficient n over k)
4104 Calculates the binomial coefficient n over k, also called the "choose"
4105 function. The result is equivalent to:
4111 This method was added in v1.84 of Math::BigInt (April 2007).
4115 print Math::BigFloat->bpi(100), "\n";
4117 Calculate PI to N digits (including the 3 before the dot). The result is
4118 rounded according to the current rounding mode, which defaults to "even".
4120 This method was added in v1.87 of Math::BigInt (June 2007).
4124 my $x = Math::BigFloat->new(1);
4125 print $x->bcos(100), "\n";
4127 Calculate the cosinus of $x, modifying $x in place.
4129 This method was added in v1.87 of Math::BigInt (June 2007).
4133 my $x = Math::BigFloat->new(1);
4134 print $x->bsin(100), "\n";
4136 Calculate the sinus of $x, modifying $x in place.
4138 This method was added in v1.87 of Math::BigInt (June 2007).
4142 my $y = Math::BigFloat->new(2);
4143 my $x = Math::BigFloat->new(3);
4144 print $y->batan2($x), "\n";
4146 Calculate the arcus tanges of C<$y> divided by C<$x>, modifying $y in place.
4147 See also L<batan()>.
4149 This method was added in v1.87 of Math::BigInt (June 2007).
4153 my $x = Math::BigFloat->new(1);
4154 print $x->batan(100), "\n";
4156 Calculate the arcus tanges of $x, modifying $x in place. See also L<batan2()>.
4158 This method was added in v1.87 of Math::BigInt (June 2007).
4164 Multiply $x by $y, and then add $z to the result.
4166 This method was added in v1.87 of Math::BigInt (June 2007).
4168 =head1 Autocreating constants
4170 After C<use Math::BigFloat ':constant'> all the floating point constants
4171 in the given scope are converted to C<Math::BigFloat>. This conversion
4172 happens at compile time.
4176 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
4178 prints the value of C<2E-100>. Note that without conversion of
4179 constants the expression 2E-100 will be calculated as normal floating point
4182 Please note that ':constant' does not affect integer constants, nor binary
4183 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
4188 Math with the numbers is done (by default) by a module called
4189 Math::BigInt::Calc. This is equivalent to saying:
4191 use Math::BigFloat lib => 'Calc';
4193 You can change this by using:
4195 use Math::BigFloat lib => 'GMP';
4197 B<Note>: General purpose packages should not be explicit about the library
4198 to use; let the script author decide which is best.
4200 Note: The keyword 'lib' will warn when the requested library could not be
4201 loaded. To suppress the warning use 'try' instead:
4203 use Math::BigFloat try => 'GMP';
4205 If your script works with huge numbers and Calc is too slow for them,
4206 you can also for the loading of one of these libraries and if none
4207 of them can be used, the code will die:
4209 use Math::BigFloat only => 'GMP,Pari';
4211 The following would first try to find Math::BigInt::Foo, then
4212 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
4214 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
4216 See the respective low-level library documentation for further details.
4218 Please note that Math::BigFloat does B<not> use the denoted library itself,
4219 but it merely passes the lib argument to Math::BigInt. So, instead of the need
4222 use Math::BigInt lib => 'GMP';
4225 you can roll it all into one line:
4227 use Math::BigFloat lib => 'GMP';
4229 It is also possible to just require Math::BigFloat:
4231 require Math::BigFloat;
4233 This will load the necessary things (like BigInt) when they are needed, and
4236 See L<Math::BigInt> for more details than you ever wanted to know about using
4237 a different low-level library.
4239 =head2 Using Math::BigInt::Lite
4241 For backwards compatibility reasons it is still possible to
4242 request a different storage class for use with Math::BigFloat:
4244 use Math::BigFloat with => 'Math::BigInt::Lite';
4246 However, this request is ignored, as the current code now uses the low-level
4247 math library for directly storing the number parts.
4251 C<Math::BigFloat> exports nothing by default, but can export the C<bpi()> method:
4253 use Math::BigFloat qw/bpi/;
4255 print bpi(10), "\n";
4259 Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
4263 Do not try to be clever to insert some operations in between switching
4266 require Math::BigFloat;
4267 my $matter = Math::BigFloat->bone() + 4; # load BigInt and Calc
4268 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
4269 my $anti_matter = Math::BigFloat->bone()+4; # now use Pari
4271 This will create objects with numbers stored in two different backend libraries,
4272 and B<VERY BAD THINGS> will happen when you use these together:
4274 my $flash_and_bang = $matter + $anti_matter; # Don't do this!
4278 =item stringify, bstr()
4280 Both stringify and bstr() now drop the leading '+'. The old code would return
4281 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
4282 reasoning and details.
4286 The following will probably not print what you expect:
4288 print $c->bdiv(123.456),"\n";
4290 It prints both quotient and reminder since print works in list context. Also,
4291 bdiv() will modify $c, so be careful. You probably want to use
4293 print $c / 123.456,"\n";
4294 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
4300 The following will probably not print what you expect:
4302 my $c = Math::BigFloat->new('3.14159');
4303 print $c->brsft(3,10),"\n"; # prints 0.00314153.1415
4305 It prints both quotient and remainder, since print calls C<brsft()> in list
4306 context. Also, C<< $c->brsft() >> will modify $c, so be careful.
4307 You probably want to use
4309 print scalar $c->copy()->brsft(3,10),"\n";
4310 # or if you really want to modify $c
4311 print scalar $c->brsft(3,10),"\n";
4315 =item Modifying and =
4319 $x = Math::BigFloat->new(5);
4322 It will not do what you think, e.g. making a copy of $x. Instead it just makes
4323 a second reference to the B<same> object and stores it in $y. Thus anything
4324 that modifies $x will modify $y (except overloaded math operators), and vice
4325 versa. See L<Math::BigInt> for details and how to avoid that.
4329 C<bpow()> now modifies the first argument, unlike the old code which left
4330 it alone and only returned the result. This is to be consistent with
4331 C<badd()> etc. The first will modify $x, the second one won't:
4333 print bpow($x,$i),"\n"; # modify $x
4334 print $x->bpow($i),"\n"; # ditto
4335 print $x ** $i,"\n"; # leave $x alone
4337 =item precision() vs. accuracy()
4339 A common pitfall is to use L<precision()> when you want to round a result to
4340 a certain number of digits:
4344 Math::BigFloat->precision(4); # does not do what you think it does
4345 my $x = Math::BigFloat->new(12345); # rounds $x to "12000"!
4346 print "$x\n"; # print "12000"
4347 my $y = Math::BigFloat->new(3); # rounds $y to "0"!
4348 print "$y\n"; # print "0"
4349 $z = $x / $y; # 12000 / 0 => NaN!
4351 print $z->precision(),"\n"; # 4
4353 Replacing L<precision> with L<accuracy> is probably not what you want, either:
4357 Math::BigFloat->accuracy(4); # enables global rounding:
4358 my $x = Math::BigFloat->new(123456); # rounded immediately to "12350"
4359 print "$x\n"; # print "123500"
4360 my $y = Math::BigFloat->new(3); # rounded to "3
4361 print "$y\n"; # print "3"
4362 print $z = $x->copy()->bdiv($y),"\n"; # 41170
4363 print $z->accuracy(),"\n"; # 4
4365 What you want to use instead is:
4369 my $x = Math::BigFloat->new(123456); # no rounding
4370 print "$x\n"; # print "123456"
4371 my $y = Math::BigFloat->new(3); # no rounding
4372 print "$y\n"; # print "3"
4373 print $z = $x->copy()->bdiv($y,4),"\n"; # 41150
4374 print $z->accuracy(),"\n"; # undef
4376 In addition to computing what you expected, the last example also does B<not>
4377 "taint" the result with an accuracy or precision setting, which would
4378 influence any further operation.
4384 L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
4385 L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
4387 The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
4388 because they solve the autoupgrading/downgrading issue, at least partly.
4390 The package at L<http://search.cpan.org/~tels/Math-BigInt> contains
4391 more documentation including a full version history, testcases, empty
4392 subclass files and benchmarks.
4396 This program is free software; you may redistribute it and/or modify it under
4397 the same terms as Perl itself.
4401 Mark Biggar, overloaded interface by Ilya Zakharevich.
4402 Completely rewritten by Tels L<http://bloodgate.com> in 2001 - 2006, and still