package Math::BigFloat;
#
-# Mike grinned. 'Two down, infinity to go' - Mike Nostrus in Before and After
+# Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
#
# The following hash values are internally used:
# _p: precision
# _f: flags, used to signal MBI not to touch our private parts
-$VERSION = '1.29';
+$VERSION = '1.32';
require 5.005;
use Exporter;
-use Math::BigInt qw/objectify/;
+use File::Spec;
+# use Math::BigInt;
@ISA = qw( Exporter Math::BigInt);
use strict;
$upgrade = undef;
$downgrade = undef;
+my $MBI = 'Math::BigInt'; # the package we are using for our private parts
+ # changable by use Math::BigFloat with => 'package'
##############################################################################
# the old code had $rnd_mode, so we need to support it, too
-$rnd_mode = 'even';
sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
sub FETCH { return $round_mode; }
sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
-BEGIN { tie $rnd_mode, 'Math::BigFloat'; }
+BEGIN
+ {
+ $rnd_mode = 'even';
+ tie $rnd_mode, 'Math::BigFloat';
+ }
##############################################################################
my %hand_ups = map { $_ => 1 }
qw / is_nan is_inf is_negative is_positive
accuracy precision div_scale round_mode fneg fabs babs fnot
- objectify
+ objectify upgrade downgrade
bone binf bnan bzero
/;
if ((ref($wanted)) && (ref($wanted) ne $class))
{
$self->{_m} = $wanted->as_number(); # get us a bigint copy
- $self->{_e} = Math::BigInt->bzero();
+ $self->{_e} = $MBI->bzero();
$self->{_m}->babs();
$self->{sign} = $wanted->sign();
return $self->bnorm();
# handle '+inf', '-inf' first
if ($wanted =~ /^[+-]?inf$/)
{
- $self->{_e} = Math::BigInt->bzero();
- $self->{_m} = Math::BigInt->bzero();
+ return $downgrade->new($wanted) if $downgrade;
+
+ $self->{_e} = $MBI->bzero();
+ $self->{_m} = $MBI->bzero();
$self->{sign} = $wanted;
$self->{sign} = '+inf' if $self->{sign} eq 'inf';
return $self->bnorm();
if (!ref $mis)
{
die "$wanted is not a number initialized to $class" if !$NaNOK;
- $self->{_e} = Math::BigInt->bzero();
- $self->{_m} = Math::BigInt->bzero();
+
+ return $downgrade->bnan() if $downgrade;
+
+ $self->{_e} = $MBI->bzero();
+ $self->{_m} = $MBI->bzero();
$self->{sign} = $nan;
}
else
{
# make integer from mantissa by adjusting exp, then convert to bigint
# undef,undef to signal MBI that we don't need no bloody rounding
- $self->{_e} = Math::BigInt->new("$$es$$ev",undef,undef); # exponent
- $self->{_m} = Math::BigInt->new("$$miv$$mfv",undef,undef); # create mant.
+ $self->{_e} = $MBI->new("$$es$$ev",undef,undef); # exponent
+ $self->{_m} = $MBI->new("$$miv$$mfv",undef,undef); # create mant.
# 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
$self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;
$self->{sign} = $$mis;
}
- # print "mbf new $self->{sign} $self->{_m} e $self->{_e}\n";
+ # if downgrade, inf, NaN or integers go down
+
+ if ($downgrade && $self->{_e}->{sign} eq '+')
+ {
+# print "downgrading $$miv$$mfv"."E$$es$$ev";
+ if ($self->{_e}->is_zero())
+ {
+ $self->{_m}->{sign} = $$mis; # negative if wanted
+ return $downgrade->new($self->{_m});
+ }
+ return $downgrade->new("$$mis$$miv$$mfv"."E$$es$$ev");
+ }
+ # print "mbf new $self->{sign} $self->{_m} e $self->{_e} ",ref($self),"\n";
$self->bnorm()->round(@r); # first normalize, then round
}
{
# used by parent class bone() to initialize number to 1
my $self = shift;
- $self->{_m} = Math::BigInt->bzero();
- $self->{_e} = Math::BigInt->bzero();
+ $self->{_m} = $MBI->bzero();
+ $self->{_e} = $MBI->bzero();
}
sub _binf
{
# used by parent class bone() to initialize number to 1
my $self = shift;
- $self->{_m} = Math::BigInt->bzero();
- $self->{_e} = Math::BigInt->bzero();
+ $self->{_m} = $MBI->bzero();
+ $self->{_e} = $MBI->bzero();
}
sub _bone
{
# used by parent class bone() to initialize number to 1
my $self = shift;
- $self->{_m} = Math::BigInt->bone();
- $self->{_e} = Math::BigInt->bzero();
+ $self->{_m} = $MBI->bone();
+ $self->{_e} = $MBI->bzero();
}
sub _bzero
{
# used by parent class bone() to initialize number to 1
my $self = shift;
- $self->{_m} = Math::BigInt->bzero();
- $self->{_e} = Math::BigInt->bone();
+ $self->{_m} = $MBI->bzero();
+ $self->{_e} = $MBI->bone();
+ }
+
+sub isa
+ {
+ my ($self,$class) = @_;
+ return if $class =~ /^Math::BigInt/; # we aren't one of these
+ UNIVERSAL::isa($self,$class);
+ }
+
+sub config
+ {
+ # return (later set?) configuration data as hash ref
+ my $class = shift || 'Math::BigFloat';
+
+ my $cfg = $MBI->config();
+
+ no strict 'refs';
+ $cfg->{class} = $class;
+ $cfg->{with} = $MBI;
+ foreach (
+ qw/upgrade downgrade precision accuracy round_mode VERSION div_scale/)
+ {
+ $cfg->{lc($_)} = ${"${class}::$_"};
+ };
+ $cfg;
}
##############################################################################
my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
- my $not_zero = !$x->is_zero();
+ my $not_zero = ! $x->is_zero();
if ($not_zero)
{
$es = $x->{_m}->bstr();
my $zeros = -$x->{_p} + $cad;
$es .= $dot.'0' x $zeros if $zeros > 0;
}
- return $es;
+ $es;
}
sub bsstr
}
my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
my $sep = 'e'.$sign;
- return $x->{_m}->bstr().$sep.$x->{_e}->bstr();
+ $x->{_m}->bstr().$sep.$x->{_e}->bstr();
}
sub numify
# Make a number from a BigFloat object
# simple return string and let Perl's atoi()/atof() handle the rest
my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
- return $x->bsstr();
+ $x->bsstr();
}
##############################################################################
# adjust so that exponents are equal
my $lxm = $x->{_m}->length();
my $lym = $y->{_m}->length();
- my $lx = $lxm + $x->{_e};
- my $ly = $lym + $y->{_e};
- my $l = $lx - $ly; $l->bneg() if $x->{sign} eq '-';
+ # the numify somewhat limits our length, but makes it much faster
+ my $lx = $lxm + $x->{_e}->numify();
+ my $ly = $lym + $y->{_e}->numify();
+ my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
return $l <=> 0 if $l != 0;
# lengths (corrected by exponent) are equal
- # so make mantissa euqal length by padding with zero (shift left)
+ # so make mantissa equal length by padding with zero (shift left)
my $diff = $lxm - $lym;
my $xm = $x->{_m}; # not yet copy it
my $ym = $y->{_m};
{
$xm = $x->{_m}->copy()->blsft(-$diff,10);
}
- my $rc = $xm->bcmp($ym);
+ my $rc = $xm->bacmp($ym);
$rc = -$rc if $x->{sign} eq '-'; # -124 < -123
$rc <=> 0;
}
# adjust so that exponents are equal
my $lxm = $x->{_m}->length();
my $lym = $y->{_m}->length();
- my $lx = $lxm + $x->{_e};
- my $ly = $lym + $y->{_e};
+ # the numify somewhat limits our length, but makes it much faster
+ my $lx = $lxm + $x->{_e}->numify();
+ my $ly = $lym + $y->{_e}->numify();
my $l = $lx - $ly;
return $l <=> 0 if $l != 0;
{
$xm = $x->{_m}->copy()->blsft(-$diff,10);
}
- $xm->bcmp($ym) <=> 0;
+ $xm->bacmp($ym) <=> 0;
}
sub badd
# return result as BFLOAT
my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
- #print "mbf badd $x $y\n";
# inf and NaN handling
if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
{
return $x if $x->{sign} eq $y->{sign};
return $x->bnan();
}
- # +-inf + something => +inf
- # something +-inf => +-inf
+ # +-inf + something => +inf; something +-inf => +-inf
$x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
return $x;
}
+ return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
+ ((!$x->isa($self)) || (!$y->isa($self)));
+
# speed: no add for 0+y or x+0
- return $x if $y->is_zero(); # x+0
+ return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
if ($x->is_zero()) # 0+y
{
# make copy, clobbering up x (modify in place!)
}
# take lower of the two e's and adapt m1 to it to match m2
- my $e = $y->{_e}; $e = Math::BigInt::bzero() if !defined $e; # if no BFLOAT
- $e = $e - $x->{_e};
+ my $e = $y->{_e};
+ $e = $MBI->bzero() if !defined $e; # if no BFLOAT ?
+ $e = $e->copy(); # make copy (didn't do it yet)
+ $e->bsub($x->{_e});
my $add = $y->{_m}->copy();
- if ($e < 0)
+ if ($e->{sign} eq '-') # < 0
{
my $e1 = $e->copy()->babs();
- $x->{_m} *= (10 ** $e1);
+ #$x->{_m} *= (10 ** $e1);
+ $x->{_m}->blsft($e1,10);
$x->{_e} += $e; # need the sign of e
}
- elsif ($e > 0)
+ elsif (!$e->is_zero()) # > 0
{
- $add *= (10 ** $e);
+ #$add *= (10 ** $e);
+ $add->blsft($e,10);
}
# else: both e are the same, so just leave them
$x->{_m}->{sign} = $x->{sign}; # fiddle with signs
# subtract second arg from first, modify first
my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
- if (!$y->is_zero()) # don't need to do anything if $y is 0
+ if ($y->is_zero()) # still round for not adding zero
{
- $y->{sign} =~ tr/+\-/-+/; # does nothing for NaN
- $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
- $y->{sign} =~ tr/+\-/-+/; # refix $y (does nothing for NaN)
+ return $x->round($a,$p,$r);
}
+
+ $y->{sign} =~ tr/+\-/-+/; # does nothing for NaN
+ $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
+ $y->{sign} =~ tr/+\-/-+/; # refix $y (does nothing for NaN)
$x; # already rounded by badd()
}
# http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
- # u = x-1, v = x +1
+ # u = x-1, v = x+1
# _ _
- # taylor: | u 1 u^3 1 u^5 |
+ # Taylor: | u 1 u^3 1 u^5 |
# ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
# |_ v 3 v^3 5 v^5 _|
+ # This takes much more steps to calculate the result:
+ # u = x-1
+ # _ _
+ # Taylor: | u 1 u^2 1 u^3 |
+ # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
+ # |_ x 2 x^2 3 x^3 _|
+
# we need to limit the accuracy to protect against overflow
my $fallback = 0;
my $scale = 0;
# we also need to disable any set A or P on $x (_find_round_parameters took
# them already into account), since these would interfere, too
delete $x->{_a}; delete $x->{_p};
- # need to disable $upgrade in BigInt, to aoid deep recursion
+ # need to disable $upgrade in BigInt, to avoid deep recursion
local $Math::BigInt::upgrade = undef;
-
- my $v = $x->copy(); $v->binc(); # v = x+1
- $x->bdec(); my $u = $x->copy(); # u = x-1; x = x-1
-
- $x->bdiv($v,$scale); # first term: u/v
-
- my $below = $v->copy();
- my $over = $u->copy();
- $u *= $u; $v *= $v; # u^2, v^2
- $below->bmul($v); # u^3, v^3
- $over->bmul($u);
- my $factor = $self->new(3); my $two = $self->new(2);
-
- my $diff = $self->bone();
- my $limit = $self->new("1E-". ($scale-1)); my $last;
- # print "diff $diff limit $limit\n";
- while ($diff->bcmp($limit) > 0)
- {
- #print "$x $over $below $factor\n";
- $diff = $x->copy()->bsub($last)->babs();
- #print "diff $diff $limit\n";
- $last = $x->copy();
- $x += $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
- $over *= $u; $below *= $v; $factor->badd($two);
- }
- $x->bmul($two);
+
+ my ($case,$limit,$v,$u,$below,$factor,$two,$next,$over,$f);
+
+ if (3 < 5)
+ #if ($x <= Math::BigFloat->new("0.5"))
+ {
+ $case = 0;
+ # print "case $case $x < 0.5\n";
+ $v = $x->copy(); $v->binc(); # v = x+1
+ $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
+ $x->bdiv($v,$scale); # first term: u/v
+ $below = $v->copy();
+ $over = $u->copy();
+ $u *= $u; $v *= $v; # u^2, v^2
+ $below->bmul($v); # u^3, v^3
+ $over->bmul($u);
+ $factor = $self->new(3); $f = $self->new(2);
+ }
+ #else
+ # {
+ # $case = 1;
+ # print "case 1 $x > 0.5\n";
+ # $v = $x->copy(); # v = x
+ # $u = $x->copy(); $u->bdec(); # u = x-1;
+ # $x->bdec(); $x->bdiv($v,$scale); # first term: x-1/x
+ # $below = $v->copy();
+ # $over = $u->copy();
+ # $below->bmul($v); # u^2, v^2
+ # $over->bmul($u);
+ # $factor = $self->new(2); $f = $self->bone();
+ # }
+ $limit = $self->new("1E-". ($scale-1));
+ #my $steps = 0;
+ while (3 < 5)
+ {
+ # we calculate the next term, and add it to the last
+ # when the next term is below our limit, it won't affect the outcome
+ # anymore, so we stop
+ $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
+ last if $next->bcmp($limit) <= 0;
+ $x->badd($next);
+ # print "step $x\n";
+ # calculate things for the next term
+ $over *= $u; $below *= $v; $factor->badd($f);
+ #$steps++;
+ }
+ $x->bmul(2) if $case == 0;
+ #print "took $steps steps\n";
# shortcut to not run trough _find_round_parameters again
if (defined $params[1])
}
# handle result = 0
return $x->bzero() if $x->is_zero() || $y->is_zero();
+
+ return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
+ ((!$x->isa($self)) || (!$y->isa($self)));
# aEb * cEd = (a*c)E(b+d)
$x->{_m}->bmul($y->{_m});
sub bdiv
{
# (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
- # (BFLOAT,BFLOAT) (quo,rem) or BINT (only rem)
+ # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
return $self->_div_inf($x,$y)
# x== 0 # also: or y == 1 or y == -1
return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
- # upgrade
- return $upgrade->bdiv($x,$y,$a,$p,$r) if defined $upgrade;
+ # upgrade ?
+ return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
# we need to limit the accuracy to protect against overflow
my $fallback = 0;
# promote BigInts and it's subclasses (except when already a BigFloat)
$y = $self->new($y) unless $y->isa('Math::BigFloat');
+ #print "bdiv $y ",ref($y),"\n";
+ # need to disable $upgrade in BigInt, to avoid deep recursion
+ local $Math::BigInt::upgrade = undef; # should be parent class vs MBI
+
# calculate the result to $scale digits and then round it
# a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
$x->{_m}->blsft($scale,10);
}
return ($x,$rem);
}
- return $x;
+ $x;
}
sub bmod
{
$x->{_m}->blsft($x->{_e},10);
}
- $x->{_e} = Math::BigInt->bzero() unless $x->{_e}->is_zero();
+ $x->{_e} = $MBI->bzero() unless $x->{_e}->is_zero();
$x->{_e}->bsub($shiftx) if $shiftx != 0;
$x->{_e}->bsub($shifty) if $shifty != 0;
# now mantissas are equalized, exponent of $x is adjusted, so calc result
-# $ym->{sign} = '-' if $neg; # bmod() will make the correction for us
$x->{_m}->bmod($ym);
}
# when user set globals, they would interfere with our calculation, so
- # disable then and later re-enable them
+ # disable them and later re-enable them
no strict 'refs';
my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
# we also need to disable any set A or P on $x (_find_round_parameters took
# them already into account), since these would interfere, too
delete $x->{_a}; delete $x->{_p};
- # need to disable $upgrade in BigInt, to aoid deep recursion
- local $Math::BigInt::upgrade = undef;
+ # need to disable $upgrade in BigInt, to avoid deep recursion
+ local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
my $xas = $x->as_number();
my $gs = $xas->copy()->bsqrt(); # some guess
+# print "guess $gs\n";
if (($x->{_e}->{sign} ne '-') # guess can't be accurate if there are
# digits after the dot
&& ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
{
# exact result
- $x->{_m} = $gs; $x->{_e} = Math::BigInt->bzero(); $x->bnorm();
+ $x->{_m} = $gs; $x->{_e} = $MBI->bzero(); $x->bnorm();
# shortcut to not run trough _find_round_parameters again
if (defined $params[1])
{
# clear a/p after round, since user did not request it
$x->{_a} = undef; $x->{_p} = undef;
}
+ # re-enable A and P, upgrade is taken care of by "local"
${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
return $x;
}
my $lx = $x->{_m}->length();
$scale = $lx if $scale < $lx;
my $e = $self->new("1E-$scale"); # make test variable
-# return $x->bnan() if $e->sign() eq 'NaN';
my $y = $x->copy();
my $two = $self->new(2);
$y = $self->new($y) unless $y->isa('Math::BigFloat');
my $rem;
- while ($diff >= $e)
+ while ($diff->bacmp($e) >= 0)
{
+ $rem = $y->copy()->bdiv($gs,$scale);
$rem = $y->copy()->bdiv($gs,$scale)->badd($gs)->bdiv($two,$scale);
- $diff = $rem->copy()->bsub($gs)->babs();
+ $diff = $rem->copy()->bsub($gs);
$gs = $rem->copy();
}
# copy over to modify $x
sub bfac
{
- # (BINT or num_str, BINT or num_str) return BINT
+ # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
# compute factorial numbers
# modifies first argument
my ($self,$x,@r) = objectify(1,@_);
- return $x->bnan() if $x->{sign} ne '+'; # inf, NnN, <0 etc => NaN
- return $x->bone(@r) if $x->is_zero() || $x->is_one(); # 0 or 1 => 1
+ return $x->bnan()
+ if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
+ ($x->{_e}->{sign} ne '+')); # digits after dot?
- return $x->bnan() if $x->{_e}->{sign} ne '+'; # digits after dot?
+ return $x->bone(@r) if $x->is_zero() || $x->is_one(); # 0 or 1 => 1
# use BigInt's bfac() for faster calc
$x->{_m}->blsft($x->{_e},10); # un-norm m
$x->{_e}->bzero(); # norm $x again
$x->{_m}->bfac(); # factorial
- $x->bnorm();
+ $x->bnorm()->round(@r);
+ }
+
+sub _pow2
+ {
+ # Calculate a power where $y is a non-integer, like 2 ** 0.5
+ my ($x,$y,$a,$p,$r) = @_;
+ my $self = ref($x);
+
+ # we need to limit the accuracy to protect against overflow
+ my $fallback = 0;
+ my $scale = 0;
+ my @params = $x->_find_round_parameters($a,$p,$r);
+
+ # no rounding at all, so must use fallback
+ if (scalar @params == 1)
+ {
+ # simulate old behaviour
+ $params[1] = $self->div_scale(); # and round to it as accuracy
+ $scale = $params[1]+4; # at least four more for proper round
+ $params[3] = $r; # round mode by caller or undef
+ $fallback = 1; # to clear a/p afterwards
+ }
+ else
+ {
+ # the 4 below is empirical, and there might be cases where it is not
+ # enough...
+ $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
+ }
+
+ # when user set globals, they would interfere with our calculation, so
+ # disable then and later re-enable them
+ no strict 'refs';
+ my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
+ my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
+ # we also need to disable any set A or P on $x (_find_round_parameters took
+ # them already into account), since these would interfere, too
+ delete $x->{_a}; delete $x->{_p};
+ # need to disable $upgrade in BigInt, to avoid deep recursion
+ local $Math::BigInt::upgrade = undef;
- #my $n = $x->copy();
- #$x->bone();
- #my $f = $self->new(2);
- #while ($f->bacmp($n) < 0)
- # {
- # $x->bmul($f); $f->binc();
- # }
- #$x->bmul($f); # last step
- $x->round(@r); # round
+ # split the second argument into its integer and fraction part
+ # we calculate the result then from these two parts, like in
+ # 2 ** 2.4 == (2 ** 2) * (2 ** 0.4)
+ my $c = $self->new($y->as_number()); # integer part
+ my $d = $y-$c; # fractional part
+ my $xc = $x->copy(); # a temp. copy
+
+ # now calculate binary fraction from the decimal fraction on the fly
+ # f.i. 0.654:
+ # 0.654 * 2 = 1.308 > 1 => 0.1 ( 1.308 - 1 = 0.308)
+ # 0.308 * 2 = 0.616 < 1 => 0.10
+ # 0.616 * 2 = 1.232 > 1 => 0.101 ( 1.232 - 1 = 0.232)
+ # and so on...
+ # The process stops when the result is exactly one, or when we have
+ # enough accuracy
+
+ # From the binary fraction we calculate the result as follows:
+ # we assume the fraction ends in 1, and we remove this one first.
+ # For each digit after the dot, assume 1 eq R and 0 eq XR, where R means
+ # take square root and X multiply with the original X.
+
+ my $i = 0;
+ while ($i++ < 50)
+ {
+ $d->badd($d); # * 2
+ last if $d->is_one(); # == 1
+ $x->bsqrt(); # 0
+ if ($d > 1)
+ {
+ $x->bsqrt(); $x->bmul($xc); $d->bdec(); # 1
+ }
+ }
+ # assume fraction ends in 1
+ $x->bsqrt(); # 1
+ if (!$c->is_one())
+ {
+ $x->bmul( $xc->bpow($c) );
+ }
+ elsif (!$c->is_zero())
+ {
+ $x->bmul( $xc );
+ }
+ # done
+
+ # shortcut to not run trough _find_round_parameters again
+ if (defined $params[1])
+ {
+ $x->bround($params[1],$params[3]); # then round accordingly
+ }
+ else
+ {
+ $x->bfround($params[2],$params[3]); # then round accordingly
+ }
+ if ($fallback)
+ {
+ # clear a/p after round, since user did not request it
+ $x->{_a} = undef; $x->{_p} = undef;
+ }
+ # restore globals
+ $$abr = $ab; $$pbr = $pb;
+ $x;
+ }
+
+sub _pow
+ {
+ # Calculate a power where $y is a non-integer, like 2 ** 0.5
+ my ($x,$y,$a,$p,$r) = @_;
+ my $self = ref($x);
+
+ # if $y == 0.5, it is sqrt($x)
+ return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
+
+ # u = y * ln x
+ # _ _
+ # Taylor: | u u^2 u^3 |
+ # x ** y = 1 + | --- + --- + * ----- + ... |
+ # |_ 1 1*2 1*2*3 _|
+
+ # we need to limit the accuracy to protect against overflow
+ my $fallback = 0;
+ my $scale = 0;
+ my @params = $x->_find_round_parameters($a,$p,$r);
+
+ # no rounding at all, so must use fallback
+ if (scalar @params == 1)
+ {
+ # simulate old behaviour
+ $params[1] = $self->div_scale(); # and round to it as accuracy
+ $scale = $params[1]+4; # at least four more for proper round
+ $params[3] = $r; # round mode by caller or undef
+ $fallback = 1; # to clear a/p afterwards
+ }
+ else
+ {
+ # the 4 below is empirical, and there might be cases where it is not
+ # enough...
+ $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
+ }
+
+ # when user set globals, they would interfere with our calculation, so
+ # disable then and later re-enable them
+ no strict 'refs';
+ my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
+ my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
+ # we also need to disable any set A or P on $x (_find_round_parameters took
+ # them already into account), since these would interfere, too
+ delete $x->{_a}; delete $x->{_p};
+ # need to disable $upgrade in BigInt, to avoid deep recursion
+ local $Math::BigInt::upgrade = undef;
+
+ my ($limit,$v,$u,$below,$factor,$next,$over);
+
+ $u = $x->copy()->blog($scale)->bmul($y);
+ $v = $self->bone(); # 1
+ $factor = $self->new(2); # 2
+ $x->bone(); # first term: 1
+
+ $below = $v->copy();
+ $over = $u->copy();
+
+ $limit = $self->new("1E-". ($scale-1));
+ #my $steps = 0;
+ while (3 < 5)
+ {
+ # we calculate the next term, and add it to the last
+ # when the next term is below our limit, it won't affect the outcome
+ # anymore, so we stop
+ $next = $over->copy()->bdiv($below,$scale);
+ last if $next->bcmp($limit) <= 0;
+ $x->badd($next);
+# print "at $x\n";
+ # calculate things for the next term
+ $over *= $u; $below *= $factor; $factor->binc();
+ #$steps++;
+ }
+
+ # shortcut to not run trough _find_round_parameters again
+ if (defined $params[1])
+ {
+ $x->bround($params[1],$params[3]); # then round accordingly
+ }
+ else
+ {
+ $x->bfround($params[2],$params[3]); # then round accordingly
+ }
+ if ($fallback)
+ {
+ # clear a/p after round, since user did not request it
+ $x->{_a} = undef; $x->{_p} = undef;
+ }
+ # restore globals
+ $$abr = $ab; $$pbr = $pb;
+ $x;
}
sub bpow
return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
return $x->bone() if $y->is_zero();
return $x if $x->is_one() || $y->is_one();
- my $y1 = $y->as_number(); # make bigint (trunc)
+
+ return $x->_pow($y,$a,$p,$r) if !$y->is_int(); # non-integer power
+
+ my $y1 = $y->as_number(); # make bigint
# if ($x == -1)
if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
{
# if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
return $y1->is_odd() ? $x : $x->babs(1);
}
- return $x if $x->is_zero() && $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
- # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
- return $x->binf() if $x->is_zero() && $y->{sign} eq '-';
+ if ($x->is_zero())
+ {
+ return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
+ # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
+ $x->binf();
+ }
# calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
$y1->babs();
my $z = $x->copy(); $x->bzero()->binc();
return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
}
- return $x->round($a,$p,$r,$y);
+ $x->round($a,$p,$r,$y);
}
###############################################################################
# if $x has digits after dot
if ($x->{_e}->{sign} eq '-')
{
- $x->{_m}->brsft(-$x->{_e},10);
- $x->{_e}->bzero();
- $x-- if $x->{sign} eq '-';
+ #$x->{_m}->brsft(-$x->{_e},10);
+ #$x->{_e}->bzero();
+ #$x-- if $x->{sign} eq '-';
+
+ $x->{_e}->{sign} = '+'; # negate e
+ $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
+ $x->{_e}->bzero(); # trunc/norm
+ $x->{_m}->binc() if $x->{sign} eq '-'; # decrement if negative
}
$x->round($a,$p,$r);
}
# if $x has digits after dot
if ($x->{_e}->{sign} eq '-')
{
- $x->{_m}->brsft(-$x->{_e},10);
- $x->{_e}->bzero();
- $x++ if $x->{sign} eq '+';
+ #$x->{_m}->brsft(-$x->{_e},10);
+ #$x->{_e}->bzero();
+ #$x++ if $x->{sign} eq '+';
+
+ $x->{_e}->{sign} = '+'; # negate e
+ $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
+ $x->{_e}->bzero(); # trunc/norm
+ $x->{_m}->binc() if $x->{sign} eq '+'; # decrement if negative
}
$x->round($a,$p,$r);
}
}
# try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
$name =~ s/^f/b/;
- return &{'Math::BigInt'."::$name"}(@_);
+ return &{"$MBI"."::$name"}(@_);
}
my $bname = $name; $bname =~ s/^f/b/;
*{$class."::$name"} = \&$bname;
sub import
{
my $self = shift;
- my $l = scalar @_; my $j = 0; my @a = @_;
- for ( my $i = 0; $i < $l ; $i++, $j++)
+ my $l = scalar @_;
+ my $lib = ''; my @a;
+ for ( my $i = 0; $i < $l ; $i++)
{
+# print "at $_[$i] (",$_[$i+1]||'undef',")\n";
if ( $_[$i] eq ':constant' )
{
# this rest causes overlord er load to step in
# print "overload @_\n";
overload::constant float => sub { $self->new(shift); };
- splice @a, $j, 1; $j--;
}
elsif ($_[$i] eq 'upgrade')
{
# this causes upgrading
- $upgrade = $_[$i+1]; # or undef to disable
- my $s = 2; $s = 1 if @a-$j < 2; # avoid "can not modify non-existant..."
- splice @a, $j, $s; $j -= $s;
+ $upgrade = $_[$i+1]; # or undef to disable
+ $i++;
+ }
+ elsif ($_[$i] eq 'downgrade')
+ {
+ # this causes downgrading
+ $downgrade = $_[$i+1]; # or undef to disable
+ $i++;
+ }
+ elsif ($_[$i] eq 'lib')
+ {
+ $lib = $_[$i+1] || ''; # default Calc
+ $i++;
+ }
+ elsif ($_[$i] eq 'with')
+ {
+ $MBI = $_[$i+1] || 'Math::BigInt'; # default Math::BigInt
+ $i++;
+ }
+ else
+ {
+ push @a, $_[$i];
}
}
+# print "mbf @a\n";
+
+ # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
+ my $mbilib = eval { Math::BigInt->config()->{lib} };
+ if ((defined $mbilib) && ($MBI eq 'Math::BigInt'))
+ {
+ # MBI already loaded
+ $MBI->import('lib',"$lib,$mbilib", 'objectify');
+ }
+ else
+ {
+ # MBI not loaded, or with ne "Math::BigInt"
+ $lib .= ",$mbilib" if defined $mbilib;
+
+# my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
+# my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
+# $file = File::Spec->catfile (@parts, $file);
+
+ if ($] < 5.006)
+ {
+ # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
+ # used in the same script, or eval inside import().
+ my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
+ my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
+ $file = File::Spec->catfile (@parts, $file);
+ eval { require $file; $MBI->import( lib => '$lib', 'objectify' ); }
+ }
+ else
+ {
+ my $rc = "use $MBI lib => '$lib', 'objectify';";
+ eval $rc;
+ }
+ }
+ die ("Couldn't load $MBI: $! $@") if $@;
+
# any non :constant stuff is handled by our parent, Exporter
# even if @_ is empty, to give it a chance
$self->SUPER::import(@a); # for subclasses
return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
- my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
- if ($zeros != 0)
- {
- $x->{_m}->brsft($zeros,10); $x->{_e} += $zeros;
- }
- # for something like 0Ey, set y to 1, and -0 => +0
- $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
+# if (!$x->{_m}->is_odd())
+# {
+ my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
+ if ($zeros != 0)
+ {
+ $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros);
+ }
+ # for something like 0Ey, set y to 1, and -0 => +0
+ $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
+# }
# this is to prevent automatically rounding when MBI's globals are set
$x->{_m}->{_f} = MB_NEVER_ROUND;
$x->{_e}->{_f} = MB_NEVER_ROUND;
# return copy as a bigint representation of this BigFloat number
my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
- my $z;
- if ($x->{_e}->is_zero())
- {
- $z = $x->{_m}->copy();
- $z->{sign} = $x->{sign};
- return $z;
- }
- $z = $x->{_m}->copy();
- if ($x->{_e} < 0)
+ my $z = $x->{_m}->copy();
+ if ($x->{_e}->{sign} eq '-') # < 0
{
- $z->brsft(-$x->{_e},10);
+ $x->{_e}->{sign} = '+'; # flip
+ $z->brsft($x->{_e},10);
+ $x->{_e}->{sign} = '-'; # flip back
}
- else
+ elsif (!$x->{_e}->is_zero()) # > 0
{
$z->blsft($x->{_e},10);
}
$len += $x->{_e} if $x->{_e}->sign() eq '+';
if (wantarray())
{
- my $t = Math::BigInt::bzero();
+ my $t = $MBI->bzero();
$t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
return ($len,$t);
}
perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
-prints the value of C<2E-100>. Note that without conversion of
+prints the value of C<2E-100>. Note that without conversion of
constants the expression 2E-100 will be calculated as normal floating point
number.
+Please note that ':constant' does not affect integer constants, nor binary
+nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
+work.
+
+=head2 Math library
+
+Math with the numbers is done (by default) by a module called
+Math::BigInt::Calc. This is equivalent to saying:
+
+ use Math::BigFloat lib => 'Calc';
+
+You can change this by using:
+
+ use Math::BigFloat lib => 'BitVect';
+
+The following would first try to find Math::BigInt::Foo, then
+Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
+
+ use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
+
+Calc.pm uses as internal format an array of elements of some decimal base
+(usually 1e7, but this might be differen for some systems) with the least
+significant digit first, while BitVect.pm uses a bit vector of base 2, most
+significant bit first. Other modules might use even different means of
+representing the numbers. See the respective module documentation for further
+details.
+
+Please note that Math::BigFloat does B<not> use the denoted library itself,
+but it merely passes the lib argument to Math::BigInt. So, instead of the need
+to do:
+
+ use Math::BigInt lib => 'GMP';
+ use Math::BigFloat;
+
+you can roll it all into one line:
+
+ use Math::BigFloat lib => 'GMP';
+
+Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details.
+
+=head2 Using Math::BigInt::Lite
+
+It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
+
+ # 1
+ use Math::BigFloat with => 'Math::BigInt::Lite';
+
+There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
+can combine these if you want. For instance, you may want to use
+Math::BigInt objects in your main script, too.
+
+ # 2
+ use Math::BigInt;
+ use Math::BigFloat with => 'Math::BigInt::Lite';
+
+Of course, you can combine this with the C<lib> parameter.
+
+ # 3
+ use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
+
+If you want to use Math::BigInt's, too, simple add a Math::BigInt B<before>:
+
+ # 4
+ use Math::BigInt;
+ use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
+
+Notice that the module with the last C<lib> will "win" and thus
+it's lib will be used if the lib is available:
+
+ # 5
+ use Math::BigInt lib => 'Bar,Baz';
+ use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
+
+That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
+words, Math::BigFloat will try to retain previously loaded libs when you
+don't specify it one.
+
+Actually, the lib loading order would be "Bar,Baz,Calc", and then
+"Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
+same as trying the latter load alone, except for the fact that Bar or Baz
+might be loaded needlessly in an intermidiate step
+
+The old way still works though:
+
+ # 6
+ use Math::BigInt lib => 'Bar,Baz';
+ use Math::BigFloat;
+
+But B<examples #3 and #4 are recommended> for usage.
+
=head1 BUGS
=over 2
projects / p5sagit/p5-mst-13.2.git / blobdiff