# _a : accuracy
# _p : precision
-$VERSION = '1.54';
+$VERSION = '1.57';
require 5.006002;
require Exporter;
{
($self,$x,$y,$a,$p,$r) = objectify(2,@_);
}
+
+ return $x if $x->modify('badd');
# inf and NaN handling
if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
# increment arg by one
my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
+ return $x if $x->modify('binc');
+
if ($x->{_es} eq '-')
{
return $x->badd($self->bone(),@r); # digits after dot
# decrement arg by one
my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
+ return $x if $x->modify('bdec');
+
if ($x->{_es} eq '-')
{
return $x->badd($self->bone('-'),@r); # digits after dot
{
my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
+ return $x if $x->modify('blog');
+
# $base > 0, $base != 1; if $base == undef default to $base == e
# $x >= 0
}
return $x->bzero(@params) if $x->is_one();
- # base not defined => base == Euler's constant e
+ # base not defined => base == Euler's number e
if (defined $base)
{
# make object, since we don't feed it through objectify() to still get the
$x;
}
+sub _len_to_steps
+ {
+ # Given D (digits in decimal), compute N so that N! (N factorial) is
+ # at least D digits long. D should be at least 50.
+ my $d = shift;
+
+ # two constants for the Ramanujan estimate of ln(N!)
+ my $lg2 = log(2 * 3.14159265) / 2;
+ my $lg10 = log(10);
+
+ # D = 50 => N => 42, so L = 40 and R = 50
+ my $l = 40; my $r = $d;
+
+ # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
+ $l = $l->numify if ref($l);
+ $r = $r->numify if ref($r);
+ $lg2 = $lg2->numify if ref($lg2);
+ $lg10 = $lg10->numify if ref($lg10);
+
+ # binary search for the right value (could this be written as the reverse of lg(n!)?)
+ while ($r - $l > 1)
+ {
+ my $n = int(($r - $l) / 2) + $l;
+ my $ramanujan =
+ int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
+ $ramanujan > $d ? $r = $n : $l = $n;
+ }
+ $l;
+ }
+
+sub bnok
+ {
+ # Calculate n over k (binomial coefficient or "choose" function) as integer.
+ # set up parameters
+ my ($self,$x,$y,@r) = (ref($_[0]),@_);
+
+ # objectify is costly, so avoid it
+ if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
+ {
+ ($self,$x,$y,@r) = objectify(2,@_);
+ }
+
+ return $x if $x->modify('bnok');
+
+ return $x->bnan() if $x->is_nan() || $y->is_nan();
+ return $x->binf() if $x->is_inf();
+
+ my $u = $x->as_int();
+ $u->bnok($y->as_int());
+
+ $x->{_m} = $u->{value};
+ $x->{_e} = $MBI->_zero();
+ $x->{_es} = '+';
+ $x->{sign} = '+';
+ $x->bnorm(@r);
+ }
+
sub bexp
{
- # Calculate e ** X (Euler's constant to the power of X)
+ # Calculate e ** X (Euler's number to the power of X)
my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
+ return $x if $x->modify('bexp');
+
return $x->binf() if $x->{sign} eq '+inf';
return $x->bzero() if $x->{sign} eq '-inf';
local $Math::BigFloat::downgrade = undef;
my $x_org = $x->copy();
- delete $x_org->{_a}; delete $x_org->{_p};
# We use the following Taylor series:
# -- + -- + -- + -- + -- + --- + --- + ---- = -----
# 1 1 2 6 24 120 720 5040 5040
- # Note that we cannot simple reduce 13700/5040 to 685/252.
+ # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!
- # So we start with A / B = 13490 / 5040 and F = 8
- my $A = $MBI->_new(13700);
- my $B = $MBI->_new(5040);
- my $F = $MBI->_new(8);
-
- # Code based on big rational math:
- my $limit = $MBI->_new('1' . '0' x ($scale - 2)); # scale=5 => 100000
-
- # while $B is not yet too big (making 1/$B too small)
- while ($MBI->_acmp($B,$limit) < 0)
+ if ($scale <= 75)
{
- # calculate $a * $f + 1
- $A = $MBI->_mul($A, $F);
- $A = $MBI->_inc($A);
- # calculate $b * $f
- $B = $MBI->_mul($B, $F);
- # increment f
- $F = $MBI->_inc($F);
+ # set $x directly from a cached string form
+ $x->{_m} = $MBI->_new(
+ "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
+ $x->{sign} = '+';
+ $x->{_es} = '-';
+ $x->{_e} = $MBI->_new(79);
}
+ else
+ {
+ # compute A and B so that e = A / B.
+
+ # After some terms we end up with this, so we use it as a starting point:
+ my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
+ my $F = $MBI->_new(42); my $step = 42;
+
+ # Compute how many steps we need to take to get $A and $B sufficiently big
+ my $steps = _len_to_steps($scale - 4);
+# print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
+ while ($step++ <= $steps)
+ {
+ # calculate $a * $f + 1
+ $A = $MBI->_mul($A, $F);
+ $A = $MBI->_inc($A);
+ # increment f
+ $F = $MBI->_inc($F);
+ }
+ # compute $B as factorial of $steps (this is faster than doing it manually)
+ my $B = $MBI->_fac($MBI->_new($steps));
+
+# print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
- $x = $self->new($MBI->_str($A))->bdiv($MBI->_str($B), $scale);
+ # compute A/B with $scale digits in the result (truncate, not round)
+ $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
+ $A = $MBI->_div( $A, $B );
- delete $x->{_a}; delete $x->{_p};
- # raise $x to the wanted power
- $x->bpow($x_org, $scale) unless $x_org->is_one();
+ $x->{_m} = $A;
+ $x->{sign} = '+';
+ $x->{_es} = '-';
+ $x->{_e} = $MBI->_new($scale);
+ }
- # shortcut to not run through _find_round_parameters again
- if (defined $params[0])
+ # $x contains now an estimate of e, with some surplus digits, so we can round
+ if (!$x_org->is_one())
{
- $x->bround($params[0],$params[2]); # then round accordingly
+ # raise $x to the wanted power and round it in one step:
+ $x->bpow($x_org, @params);
}
else
{
- $x->bfround($params[1],$params[2]); # then round accordingly
+ # else just round the already computed result
+ delete $x->{_a}; delete $x->{_p};
+ # shortcut to not run through _find_round_parameters again
+ if (defined $params[0])
+ {
+ $x->bround($params[0],$params[2]); # then round accordingly
+ }
+ else
+ {
+ $x->bfround($params[1],$params[2]); # then round accordingly
+ }
}
if ($fallback)
{
($self,$x,$y,$a,$p,$r) = objectify(2,@_);
}
+ return $x if $x->modify('bmul');
+
return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
# inf handling
($self,$x,$y,$a,$p,$r) = objectify(2,@_);
}
+ return $x if $x->modify('bdiv');
+
return $self->_div_inf($x,$y)
if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
($self,$x,$y,$a,$p,$r) = objectify(2,@_);
}
+ return $x if $x->modify('bmod');
+
# handle NaN, inf, -inf
if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
{
($self,$x,$y,$a,$p,$r) = objectify(2,@_);
}
+ return $x if $x->modify('broot');
+
# NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
$y->{sign} !~ /^\+$/;
# calculate square root
my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
+ return $x if $x->modify('bsqrt');
+
return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
# objectify is costly, so avoid it
($self,$x,@r) = objectify(1,@_) if !ref($x);
- return $x if $x->{sign} eq '+inf'; # inf => inf
+ # inf => inf
+ return $x if $x->modify('bfac') || $x->{sign} eq '+inf';
+
return $x->bnan()
if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
($x->{_es} ne '+')); # digits after dot?
($self,$x,$y,$a,$p,$r) = objectify(2,@_);
}
+ return $x if $x->modify('bpow');
+
return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
return $x if $x->{sign} =~ /^[+-]inf$/;
my $self = shift;
my $l = scalar @_;
my $lib = ''; my @a;
+ my $lib_kind = 'try';
$IMPORT=1;
for ( my $i = 0; $i < $l ; $i++)
{
$downgrade = $_[$i+1]; # or undef to disable
$i++;
}
- elsif ($_[$i] eq 'lib')
+ elsif ($_[$i] =~ /^(lib|try|only)\z/)
{
# alternative library
$lib = $_[$i+1] || ''; # default Calc
+ $lib_kind = $1; # lib, try or only
$i++;
}
elsif ($_[$i] eq 'with')
if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
{
# MBI already loaded
- Math::BigInt->import('try',"$lib,$mbilib", 'objectify');
+ Math::BigInt->import( $lib_kind, "$lib,$mbilib", 'objectify');
}
else
{
# Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
# used in the same script, or eval inside import(). So we require MBI:
require Math::BigInt;
- Math::BigInt->import( try => $lib, 'objectify' );
+ Math::BigInt->import( $lib_kind => $lib, 'objectify' );
}
if ($@)
{
# return copy as a bigint representation of this BigFloat number
my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
+ return $x if $x->modify('as_number');
+
my $z = $MBI->_copy($x->{_m});
if ($x->{_es} eq '-') # < 0
{
my $class = "Math::BigInt";
use 5.006002;
-$VERSION = '1.82';
+$VERSION = '1.86';
@ISA = qw(Exporter);
@EXPORT_OK = qw(objectify bgcd blcm);
$x->round(@r);
}
+sub bnok
+ {
+ # Calculate n over k (binomial coefficient or "choose" function) as integer.
+ # set up parameters
+ my ($self,$x,$y,@r) = (ref($_[0]),@_);
+
+ # objectify is costly, so avoid it
+ if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
+ {
+ ($self,$x,$y,@r) = objectify(2,@_);
+ }
+
+ return $x if $x->modify('bnok');
+ return $x->bnan() if $x->{sign} eq 'NaN' || $y->{sign} eq 'NaN';
+ return $x->binf() if $x->{sign} eq '+inf';
+
+ # k > n or k < 0 => 0
+ my $cmp = $x->bacmp($y);
+ return $x->bzero() if $cmp < 0 || $y->{sign} =~ /^-/;
+ # k == n => 1
+ return $x->bone(@r) if $cmp == 0;
+
+ if ($CALC->can('_nok'))
+ {
+ $x->{value} = $CALC->_nok($x->{value},$y->{value});
+ }
+ else
+ {
+ # ( 7 ) 7! 7*6*5 * 4*3*2*1 7 * 6 * 5
+ # ( - ) = --------- = --------------- = ---------
+ # ( 3 ) 3! (7-3)! 3*2*1 * 4*3*2*1 3 * 2 * 1
+
+ # compute n - k + 2 (so we start with 5 in the example above)
+ my $z = $x - $y;
+ if (!$z->is_one())
+ {
+ $z->binc();
+ my $r = $z->copy(); $z->binc();
+ my $d = $self->new(2);
+ while ($z->bacmp($x) <= 0) # f < x ?
+ {
+ $r->bmul($z); $r->bdiv($d);
+ $z->binc(); $d->binc();
+ }
+ $x->{value} = $r->{value}; $x->{sign} = '+';
+ }
+ else { $x->bone(); }
+ }
+ $x->round(@r);
+ }
+
sub bexp
{
# Calculate e ** $x (Euler's number to the power of X), truncated to
my $u;
{
# run through Math::BigFloat unless told otherwise
+ require Math::BigFloat unless defined $upgrade;
local $upgrade = 'Math::BigFloat' unless defined $upgrade;
# calculate result, truncate it to integer
$u = $upgrade->bexp($upgrade->new($x),@r);
$x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
$x->bfac(); # factorial of $x (1*2*3*4*..$x)
+ $x->bnok($y); # x over y (binomial coefficient n over k)
+
$x->blog(); # logarithm of $x to base e (Euler's number)
$x->blog($base); # logarithm of $x to base $base (f.i. 2)
$x->bexp(); # calculate e ** $x where e is Euler's number
See also L<blog()>.
+=head2 bnok()
+
+ $x->bnok($y); # x over y (binomial coefficient n over k)
+
+Calculates the binomial coefficient n over k, also called the "choose"
+function. The result is equivalent to:
+
+ ( n ) n!
+ | - | = -------
+ ( k ) k!(n-k)!
+
+This method was added in v1.84 of Math::BigInt (April 2007).
+
=head2 blsft()
$x->blsft($y); # left shift in base 2
use strict;
# use warnings; # dont use warnings for older Perls
-use vars qw/$VERSION/;
-
-$VERSION = '0.48';
+our $VERSION = '0.50';
# Package to store unsigned big integers in decimal and do math with them
# global constants, flags and accessory
# announce that we are compatible with MBI v1.70 and up
-sub api_version () { 1; }
+sub api_version () { 2; }
# constants for easier life
my ($BASE,$BASE_LEN,$MBASE,$RBASE,$MAX_VAL,$BASE_LEN_SMALL);
sub _base_len
{
- # set/get the BASE_LEN and assorted other, connected values
- # used only be the testsuite, set is used only by the BEGIN block below
+ # Set/get the BASE_LEN and assorted other, connected values.
+ # Used only by the testsuite, the set variant is used only by the BEGIN
+ # block below:
shift;
my $b = shift;
undef &_mul;
undef &_div;
- # $caught & 1 != 0 => cannot use MUL
- # $caught & 2 != 0 => cannot use DIV
- # The parens around ($caught & 1) were important, indeed, if we would use
- # & here.
+ # ($caught & 1) != 0 => cannot use MUL
+ # ($caught & 2) != 0 => cannot use DIV
if ($caught == 2) # 2
{
# must USE_MUL since we cannot use DIV
$AND_MASK = __PACKAGE__->_new( ( 2 ** $AND_BITS ));
$XOR_MASK = __PACKAGE__->_new( ( 2 ** $XOR_BITS ));
$OR_MASK = __PACKAGE__->_new( ( 2 ** $OR_BITS ));
+
+ # We can compute the approximate lenght no faster than the real length:
+ *_alen = \&_len;
}
###############################################################################
[ 10 ];
}
+sub _1ex
+ {
+ # create a 1Ex
+ my $rem = $_[1] % $BASE_LEN; # remainder
+ my $parts = $_[1] / $BASE_LEN; # parts
+
+ # 000000, 000000, 100
+ [ (0) x $parts, '1' . ('0' x $rem) ];
+ }
+
sub _copy
{
# make a true copy
$xi = shift @prod || 0; # || 0 makes v5.005_3 happy
}
push @$xv, @prod;
- __strip_zeros($xv);
+ # can't have leading zeros
+# __strip_zeros($xv);
$xv;
}
# multiply two numbers in internal representation
# modifies first arg, second need not be different from first
my ($c,$xv,$yv) = @_;
-
+
if (@$yv == 1)
{
# shortcut for two small numbers, also handles $x == 0
my $y = $yv->[0]; my $car = 0;
foreach my $i (@$xv)
{
+ # old, slower code (before use integer;)
$i = $i * $y + $car; $car = int($i / $MBASE); $i -= $car * $MBASE;
+ #$i = $i * $y + $car; $i -= ($car = $i / $MBASE) * $MBASE;
}
push @$xv, $car if $car != 0;
return $xv;
for $yi (@$yv)
{
$prod = $xi * $yi + ($prod[$cty] || 0) + $car;
- $prod[$cty++] =
- $prod - ($car = int($prod / $MBASE)) * $MBASE;
+ $prod[$cty++] = $prod - ($car = int($prod / $MBASE)) * $MBASE;
}
$prod[$cty] += $car if $car; # need really to check for 0?
$xi = shift @prod || 0; # || 0 makes v5.005_3 happy
}
push @$xv, @prod;
- __strip_zeros($xv);
+ # can't have leading zeros
+# __strip_zeros($xv);
$xv;
}
sub _len
{
- # compute number of digits
+ # compute number of digits in base 10
# int() because add/sub sometimes leaves strings (like '00005') instead of
# '5' in this place, thus causing length() to report wrong length
$cx;
}
+sub _nok
+ {
+ # n over k
+ # ref to array, return ref to array
+ my ($c,$n,$k) = @_;
+
+ # ( 7 ) 7! 7*6*5 * 4*3*2*1 7 * 6 * 5
+ # ( - ) = --------- = --------------- = ---------
+ # ( 3 ) 3! (7-3)! 3*2*1 * 4*3*2*1 3 * 2 * 1
+
+ # compute n - k + 2 (so we start with 5 in the example above)
+ my $x = _copy($c,$n);
+
+ _sub($c,$n,$k);
+ if (!_is_one($c,$n))
+ {
+ _inc($c,$n);
+ my $f = _copy($c,$n); _inc($c,$f); # n = 5, f = 6, d = 2
+ my $d = _two($c);
+ while (_acmp($c,$f,$x) <= 0) # f < n ?
+ {
+ # n = (n * f / d) == 5 * 6 / 2 => n == 3
+ $n = _mul($c,$n,$f); $n = _div($c,$n,$d);
+ # f = 7, d = 3
+ _inc($c,$f); _inc($c,$d);
+ }
+ }
+ else
+ {
+ # keep ref to $n and set it to 1
+ splice (@$n,1); $n->[0] = 1;
+ }
+ $n;
+ }
+
+my @factorials = (
+ 1,
+ 1,
+ 2,
+ 2*3,
+ 2*3*4,
+ 2*3*4*5,
+ 2*3*4*5*6,
+ 2*3*4*5*6*7,
+);
+
sub _fac
{
# factorial of $x
# ref to array, return ref to array
my ($c,$cx) = @_;
- if ((@$cx == 1) && ($cx->[0] <= 2))
+ if ((@$cx == 1) && ($cx->[0] <= 7))
{
- $cx->[0] ||= 1; # 0 => 1, 1 => 1, 2 => 2
+ $cx->[0] = $factorials[$cx->[0]]; # 0 => 1, 1 => 1, 2 => 2 etc.
return $cx;
}
+ if ((@$cx == 1) && # we do this only if $x >= 12 and $x <= 7000
+ ($cx->[0] >= 12 && $cx->[0] < 7000))
+ {
+
+ # Calculate (k-j) * (k-j+1) ... k .. (k+j-1) * (k + j)
+ # See http://blogten.blogspot.com/2007/01/calculating-n.html
+ # The above series can be expressed as factors:
+ # k * k - (j - i) * 2
+ # We cache k*k, and calculate (j * j) as the sum of the first j odd integers
+
+ # This will not work when N exceeds the storage of a Perl scalar, however,
+ # in this case the algorithm would be way to slow to terminate, anyway.
+
+ # As soon as the last element of $cx is 0, we split it up and remember
+ # how many zeors we got so far. The reason is that n! will accumulate
+ # zeros at the end rather fast.
+ my $zero_elements = 0;
+
+ # If n is even, set n = n -1
+ my $k = _num($c,$cx); my $even = 1;
+ if (($k & 1) == 0)
+ {
+ $even = $k; $k --;
+ }
+ # set k to the center point
+ $k = ($k + 1) / 2;
+# print "k $k even: $even\n";
+ # now calculate k * k
+ my $k2 = $k * $k;
+ my $odd = 1; my $sum = 1;
+ my $i = $k - 1;
+ # keep reference to x
+ my $new_x = _new($c, $k * $even);
+ @$cx = @$new_x;
+ if ($cx->[0] == 0)
+ {
+ $zero_elements ++; shift @$cx;
+ }
+# print STDERR "x = ", _str($c,$cx),"\n";
+ my $BASE2 = int(sqrt($BASE))-1;
+ my $j = 1;
+ while ($j <= $i)
+ {
+ my $m = ($k2 - $sum); $odd += 2; $sum += $odd; $j++;
+ while ($j <= $i && ($m < $BASE2) && (($k2 - $sum) < $BASE2))
+ {
+ $m *= ($k2 - $sum);
+ $odd += 2; $sum += $odd; $j++;
+# print STDERR "\n k2 $k2 m $m sum $sum odd $odd\n"; sleep(1);
+ }
+ if ($m < $BASE)
+ {
+ _mul($c,$cx,[$m]);
+ }
+ else
+ {
+ _mul($c,$cx,$c->_new($m));
+ }
+ if ($cx->[0] == 0)
+ {
+ $zero_elements ++; shift @$cx;
+ }
+# print STDERR "Calculate $k2 - $sum = $m (x = ", _str($c,$cx),")\n";
+ }
+ # multiply in the zeros again
+ unshift @$cx, (0) x $zero_elements;
+ return $cx;
+ }
+
# go forward until $base is exceeded
# limit is either $x steps (steps == 100 means a result always too high) or
# $base.
$n = _copy($c,$cx);
}
- $cx->[0] = $last; splice (@$cx,1); # keep ref to $x
+ # Set $cx to the last result below $BASE (but keep ref to $x)
+ $cx->[0] = $last; splice (@$cx,1);
+ # As soon as the last element of $cx is 0, we split it up and remember
+ # how many zeors we got so far. The reason is that n! will accumulate
+ # zeros at the end rather fast.
my $zero_elements = 0;
# do left-over steps fit into a scalar?
if (ref $n eq 'ARRAY')
{
# No, so use slower inc() & cmp()
+ # ($n is at least $BASE here)
+ my $base_2 = int(sqrt($BASE)) - 1;
+ #print STDERR "base_2: $base_2\n";
+ while ($step < $base_2)
+ {
+ if ($cx->[0] == 0)
+ {
+ $zero_elements ++; shift @$cx;
+ }
+ my $b = $step * ($step + 1); $step += 2;
+ _mul($c,$cx,[$b]);
+ }
$step = [$step];
- while (_acmp($step,$n) <= 0)
+ while (_acmp($c,$step,$n) <= 0)
{
- # as soon as the last element of $cx is 0, we split it up and remember
- # how many zeors we got so far. The reason is that n! will accumulate
- # zeros at the end rather fast.
if ($cx->[0] == 0)
{
$zero_elements ++; shift @$cx;
else
{
# Yes, so we can speed it up slightly
- while ($step <= $n)
+
+# print "# left over steps $n\n";
+
+ my $base_4 = int(sqrt(sqrt($BASE))) - 2;
+ #print STDERR "base_4: $base_4\n";
+ my $n4 = $n - 4;
+ while ($step < $n4 && $step < $base_4)
{
- # When the last element of $cx is 0, we split it up and remember
- # how many we got so far. The reason is that n! will accumulate
- # zeros at the end rather fast.
if ($cx->[0] == 0)
{
$zero_elements ++; shift @$cx;
}
+ my $b = $step * ($step + 1); $step += 2; $b *= $step * ($step + 1); $step += 2;
+ _mul($c,$cx,[$b]);
+ }
+ my $base_2 = int(sqrt($BASE)) - 1;
+ my $n2 = $n - 2;
+ #print STDERR "base_2: $base_2\n";
+ while ($step < $n2 && $step < $base_2)
+ {
+ if ($cx->[0] == 0)
+ {
+ $zero_elements ++; shift @$cx;
+ }
+ my $b = $step * ($step + 1); $step += 2;
+ _mul($c,$cx,[$b]);
+ }
+ # do what's left over
+ while ($step <= $n)
+ {
_mul($c,$cx,[$step]); $step++;
+ if ($cx->[0] == 0)
+ {
+ $zero_elements ++; shift @$cx;
+ }
}
}
# multiply in the zeros again
- while ($zero_elements-- > 0)
- {
- unshift @$cx, 0;
- }
+ unshift @$cx, (0) x $zero_elements;
$cx; # return result
}
if (scalar @$x == 1)
{
- # fit's into one Perl scalar, so result can be computed directly
+ # fits into one Perl scalar, so result can be computed directly
$x->[0] = int(sqrt($x->[0]));
return $x;
}
}
else
{
- # fit's into one Perl scalar, so result can be computed directly
+ # fits into one Perl scalar, so result can be computed directly
# cannot use int() here, because it rounds wrongly (try
# (81 ** 3) ** (1/3) to see what I mean)
#$x->[0] = int( $x->[0] ** (1 / $n->[0]) );
# convert a decimal number to hex (ref to array, return ref to string)
my ($c,$x) = @_;
- # fit's into one element (handle also 0x0 case)
+ # fits into one element (handle also 0x0 case)
return sprintf("0x%x",$x->[0]) if @$x == 1;
my $x1 = _copy($c,$x);
# convert a decimal number to bin (ref to array, return ref to string)
my ($c,$x) = @_;
- # fit's into one element (and Perl recent enough), handle also 0b0 case
+ # fits into one element (and Perl recent enough), handle also 0b0 case
# handle zero case for older Perls
if ($] <= 5.005 && @$x == 1 && $x->[0] == 0)
{
# convert a decimal number to octal (ref to array, return ref to string)
my ($c,$x) = @_;
- # fit's into one element (handle also 0 case)
+ # fits into one element (handle also 0 case)
return sprintf("0%o",$x->[0]) if @$x == 1;
my $x1 = _copy($c,$x);
sub _from_oct
{
- # convert a octal number to decimal (ref to string, return ref to array)
+ # convert a octal number to decimal (string, return ref to array)
my ($c,$os) = @_;
# for older Perls, play safe
sub _from_hex
{
- # convert a hex number to decimal (ref to string, return ref to array)
+ # convert a hex number to decimal (string, return ref to array)
my ($c,$hs) = @_;
my $m = _new($c, 0x10000000); # 28 bit at a time (<32 bit!)
sub _from_bin
{
- # convert a hex number to decimal (ref to string, return ref to array)
+ # convert a hex number to decimal (string, return ref to array)
my ($c,$bs) = @_;
# instead of converting X (8) bit at a time, it is faster to "convert" the
The following functions MUST be defined in order to support the use by
Math::BigInt v1.70 or later:
- api_version() return API version, minimum 1 for v1.70
+ api_version() return API version, 1 for v1.70, 2 for v1.83
_new(string) return ref to new object from ref to decimal string
_zero() return a new object with value 0
_one() return a new object with value 1
_check(obj) check whether internal representation is still intact
return 0 for ok, otherwise error message as string
- _from_hex(str) return ref to new object from ref to hexadecimal string
- _from_bin(str) return ref to new object from ref to binary string
- _from_oct(str) return ref to new object from ref to octal string
+ _from_hex(str) return new object from a hexadecimal string
+ _from_bin(str) return new object from a binary string
+ _from_oct(str) return new object from an octal string
_as_hex(str) return string containing the value as
unsigned hex string, with the '0x' prepended.
_and(obj1,obj2) AND (bit-wise) object 1 with object 2
_or(obj1,obj2) OR (bit-wise) object 1 with object 2
- _mod(obj,obj) Return remainder of div of the 1st by the 2nd object
+ _mod(obj1,obj2) Return remainder of div of the 1st by the 2nd object
_sqrt(obj) return the square root of object (truncated to int)
_root(obj) return the n'th (n >= 3) root of obj (truncated to int)
_fac(obj) return factorial of object 1 (1*2*3*4..)
- _pow(obj,obj) return object 1 to the power of object 2
+ _pow(obj1,obj2) return object 1 to the power of object 2
return undef for NaN
_zeros(obj) return number of trailing decimal zeros
_modinv return inverse modulus
undef : unknown whether result is exactly RESULT
_gcd(obj,obj) return Greatest Common Divisor of two objects
+The following functions are REQUIRED for an api_version of 2 or greater:
+
+ _1ex($x) create the number 1Ex where x >= 0
+ _alen(obj) returns approximate count of the decimal digits of the
+ object. This estimate MUST always be greater or equal
+ to what _len() returns.
+ _nok(n,k) calculate n over k (binomial coefficient)
+
The following functions are optional, and can be defined if the underlying lib
has a fast way to do them. If undefined, Math::BigInt will use pure Perl (hence
slow) fallback routines to emulate these:
_signed_and
_signed_xor
-
Input strings come in as unsigned but with prefix (i.e. as '123', '0xabc'
or '0b1101').
The first parameter can be modified, that includes the possibility that you
return a reference to a completely different object instead. Although keeping
-the reference and just changing it's contents is preferred over creating and
+the reference and just changing its contents is preferred over creating and
returning a different reference.
Return values are always references to objects, strings, or true/false for
=head1 SEE ALSO
-L<Math::BigInt>, L<Math::BigFloat>, L<Math::BigInt::BitVect>,
+L<Math::BigInt>, L<Math::BigFloat>,
L<Math::BigInt::GMP>, L<Math::BigInt::FastCalc> and L<Math::BigInt::Pari>.
=cut
}
print "# INC = @INC\n";
- plan tests => 2042;
+ plan tests => 2064;
}
use Math::BigFloat lib => 'BareCalc';
}
print "# INC = @INC\n";
- plan tests => 3055;
+ plan tests => 3093;
}
use Math::BigInt lib => 'BareCalc';
$class = "Math::BigInt";
$CL = "Math::BigInt::BareCalc";
-my $version = '1.61'; # for $VERSION tests, match current release (by hand!)
+my $version = '1.84'; # for $VERSION tests, match current release (by hand!)
require 'bigintpm.inc'; # perform same tests as bigintpm
$try .= '$x->facmp($y);';
} elsif ($f eq "fpow") {
$try .= '$x ** $y;';
+ } elsif ($f eq "bnok") {
+ $try .= '$x->bnok($y);';
} elsif ($f eq "froot") {
$try .= "$setup; \$x->froot(\$y);";
} elsif ($f eq "fadd") {
+27:+90:270
+1034:+804:415668
$div_scale = 40;
+&bnok
++inf:10:inf
+NaN:NaN:NaN
+NaN:1:NaN
+1:NaN:NaN
+1:1:1
+# k > n
+1:2:0
+2:3:0
+# k < 0
+1:-2:0
+# 7 over 3 = 35
+7:3:35
+7:6:1
+100:90:17310309456440
&flog
0::NaN
-1::NaN
}
print "# INC = @INC\n";
- plan tests => 2042
+ plan tests => 2064
+ 3; # own tests
}
print "1..0\n";
exit(0);
}
- plan tests => 322;
+ plan tests => 375;
}
use Math::BigInt::Calc;
ok ($C->_is_even($C->_one()) || 0,0); ok ($C->_is_even($C->_zero()),1);
# _len
-$x = $C->_new("1"); ok ($C->_len($x),1);
-$x = $C->_new("12"); ok ($C->_len($x),2);
-$x = $C->_new("123"); ok ($C->_len($x),3);
-$x = $C->_new("1234"); ok ($C->_len($x),4);
-$x = $C->_new("12345"); ok ($C->_len($x),5);
-$x = $C->_new("123456"); ok ($C->_len($x),6);
-$x = $C->_new("1234567"); ok ($C->_len($x),7);
-$x = $C->_new("12345678"); ok ($C->_len($x),8);
-$x = $C->_new("123456789"); ok ($C->_len($x),9);
-
-$x = $C->_new("8"); ok ($C->_len($x),1);
-$x = $C->_new("21"); ok ($C->_len($x),2);
-$x = $C->_new("321"); ok ($C->_len($x),3);
-$x = $C->_new("4321"); ok ($C->_len($x),4);
-$x = $C->_new("54321"); ok ($C->_len($x),5);
-$x = $C->_new("654321"); ok ($C->_len($x),6);
-$x = $C->_new("7654321"); ok ($C->_len($x),7);
-$x = $C->_new("87654321"); ok ($C->_len($x),8);
-$x = $C->_new("987654321"); ok ($C->_len($x),9);
-
-for (my $i = 1; $i < 9; $i++)
+for my $method (qw/_alen _len/)
{
- my $a = "$i" . '0' x ($i-1);
- $x = $C->_new($a);
- print "# Tried len '$a'\n" unless ok ($C->_len($x),$i);
+ $x = $C->_new("1"); ok ($C->$method($x),1);
+ $x = $C->_new("12"); ok ($C->$method($x),2);
+ $x = $C->_new("123"); ok ($C->$method($x),3);
+ $x = $C->_new("1234"); ok ($C->$method($x),4);
+ $x = $C->_new("12345"); ok ($C->$method($x),5);
+ $x = $C->_new("123456"); ok ($C->$method($x),6);
+ $x = $C->_new("1234567"); ok ($C->$method($x),7);
+ $x = $C->_new("12345678"); ok ($C->$method($x),8);
+ $x = $C->_new("123456789"); ok ($C->$method($x),9);
+
+ $x = $C->_new("8"); ok ($C->$method($x),1);
+ $x = $C->_new("21"); ok ($C->$method($x),2);
+ $x = $C->_new("321"); ok ($C->$method($x),3);
+ $x = $C->_new("4321"); ok ($C->$method($x),4);
+ $x = $C->_new("54321"); ok ($C->$method($x),5);
+ $x = $C->_new("654321"); ok ($C->$method($x),6);
+ $x = $C->_new("7654321"); ok ($C->$method($x),7);
+ $x = $C->_new("87654321"); ok ($C->$method($x),8);
+ $x = $C->_new("987654321"); ok ($C->$method($x),9);
+
+ $x = $C->_new("0"); ok ($C->$method($x),1);
+ $x = $C->_new("20"); ok ($C->$method($x),2);
+ $x = $C->_new("320"); ok ($C->$method($x),3);
+ $x = $C->_new("4320"); ok ($C->$method($x),4);
+ $x = $C->_new("54320"); ok ($C->$method($x),5);
+ $x = $C->_new("654320"); ok ($C->$method($x),6);
+ $x = $C->_new("7654320"); ok ($C->$method($x),7);
+ $x = $C->_new("87654320"); ok ($C->$method($x),8);
+ $x = $C->_new("987654320"); ok ($C->$method($x),9);
+
+ for (my $i = 1; $i < 9; $i++)
+ {
+ my $a = "$i" . '0' x ($i-1);
+ $x = $C->_new($a);
+ print "# Tried len '$a'\n" unless ok ($C->_len($x),$i);
+ }
}
# _digit
print "# _pow( ", '9' x $i, ", 2) \n" unless
ok ($C->_str($C->_pow($x,$n)),$rc);
- if ($i <= 7)
+ # if $i > $BASE_LEN, the test takes a really long time:
+ if ($i <= $BASE_LEN)
{
$x = '9' x $i; $x = $C->_new($x);
$n = '9' x $i; $n = $C->_new($n);
+ print "# _root( ", '9' x $i, ", ", 9 x $i, ") \n";
print "# _root( ", '9' x $i, ", ", 9 x $i, ") \n" unless
ok ($C->_str($C->_root($x,$n)),'1');
print "# _root( ", '9' x $i, ", ", 9 x $i, ") \n" unless
ok ($C->_str($C->_root($x,$n)), $res->[$i-2]);
}
+ else
+ {
+ ok ("skipped $i", "skipped $i");
+ ok ("skipped $i", "skipped $i");
+ }
}
##############################################################################
ok ($C->_as_bin( $C->_new("12")), '0b1100');
ok ($C->_as_oct( $C->_new("64")), '0100');
+# _1ex
+ok ($C->_str($C->_1ex(0)), "1");
+ok ($C->_str($C->_1ex(1)), "10");
+ok ($C->_str($C->_1ex(2)), "100");
+ok ($C->_str($C->_1ex(12)), "1000000000000");
+ok ($C->_str($C->_1ex(16)), "10000000000000000");
+
# _check
$x = $C->_new("123456789");
ok ($C->_check($x),0);
$try .= '$m = $m->bstr(); $m = "NaN" if !defined $m;';
$try .= '$e = $e->bstr(); $e = "NaN" if !defined $e;';
$try .= '"$m,$e";';
+ }elsif ($f eq "bexp"){
+ $try .= "\$x->bexp();";
} else {
# binary ops
$try .= "\$y = $class->new('$args[1]');";
{
$try .= "\$x >> \$y;";
}
+ }elsif ($f eq "bnok"){
+ $try .= "\$x->bnok(\$y);";
}elsif ($f eq "broot"){
$try .= "\$x->broot(\$y);";
}elsif ($f eq "blog"){
10:3628800
11:39916800
12:479001600
+20:2432902008176640000
+22:1124000727777607680000
+69:171122452428141311372468338881272839092270544893520369393648040923257279754140647424000000000000000
&bpow
abc:12:NaN
12:abc:NaN
18446744073709551616:128:1
# 213 ** 15
84274086103068221283760416414557757:15:213
-# see t/bigroot for more tests
+# see t/bigroot.t for more tests
&bsqrt
145:12
144:12
Nan:NaN
+inf:inf
-inf:NaN
+# see t/biglog.t for more tests
+&bexp
+NaN:NaN
+inf:inf
+1:2
+2:7
+# see t/bignok.t for more tests
+&bnok
++inf:10:inf
+NaN:NaN:NaN
+NaN:1:NaN
+1:NaN:NaN
+1:1:1
+# k > n
+1:2:0
+2:3:0
+# k < 0
+1:-2:0
+# 7 over 3 = 35
+7:3:35
+7:6:1
+100:90:17310309456440
+100:95:75287520
&bround
$round_mode('trunc')
0:12:0
my $location = $0; $location =~ s/bigintpm.t//;
unshift @INC, $location; # to locate the testing files
chdir 't' if -d 't';
- plan tests => 3055;
+ plan tests => 3093;
}
use Math::BigInt lib => 'Calc';
}
print "# INC = @INC\n";
- plan tests => 66;
+ plan tests => 68;
}
use Math::BigFloat;
test_bpow ('0.2','0.41',10, '0.5169187652');
#############################################################################
-# test bexp()
+# test bexp() with cached results
is ($cl->new(1)->bexp(), '2.718281828459045235360287471352662497757', 'bexp(1)');
is ($cl->new(2)->bexp(40), $cl->new(1)->bexp(45)->bpow(2,40), 'bexp(2)');
is ($cl->new("12.5")->bexp(61), $cl->new(1)->bexp(65)->bpow(12.5,61), 'bexp(12.5)');
+#############################################################################
+# test bexp() with big values (non-cached)
+
+is ($cl->new(1)->bexp(100),
+ '2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427',
+ 'bexp(100)');
+
+is ($cl->new("12.5")->bexp(91), $cl->new(1)->bexp(95)->bpow(12.5,91),
+ 'bexp(12.5) to 91 digits');
+
# all done
1;
# But it is better to test the numerical functionality, instead of not testing
# it at all.
-use Test;
+use Test::More;
use strict;
BEGIN
my ($x,$n,$y,$scale,$result) = @_;
my $s = $scale || 'undef';
- print "# Try: $cl $x->bpow($n)->broot($y,$s) == $result:\n";
- ok ($cl->new($x)->bpow($n)->broot($y,$scale),$result);
+ is ($cl->new($x)->bpow($n)->broot($y,$scale),$result, "Try: $cl $x->bpow($n)->broot($y,$s) == $result");
$result =~ s/\..*//;
- print "# Try: $c $x->bpow($n)->broot($y,$s) == $result:\n";
- ok ($c->new($x)->bpow($n)->broot($y,$scale),$result);
+ is ($c->new($x)->bpow($n)->broot($y,$scale),$result, "Try: $c $x->bpow($n)->broot($y,$s) == $result");
}
}
print "# INC = @INC\n";
- plan tests => 2;
+ plan tests => 4;
}
# first load BigInt with Calc
is (Math::BigFloat::config()->{lib}, 'Math::BigInt::BareCalc', 'BigFloat was notified');
+# See that Math::BigFloat supports "only"
+eval "Math::BigFloat->import('only' => 'Calc')";
+is (Math::BigFloat::config()->{lib}, 'Math::BigInt::Calc', '"only" worked');
+
+# See that Math::BigFloat supports "try"
+eval "Math::BigFloat->import('try' => 'BareCalc')";
+is (Math::BigFloat::config()->{lib}, 'Math::BigInt::BareCalc', '"try" worked');
+
#!/usr/bin/perl -w
-use Test;
+use Test::More;
use strict;
my $count;
print "# A $A\n# B $B\n";
if ($A->is_zero() || $B->is_zero())
{
- for (1..4) { ok (1,1); } next;
+ for (1..4) { is (1,1, 'skipped this test'); } next;
}
# check that int(A/B)*B + A % B == A holds for all inputs
print "# seed $seed, ". join(' ',Math::BigInt::Calc->_base_len()),"\n".
"# tried $ADB * $B + $two*$AMB - $AMB\n"
- unless ok ($ADB*$B+$two*$AMB-$AMB,$As);
- if (ok ($ADB*$B/$B,$ADB))
+ unless is ($ADB*$B+$two*$AMB-$AMB,$As, "ADB * B + 2 * AMB - AMB == A");
+ if (is ($ADB*$B/$B,$ADB, "ADB * B / B == ADB"))
{
print "# seed: $seed, \$ADB * \$B / \$B = ", $ADB * $B / $B, " != $ADB (\$B=$B)\n";
if (Math::BigInt->config()->{lib} =~ /::Calc/)
# print "check: $ADB $AMB";
print "# seed $seed, ". join(' ',Math::BigInt::Calc->_base_len()),"\n".
"# tried $ADB * $A + $two*$AMB - $AMB\n"
- unless ok ($ADB*$A+$two*$AMB-$AMB,$Bs);
+ unless is ($ADB*$A+$two*$AMB-$AMB,$Bs, "ADB * A + 2 * AMB - AMB == B");
print "# +$two * $AMB = ",$ADB * $A + $two * $AMB,"\n";
print "# -$AMB = ",$ADB * $A + $two * $AMB - $AMB,"\n";
print "# seed $seed, \$ADB * \$A / \$A = ", $ADB * $A / $A, " != $ADB (\$A=$A)\n"
- unless ok ($ADB*$A/$A,$ADB);
+ unless is ($ADB*$A/$A,$ADB, "ADB * A/A == ADB");
}
}
print "# INC = @INC\n";
- plan tests => 2042
+ plan tests => 2064
+ 6; # + our own tests
}
}
print "# INC = @INC\n";
- plan tests => 3055
+ plan tests => 3093
+ 5; # +5 own tests
}
}
print "# INC = @INC\n";
- plan tests => 2042
+ plan tests => 2064
+ 1;
}
projects / p5sagit/p5-mst-13.2.git / commitdiff