This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.14
in repository libmath-prime-util-perl.
commit 0a75888e02a5d59db80296fb3d2e27ec437331ce
Author: Dana Jacobsen <d...@acm.org>
Date: Mon Nov 26 01:01:17 2012 -0800
Fixup 5.6.2, and some li and Ei range cases
---
lib/Math/Prime/Util.pm | 14 +++++++++++---
lib/Math/Prime/Util/PP.pm | 49 ++++++++++++++++++++++++++++++++---------------
t/16-randomprime.t | 8 ++++++++
t/18-functions.t | 12 ++++++++----
4 files changed, 61 insertions(+), 22 deletions(-)
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 5eab9a2..e40c547 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -541,7 +541,10 @@ sub primes {
$_random_ndigit_ranges[$digits] = [next_prime($low),
prev_prime($high)];
} else {
my $low = int(10 ** ($digits-1));
- my $high = int(10 ** $digits);
+ my $high = int(10 ** int($digits));
+ # Note: Perl 5.6.2 cannot represent 10**15 as an integer, so things
will
+ # crash all over the place if you try. We can stringify it, but then
it
+ # starts failing math tests later.
$high = ~0 if $high > ~0;
$_random_ndigit_ranges[$digits] = [next_prime($low),
prev_prime($high)];
}
@@ -1153,7 +1156,8 @@ sub prime_count_approx {
# my $result = int(LogarithmicIntegral($x) -
LogarithmicIntegral(sqrt($x))/2);
- my $tol = 10**-(length(int($x))+1);
+ my $xlen = (ref($x) eq 'Math::BigFloat') ? length($x->bfloor->bstr()) :
length(int($x));
+ my $tol = 10**-$xlen;
my $result = RiemannR($x, $tol) + 0.5;
return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq
'Math::BigFloat';
@@ -1414,7 +1418,9 @@ sub RiemannR {
sub ExponentialIntegral {
my($n) = @_;
- croak "Invalid input to ExponentialIntegral: x must be != 0" if $n == 0;
+ return 0+'-inf' if $n == 0;
+ return 0 if $n == (0 + '-inf');
+ return 0+'inf' if $n == (0 + 'inf');
return Math::Prime::Util::PP::ExponentialIntegral($n) if defined
$bignum::VERSION || ref($n) eq 'Math::BigFloat';
return Math::Prime::Util::PP::ExponentialIntegral($n) if !$_Config{'xs'};
@@ -1424,6 +1430,8 @@ sub ExponentialIntegral {
sub LogarithmicIntegral {
my($n) = @_;
return 0 if $n == 0;
+ return 0+'-inf' if $n == 1;
+ return 0+'inf' if $n == (0 + 'inf');
croak("Invalid input to LogarithmicIntegral: x must be >= 0") if $n <= 0;
if ( defined $bignum::VERSION || ref($n) eq 'Math::BigFloat' ) {
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index d047d9f..5689490 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1296,13 +1296,14 @@ my $_const_li2 =
1.045163780117492784844588889194613136522615578151;
sub ExponentialIntegral {
my($x, $tol) = @_;
+ return 0+'-inf' if $x == 0;
+ return 0 if $x == (0 + '-inf');
+ return 0+'inf' if $x == (0 + 'inf');
$tol = 1e-16 unless defined $tol;
my $sum = 0.0;
my($y, $t);
my $c = 0.0;
- croak "Invalid input to ExponentialIntegral: x must be != 0" if $x == 0;
-
$x = new Math::BigFloat "$x" if defined $bignum::VERSION && ref($x) ne
'Math::BigFloat';
my $val; # The result from one of the four methods
@@ -1353,12 +1354,11 @@ sub ExponentialIntegral {
} else {
# Asymptotic divergent series
my $invx = 1.0 / $x;
- $val = exp($x) * $invx;
- $sum = 1.0;
- my $term = 1.0;
- for my $n (1 .. 200) {
+ my $term = $invx;
+ $sum = 1.0 + $term;
+ for my $n (2 .. 200) {
my $last_term = $term;
- $term *= $n*$invx;
+ $term *= $n * $invx;
last if $term < $tol;
if ($term < $last_term) {
$y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
@@ -1367,7 +1367,7 @@ sub ExponentialIntegral {
last;
}
}
- $val *= $sum;
+ $val = exp($x) * $invx * $sum;
}
$val;
}
@@ -1375,7 +1375,8 @@ sub ExponentialIntegral {
sub LogarithmicIntegral {
my($x) = @_;
return 0 if $x == 0;
- return 0+(-Infinity) if $x == 1;
+ return 0+'-inf' if $x == 1;
+ return 0+'inf' if $x == (0 + 'inf');
return $_const_li2 if $x == 2;
croak "Invalid input to LogarithmicIntegral: x must be > 0" if $x <= 0;
@@ -1384,12 +1385,13 @@ sub LogarithmicIntegral {
# Do divergent series here for big inputs. Common for big pc approximations.
if ($x > 1e16) {
- my $tol = 1e-16;
+ my $tol = 1e-20;
my $invx = 1.0 / $logx;
- my $val = $x * $invx;
- my $sum = 1.0;
- my $term = 1.0;
- for my $n (1 .. 200) {
+ # n = 0 => 0!/(logx)^0 = 1/1 = 1
+ # n = 1 => 1!/(logx)^1 = 1/logx
+ my $term = $invx;
+ my $sum = 1.0 + $term;
+ for my $n (2 .. 200) {
my $last_term = $term;
$term *= $n * $invx;
last if $term < $tol;
@@ -1400,9 +1402,26 @@ sub LogarithmicIntegral {
last;
}
}
- $val *= $sum;
+ my $val = $x * $invx * $sum;
return $val;
}
+ # Convergent series.
+ if ($x >= 1) {
+ my $tol = 1e-20;
+ my $fact_n = 1.0;
+ my $nfac = 1.0;
+ my $sum = 0.0;
+ for my $n (1 .. 200) {
+ $fact_n *= $logx/$n;
+ my $term = $fact_n / $n;
+ $sum += $term;
+ last if $term < $tol;
+ }
+ my $eulerconst = (ref($x) eq 'Math::BigFloat') ?
Math::BigFloat->new('0.577215664901532860606512090082402431042159335939923598805767')
: $_const_euler;
+ my $val = $eulerconst + log($logx) + $sum;
+ return $val;
+ }
+
ExponentialIntegral($logx);
}
diff --git a/t/16-randomprime.t b/t/16-randomprime.t
index ad493ba..8292294 100644
--- a/t/16-randomprime.t
+++ b/t/16-randomprime.t
@@ -150,6 +150,14 @@ foreach my $digits ( @random_ndigit_tests ) {
"$digits-digit random prime is in range and prime");
}
+SKIP: {
+ if ($use64 && $broken64) {
+ my $num_nbit_tests = scalar @random_nbit_tests;
+ @random_nbit_tests = grep { $_ < 50 } @random_nbit_tests;
+ my $nskip = $num_nbit_tests - scalar @random_nbit_tests;
+ skip "Skipping random 50+ bit primes on broken 64-bit Perl", 2*$nskip;
+ }
+}
foreach my $bits ( @random_nbit_tests ) {
check_bits( random_nbit_prime($bits), $bits, "nbit" );
check_bits( random_maurer_prime($bits), $bits, "Maurer" );
diff --git a/t/18-functions.t b/t/18-functions.t
index a6415eb..f597652 100644
--- a/t/18-functions.t
+++ b/t/18-functions.t
@@ -8,10 +8,8 @@ use Math::Prime::Util qw/prime_count ExponentialIntegral
LogarithmicIntegral Rie
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
-plan tests => 4 + 1 + 1 + 16 + 11 + 9;
+plan tests => 3 + 6 + 1 + 16 + 11 + 9;
-eval { ExponentialIntegral(0); };
-like($@, qr/invalid/i, "Ei(0) is invalid");
eval { LogarithmicIntegral(-1); };
like($@, qr/invalid/i, "li(-1) is invalid");
eval { RiemannR(0); };
@@ -19,7 +17,13 @@ like($@, qr/invalid/i, "R(0) is invalid");
eval { RiemannR(-1); };
like($@, qr/invalid/i, "R(-1) is invalid");
-cmp_ok( LogarithmicIntegral(1), '<=', 0 - (~0 * ~0), "li(1) is -infinity" );
+cmp_ok( ExponentialIntegral(0), '<=', 0 - (~0 * ~0), "Ei(0) is -infinity"
);
+cmp_ok( ExponentialIntegral('-inf'),'==', 0, "Ei(-inf) is 0" );
+cmp_ok( ExponentialIntegral('inf'), '>=', 0 + (~0 * ~0), "Ei(inf) is infinity"
);
+
+cmp_ok( LogarithmicIntegral(0), '==', 0, "li(0) is 0" );
+cmp_ok( LogarithmicIntegral(1), '<=', 0 - (~0 * ~0), "li(1) is -infinity" );
+cmp_ok( LogarithmicIntegral('inf'),'>=', 0 + (~0 * ~0), "li(inf) is infinity"
);
# Example used in Math::Cephes
cmp_closeto( ExponentialIntegral(2.2), 5.732614700, 1e-06, "Ei(2.2)");
--
Alioth's /usr/local/bin/git-commit-notice on
/srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
