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