[libmath-prime-util-perl] 01/59: PP sieve benchmark, update tests, MR always returns 0 for even input
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:50 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.10
in repository libmath-prime-util-perl.
commit 8c1a02c2ef678b37751d706e44e22b2663a9e63c
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Jun 25 19:18:53 2012 -0600
PP sieve benchmark, update tests, MR always returns 0 for even input
---
Changes | 14 +-
MANIFEST | 1 +
TODO | 6 +-
XS.xs | 1 +
examples/bench-sieve-pp.pl | 274 ++++++++++++++++++++++++++++
lib/Math/Prime/Util.pm | 29 ++-
lib/Math/Prime/Util/PP.pm | 1 +
t/02-can.t | 3 +-
t/13-primecount.t | 88 +++++----
t/14-nthprime.t | 48 +++--
t/17-pseudoprime.t | 121 +++++++------
t/80-pp.t | 438 +++++++++++++++++++++++++++++++++++++++++++++
12 files changed, 900 insertions(+), 124 deletions(-)
diff --git a/Changes b/Changes
index 4d444fa..b13078c 100644
--- a/Changes
+++ b/Changes
@@ -1,14 +1,22 @@
Revision history for Perl extension Math::Prime::Util.
+0.10 30 June 2012
+ - Make miller_rabin return 0 if given even number.
+ - Add tests for PP.
+ - Cleanup of some tests.
+
0.09 25 June 2012
- - Pure Perl code. Passes all tests, but 1 to 120x slower. Test suite
- as a whole is 38x slower.
+ - Pure Perl code added. Passes all tests. Used only if the XSLoader
+ fails. It's 1-120x slower than the C code. When forced to use the
+ PP code, the test suite is 38x slower on 64-bit, 16x slower on 32-bit
+ (in 64-bit, the test suite runs some large numbers through routines
+ like prime_count and nth_prime that are much faster in C).
- Modifications to threading test:
- some machines were failing because they use non-TS rand. Fix by
making our own rand.
- Win32 was failing because of unique threading issues. It barfs
if you free memory on a different thread than allocated it.
- - is_prime could return 1 in some cases. Fixed.
+ - is_prime could return 1 in some cases. Fixed to only return 0 or 2.
0.08 22 June 2012
- Added thread safety and tested good concurrency.
diff --git a/MANIFEST b/MANIFEST
index 49ad6b1..c633711 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -53,6 +53,7 @@ t/30-relations.t
t/31-threading.t
t/50-factoring.t
t/51-primearray.t
+t/80-pp.t
t/90-release-perlcritic.t
t/91-release-pod-syntax.t
t/92-release-pod-coverage.t
diff --git a/TODO b/TODO
index b78481f..da5c284 100644
--- a/TODO
+++ b/TODO
@@ -9,7 +9,9 @@
- Add test to check maxbits in compiled library vs. Perl
-- Pure perl implementations
+- Figure out documentation solution for PP.pm
+
+- Is the current PP.pm setup the way we want to do it?
- input validation (in XS, or do we need to make Perl wrappers for everything?)
We can do inpuut val in XS by looking at the NV. But I think long term we'll
@@ -24,5 +26,3 @@
- Move .c / .h files into separate directory.
version does it in a painful way. Something simpler to be had?
-
-- threading and Windows. Ugh.
diff --git a/XS.xs b/XS.xs
index 1e0e0f2..d7f2f32 100644
--- a/XS.xs
+++ b/XS.xs
@@ -351,6 +351,7 @@ miller_rabin(IN UV n, ...)
croak("No bases given to miller_rabin");
if ( (n == 0) || (n == 1) ) XSRETURN_IV(0); /* 0 and 1 are composite */
if ( (n == 2) || (n == 3) ) XSRETURN_IV(2); /* 2 and 3 are prime */
+ if (( n % 2 ) == 0) XSRETURN_IV(0); /* MR works with odd n */
while (c < items) {
int b = 0;
while (c < items) {
diff --git a/examples/bench-sieve-pp.pl b/examples/bench-sieve-pp.pl
new file mode 100644
index 0000000..6eb84d8
--- /dev/null
+++ b/examples/bench-sieve-pp.pl
@@ -0,0 +1,274 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Benchmark qw/:all/;
+
+my $upper = shift || 8192;
+my $count = shift || -1;
+my $countarg;
+
+my $pc_subs = {
+ "Scriptol" => sub { scriptol($countarg) },
+ "Shootout" => sub { shootout($countarg) },
+ "In Many" => sub { inmany($countarg) },
+ "Rosetta 1" => sub { rosetta1($countarg) },
+ "Rosetta 2" => sub { rosetta2($countarg) },
+ "Rosetta 3" => sub { rosetta3($countarg) },
+ "Rosetta 4" => sub { rosetta4($countarg) },
+ "DJ Vec" => sub { dj1($countarg) },
+ "DJ Array" => sub { dj2($countarg) },
+ "DJ String1" => sub { dj3($countarg) },
+ "DJ String2" => sub { dj4($countarg) },
+};
+
+my %verify = (
+ 10 => 4,
+ 100 => 25,
+ 1000 => 168,
+ 10000 => 1229,
+ 100000 => 9592,
+);
+
+# Verify
+while (my($name, $sub) = each (%$pc_subs)) {
+ while (my($n, $pin) = each (%verify)) {
+ $countarg = $n;
+ my $picount = $sub->();
+ die "$name ($n) = $picount, should be $pin" unless $picount == $pin;
+ }
+}
+print "Done with verification, starting benchmark\n";
+
+$countarg = $upper;
+cmpthese($count, $pc_subs);
+
+
+
+# www.scriptol.com/programming/sieve.php
+sub scriptol {
+ my($max) = @_;
+ return 0 if $upper < 2;
+ return 1 if $upper < 3;
+
+ my @flags = (0 .. $max);
+ for my $i (2 .. int(sqrt($max)) + 1)
+ {
+ next unless defined $flags[$i];
+ for (my $k=$i+$i; $k <= $max; $k+=$i)
+ {
+ undef $flags[$k];
+ }
+ }
+ my $count = 0;
+ for my $j (2 .. $max) {
+ $count++ if defined $flags[$j];
+ }
+ $count;
+}
+
+# http://dada.perl.it/shootout/sieve.perl.html
+sub shootout {
+ my($max) = @_;
+ return 0 if $upper < 2;
+ return 1 if $upper < 3;
+
+ my $count = 0;
+ my @flags = (0 .. $max);
+ for my $i (2 .. $max) {
+ next unless defined $flags[$i];
+ for (my $k=$i+$i; $k <= $max; $k+=$i) {
+ undef $flags[$k];
+ }
+ $count++;
+ }
+ $count;
+}
+
+# http://c2.com/cgi/wiki?SieveOfEratosthenesInManyProgrammingLanguages
+sub inmany {
+ my($max) = @_;
+ return 0 if $upper < 2;
+ return 1 if $upper < 3;
+
+ my @c;
+ for(my $t=3; $t*$t<$max; $t+=2) {
+ if (!$c[$t]) {
+ for(my $s=$t*$t; $s<$max; $s+=$t*2) { $c[$s]++ }
+ }
+ }
+ my $count = 1;
+ for(my $t=3; $t<$max; $t+=2) {
+ $c[$t] || $count++;
+ }
+ $count;
+}
+
+# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl
+sub rosetta1 {
+ my($max) = @_;
+ return 0 if $upper < 2;
+ return 1 if $upper < 3;
+
+ my $count = 0; #my @primes;
+ my @tested = (1);
+ my $j = 1;
+ while ($j < $max) {
+ next if $tested[$j++];
+ $count++; #push @primes, $j;
+ for (my $k= $j; $k <= $max; $k+=$j) {
+ $tested[$k-1]= 1;
+ }
+ }
+ $count; #scalar @primes;
+}
+
+# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl
+sub rosetta2 {
+ my($max) = @_;
+ return 0 if $upper < 2;
+ return 1 if $upper < 3;
+
+ my $count = 0; #my @primes;
+ my $nonPrimes = '';
+ foreach my $p (2 .. $max) {
+ unless (vec($nonPrimes, $p, 1)) {
+ for (my $i = $p * $p; $i <= $max; $i += $p) {
+ vec($nonPrimes, $i, 1) = 1;
+ }
+ $count++; #push @primes, $p;
+ }
+ }
+ $count; #scalar @primes;
+}
+
+# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl
+sub rosetta3 {
+ my($max) = @_;
+ return 0 if $upper < 2;
+ return 1 if $upper < 3;
+
+ my $i;
+ my @s;
+ my $count = scalar
+ grep { not $s[ $i = $_ ] and do
+ { $s[ $i += $_ ]++ while $i <= $max; 1 }
+ } 2 .. $max;
+ $count; #scalar @primes;
+}
+
+# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl
+sub rosetta4 {
+ my($max) = @_;
+ return 0 if $upper < 2;
+ return 1 if $upper < 3;
+
+ my $i;
+ my $s = '';
+ my $count = scalar
+ grep { not vec $s, $i = $_, 1 and do
+ { (vec $s, $i += $_, 1) = 1 while $i <= $max; 1 }
+ } 2 .. $max;
+ $count; #scalar @primes;
+}
+
+sub dj1 {
+ my($end) = @_;
+ return 0 if $upper < 2;
+ return 1 if $upper < 3;
+
+ # vector
+ my $sieve = '';
+ my $n = 3;
+ while ( ($n*$n) <= $end ) {
+ my $s = $n*$n;
+ while ($s <= $end) {
+ vec($sieve, $s >> 1, 1) = 1;
+ $s += 2*$n;
+ }
+ do { $n += 2 } while vec($sieve, $n >> 1, 1) != 0;
+ }
+ my $count = 1;
+ $n = 3;
+ while ($n <= $end) {
+ $count++ if !vec($sieve, $n >> 1, 1);
+ $n += 2;
+ }
+ $count;
+}
+
+sub dj2 {
+ my($end) = @_;
+ return 0 if $upper < 2;
+ return 1 if $upper < 3;
+
+ # array
+ my @sieve;
+ my $n = 3;
+ while ( ($n*$n) <= $end ) {
+ my $s = $n*$n;
+ while ($s <= $end) {
+ $sieve[$s>>1] = 1;
+ $s += 2*$n;
+ }
+ do { $n += 2 } while $sieve[$n>>1];
+ }
+ my $count = 1;
+ $n = 3;
+ while ($n <= $end) {
+ $count++ if !$sieve[$n>>1];
+ $n += 2;
+ }
+ $count;
+}
+
+# ~2x faster than inmany, lots faster than the others. Only loses to dj4,
+# which is just this code with a presieve added.
+sub dj3 {
+ my($end) = @_;
+ return 0 if $upper < 2;
+ return 1 if $upper < 3;
+ $end-- if ($end & 1) == 0;
+
+ # string
+ my $sieve = '1' . '0' x ($end>>1);
+ my $n = 3;
+ while ( ($n*$n) <= $end ) {
+ my $s = $n*$n;
+ my $filter_s = $s >> 1;
+ my $filter_end = $end >> 1;
+ while ($filter_s <= $filter_end) {
+ substr($sieve, $filter_s, 1) = '1';
+ $filter_s += $n;
+ }
+ do { $n += 2 } while substr($sieve, $n>>1, 1);
+ }
+ my $count = 1 + $sieve =~ tr/0//;
+ $count;
+}
+
+# 2-3x faster than inmany, 6-7x faster than any of the other non-DJ methods.
+sub dj4 {
+ my($end) = @_;
+ return 0 if $upper < 2;
+ return 1 if $upper < 3;
+ $end-- if ($end & 1) == 0;
+
+ # string with prefill
+ my $whole = int( ($end>>1) / 15);
+ my $sieve = '100010010010110' . '011010010010110' x $whole;
+ substr($sieve, ($end>>1)+1) = '';
+ my $n = 7;
+ while ( ($n*$n) <= $end ) {
+ my $s = $n*$n;
+ my $filter_s = $s >> 1;
+ my $filter_end = $end >> 1;
+ while ($filter_s <= $filter_end) {
+ substr($sieve, $filter_s, 1) = '1';
+ $filter_s += $n;
+ }
+ do { $n += 2 } while substr($sieve, $n>>1, 1);
+ }
+ my $count = 1 + $sieve =~ tr/0//;
+ $count;
+}
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 6a8f3b7..bb5762a 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -526,15 +526,25 @@ polynomials, plus a correction term for small values to reduce the error.
Takes a positive number as input and one or more bases. The bases must be
between C<2> and C<n - 2>. Returns 2 is C<n> is definitely prime, 1 if C<n>
-is probably prime, and 0 if C<n> is definitely composite. Since this is
-just the Miller-Rabin test, a value of 2 is only returned for inputs of
-2 and 3, which are shortcut. If 0 is returned, then the number really is a
-composite. If 1 is returned, there is a good chance the number is prime
-(depending on the input and the bases), but we cannot be sure.
+is probably prime, and 0 if C<n> is definitely composite. A value of 2 will
+only be returned for the inputs of 2 and 3, which are shortcut.
+
+If 0 is returned, then the number really is a composite. If 1 is returned,
+then it is either a prime or a strong pseudoprime to all the given bases.
+Given enough bases, the chances become very, very strong that the number is
+actually prime.
This is usually used in combination with other tests to make either stronger
tests (e.g. the strong BPSW test) or deterministic results for numbers less
than some verified limit (such as the C<is_prob_prime> function in this module).
+However, given the chances of passing multiple bases, there are some math
+packages that just use multiple MR tests for primality testing.
+
+Even numbers other than 2 will always return 0 (composite). While the
+algorithm does run with even input, most sources define it only on odd input.
+Returning composite for all non-2 even input makes the function match most
+other implementations including L<Math::Primality>'s C<is_strong_pseudoprime>
+function.
=head2 is_prob_prime
@@ -585,9 +595,10 @@ Examples:
# Use Mersenne Twister
use Math::Random::MT::Auto qw/rand/;
- # Use my custom random function
+ # Use a custom random function
sub rand { ... }
+
=head2 random_ndigit_prime
say "My 4-digit prime number is: ", random_ndigit_prime(4);
@@ -774,6 +785,12 @@ approximation to the prime counting function.
Accuracy should be at least 14 digits.
+=head1 EXAMPLES
+
+Print pseudoprimes base 17:
+
+ perl -MMath::Prime::Util=:all -E 'my $n=$base|1; while(1) { print "$n " if miller_rabin($n,$base) && !is_prime($n); $n+=2; } BEGIN {$|=1; $base=17}'
+
=head1 LIMITATIONS
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 6f8e94e..d49901d 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -603,6 +603,7 @@ sub miller_rabin {
return 0 if ($n == 0) || ($n == 1);
return 2 if ($n == 2) || ($n == 3);
+ return 0 if ($n % 2) == 0;
# I was using bignum here for a while, but doing "$a ** $d" with a
# big $d is **ridiculously** slow. It's thousands of times faster
diff --git a/t/02-can.t b/t/02-can.t
index f3e0aea..c2dd783 100644
--- a/t/02-can.t
+++ b/t/02-can.t
@@ -13,6 +13,7 @@ my @functions = qw(
prime_count prime_count_lower prime_count_upper prime_count_approx
nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
random_prime random_ndigit_prime
- factor
+ factor all_factors
+ ExponentialIntegral LogarithmicIntegral RiemannR
);
can_ok( 'Math::Prime::Util', @functions);
diff --git a/t/13-primecount.t b/t/13-primecount.t
index c2860c3..a20103e 100644
--- a/t/13-primecount.t
+++ b/t/13-primecount.t
@@ -8,10 +8,6 @@ use Math::Prime::Util qw/prime_count prime_count_lower prime_count_upper prime_c
my $use64 = Math::Prime::Util::_maxbits > 32;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
-plan tests => 1 + 15*3 + 9 + 6*$extra + 3*18*$use64 + 11 + 6*$use64;
-
-ok( eval { prime_count(13); 1; }, "prime_count in void context");
-
# Powers of 2: http://oeis.org/A007053/b007053.txt
# Powers of 10: http://oeis.org/A006880/b006880.txt
my %pivals32 = (
@@ -51,25 +47,9 @@ my %pivals64 = (
1152921504606846975 => 28423094496953330,
18446744073709551615 => 425656284035217743,
);
-while (my($n, $pin) = each (%pivals32)) {
- cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
- cmp_ok( prime_count_lower($n), '<=', $pin, "Pi($n) >= lower estimate" );
- if ( ($n <= 2000000) || $extra ) {
- is( prime_count($n), $pin, "Pi($n) = $pin" );
- }
- my $approx_range = abs($pin - prime_count_approx($n));
- my $range_limit = ($n <= 100000000) ? 100 : 500;
- cmp_ok( $approx_range, '<=', $range_limit, "prime_count_approx($n) within $range_limit");
-}
-if ($use64) {
- while (my($n, $pin) = each (%pivals64)) {
- cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
- cmp_ok( prime_count_lower($n), '<=', $pin, "Pi($n) >= lower estimate" );
- my $approx = prime_count_approx($n);
- my $percent_limit = 0.0005;
- cmp_ok( abs($pin - $approx) / $pin, '<=', $percent_limit/100.0, "prime_count approx($n) within $percent_limit\% of Pi($n)");
- }
-}
+my %pivals_small = map { $_ => $pivals32{$_} }
+ grep { ($_ <= 2000000) || $extra }
+ keys %pivals32;
# ./primesieve 1e10 -o2**32 -c1
# ./primesieve 24689 7973249 -c1
@@ -93,31 +73,63 @@ my %intervals = (
"127976334671 +467" => 1,
"127976334672 +466" => 0,
);
+delete @intervals{ grep { (parse_range($_))[1] > ~0 } keys %intervals };
+
+plan tests => 0 + 1
+ + 3*scalar(keys %pivals32)
+ + scalar(keys %pivals_small)
+ + $use64 * scalar(keys %pivals64)
+ + scalar(keys %intervals);
+
+ok( eval { prime_count(13); 1; }, "prime_count in void context");
+
+while (my($n, $pin) = each (%pivals32)) {
+ cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
+ cmp_ok( prime_count_lower($n), '<=', $pin, "Pi($n) >= lower estimate" );
+ my $approx_range = abs($pin - prime_count_approx($n));
+ my $range_limit = ($n <= 100000000) ? 100 : 500;
+ cmp_ok( $approx_range, '<=', $range_limit, "prime_count_approx($n) within $range_limit");
+}
+while (my($n, $pin) = each (%pivals_small)) {
+ is( prime_count($n), $pin, "Pi($n) = $pin" );
+}
+if ($use64) {
+ while (my($n, $pin) = each (%pivals64)) {
+ cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
+ cmp_ok( prime_count_lower($n), '<=', $pin, "Pi($n) >= lower estimate" );
+ my $approx = prime_count_approx($n);
+ my $percent_limit = 0.0005;
+ cmp_ok( abs($pin - $approx) / $pin, '<=', $percent_limit/100.0, "prime_count approx($n) within $percent_limit\% of Pi($n)");
+ }
+}
while (my($range, $expect) = each (%intervals)) {
+ my($low,$high) = parse_range($range);
+ is( prime_count($low,$high), $expect, "prime_count($range) = $expect");
+}
+
+sub parse_range {
+ my($range) = @_;
my($low,$high);
+ my $fixnum = sub {
+ my $nstr = shift;
+ $nstr =~ s/^(\d+)e(\d+)$/$1*(10**$2)/e;
+ $nstr =~ s/^(\d+)\*\*(\d+)$/$1**$2/e;
+ die "Unknown string in test" unless $nstr =~ /^\d+$/;
+ $nstr;
+ };
if ($range =~ /(\S+)\s+to\s+(\S+)/) {
- $low = fixnum($1);
- $high = fixnum($2);
+ $low = $fixnum->($1);
+ $high = $fixnum->($2);
} elsif ($range =~ /(\S+)\s*\+\s*(\S+)/) {
- $low = fixnum($1);
- $high = $low + fixnum($2);
+ $low = $fixnum->($1);
+ $high = $low + $fixnum->($2);
} else {
die "Can't parse test data";
}
- next if $high > ~0;
- is( prime_count($low,$high), $expect, "prime_count($range) = $expect");
+ ($low,$high);
}
-exit(0);
-sub fixnum {
- my $nstr = shift;
- $nstr =~ s/^(\d+)e(\d+)$/$1*(10**$2)/e;
- $nstr =~ s/^(\d+)\*\*(\d+)$/$1**$2/e;
- die "Unknown string in test" unless $nstr =~ /^\d+$/;
- $nstr;
-}
-
# TODO: intervals. From primesieve:
# 155428406, // prime count 2^32 interval starting at 10^12
# 143482916, // prime count 2^32 interval starting at 10^13
diff --git a/t/14-nthprime.t b/t/14-nthprime.t
index e777ca9..19f660d 100644
--- a/t/14-nthprime.t
+++ b/t/14-nthprime.t
@@ -12,8 +12,6 @@ my $broken64 = (18446744073709550592 == ~0);
my $nsmallprimes = 1000;
my $nth_small_prime = 7919; # nth_prime(1000)
-plan tests => 7*2 + $nsmallprimes+1 + 9*3 + 7 + ($extra ? 9 : 7) + ($use64 ? 9*3 : 0);
-
my %pivals32 = (
1 => 0,
10 => 4,
@@ -24,19 +22,6 @@ my %pivals32 = (
1000000 => 78498,
);
-while (my($n, $pin) = each (%pivals32)) {
- my $next = $pin+1;
- cmp_ok( nth_prime($pin), '<=', $n, "nth_prime($pin) <= $n");
- cmp_ok( nth_prime($next), '>=', $n, "nth_prime($next) >= $n");
-}
-
-my @small_primes = (0);
-push @small_primes, @{primes($nth_small_prime)};
-foreach my $n (0 .. $nsmallprimes) {
- is(nth_prime($n), $small_primes[$n], "The ${n}th prime is $small_primes[$n]");
-}
-
-
# Powers of 10: http://oeis.org/A006988/b006988.txt
# Powers of 2: http://oeis.org/A033844/b033844.txt
my %nthprimes32 = (
@@ -61,19 +46,44 @@ my %nthprimes64 = (
10000000000000000 => 394906913903735329,
100000000000000000 => 4185296581467695669,
);
+my %nthprimes_small = map { $_ => $nthprimes32{$_} }
+ grep { ($_ <= 2000000) || $extra }
+ keys %nthprimes32;
+
+my @small_primes = (0, @{primes($nth_small_prime)});
+
+#plan tests => 7*2 + $nsmallprimes+1 + 9*3 + 7 + ($extra ? 9 : 7) + ($use64 ? 9*3 : 0);
+
+plan tests => 0 + 2*scalar(keys %pivals32)
+ + scalar @small_primes
+ + 3*scalar(keys %nthprimes32)
+ + scalar(keys %nthprimes_small)
+ + $use64 * 3 * scalar(keys %nthprimes64)
+ + 4
+ + 3;
+
+
+while (my($n, $pin) = each (%pivals32)) {
+ my $next = $pin+1;
+ cmp_ok( nth_prime($pin), '<=', $n, "nth_prime($pin) <= $n");
+ cmp_ok( nth_prime($next), '>=', $n, "nth_prime($next) >= $n");
+}
+
+foreach my $n (0 .. $#small_primes) {
+ is(nth_prime($n), $small_primes[$n], "The ${n}th prime is $small_primes[$n]");
+}
while (my($n, $nth) = each (%nthprimes32)) {
cmp_ok( nth_prime_upper($n), '>=', $nth, "nth_prime($n) <= upper estimate" );
cmp_ok( nth_prime_lower($n), '<=', $nth, "nth_prime($n) >= lower estimate" );
- if ( ($n <= 2000000) || $extra ) {
- is( nth_prime($n), $nth, "nth_prime($n) = $nth" );
- }
-
my $approx = nth_prime_approx($n);
my $percent_limit = ($n >= 775) ? 1 : 2;
cmp_ok( abs($nth - $approx) / $nth, '<=', $percent_limit/100.0, "nth_prime_approx($n) = $approx within $percent_limit\% of $nth");
}
+while (my($n, $nth) = each (%nthprimes_small)) {
+ is( nth_prime($n), $nth, "nth_prime($n) = $nth" );
+}
if ($use64) {
while (my($n, $nth) = each (%nthprimes64)) {
diff --git a/t/17-pseudoprime.t b/t/17-pseudoprime.t
index e78d15e..ae07392 100644
--- a/t/17-pseudoprime.t
+++ b/t/17-pseudoprime.t
@@ -8,7 +8,68 @@ use Math::Prime::Util qw/is_prime miller_rabin/;
my $use64 = Math::Prime::Util::_maxbits > 32;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
-plan tests => 3 + 4 + 295 + 4 + 4*$use64 + 1 + 1*$extra + 161;
+# small primes
+my @sp = qw/2 3 5 7 11 13 17 19 23 29 31 37/;
+# strong pseudoprimes for all prime bases 2 .. pn
+my @phis = qw/2047 1373653 25326001 3215031751 2152302898747 3474749660383 341550071728321 341550071728321/;
+$#phis = 3 unless $use64;
+
+# pseudoprimes from 2-100k for each prime base
+#
+# Using a different codebase to get reference values:
+# perl -E 'use Math::Primality ":all"; for (2 .. 100000) { print "$_ " if is_strong_pseudoprime($_,17) && !is_prime($_); } print "\n"'
+#
+# With large values, one of:
+#
+# perl -MMath::Primality=:all -E 'my $_=$base|1; while(1) {print "$_ " if is_strong_pseudoprime($_,$base) && !is_prime($_); $_+=2; } print "\n"; BEGIN {$|=1; $base=553174392}'
+#
+# perl -MMath::Primality=is_strong_pseudoprime -MMath::Prime::Util=is_prime -E 'my $_=$base|1; while(1) {print "$_ " if is_strong_pseudoprime($_,$base) && !is_prime($_); $_+=2; } print "\n"; BEGIN {$|=1; $base=553174392}'
+#
+# ~30x faster than Math::Primality:
+# perl -MMath::Prime::Util=:all -E 'my $_=$base|1; while(1) {print "$_ " if miller_rabin($_,$base) && !is_prime($_); $_+=2; } print "\n"; BEGIN {$|=1; $base=553174392}'
+
+my %pseudoprimes = (
+ 2 => [ qw/2047 3277 4033 4681 8321 15841 29341 42799 49141 52633 65281 74665 80581 85489 88357 90751/ ],
+ 3 => [ qw/121 703 1891 3281 8401 8911 10585 12403 16531 18721 19345 23521 31621 44287 47197 55969 63139 74593 79003 82513 87913 88573 97567/ ],
+ 5 => [ qw/781 1541 5461 5611 7813 13021 14981 15751 24211 25351 29539 38081 40501 44801 53971 79381/ ],
+ 7 => [ qw/25 325 703 2101 2353 4525 11041 14089 20197 29857 29891 39331 49241 58825 64681 76627 78937 79381 87673 88399 88831/ ],
+ 11 => [ qw/133 793 2047 4577 5041 12403 13333 14521 17711 23377 43213 43739 47611 48283 49601 50737 50997 56057 58969 68137 74089 85879 86347 87913 88831/ ],
+ 13 => [ qw/85 1099 5149 7107 8911 9637 13019 14491 17803 19757 20881 22177 23521 26521 35371 44173 45629 54097 56033 57205 75241 83333 85285 86347/ ],
+ 17 => [ qw/9 91 145 781 1111 2821 4033 4187 5365 5833 6697 7171 15805 19729 21781 22791 24211 26245 31621 33001 33227 34441 35371 38081 42127 49771 71071 74665 77293 78881 88831 96433 97921 98671/ ],
+ 19 => [ qw/9 49 169 343 1849 2353 2701 4033 4681 6541 6697 7957 9997 12403 13213 13747 15251 16531 18769 19729 24761 30589 31621 31861 32477 41003 49771 63139 64681 65161 66421 68257 73555 96049/ ],
+ 23 => [ qw/169 265 553 1271 2701 4033 4371 4681 6533 6541 7957 8321 8651 8911 9805 14981 18721 25201 31861 34133 44173 47611 47783 50737 57401 62849 82513 96049/ ],
+ 29 => [ qw/15 91 341 469 871 2257 4371 4411 5149 6097 8401 11581 12431 15577 16471 19093 25681 28009 29539 31417 33001 48133 49141 54913 79003/ ],
+ 31 => [ qw/15 49 133 481 931 6241 8911 9131 10963 11041 14191 17767 29341 56033 58969 68251 79003 83333 87061 88183/ ],
+ 37 => [ qw/9 451 469 589 685 817 1333 3781 8905 9271 18631 19517 20591 25327 34237 45551 46981 47587 48133 59563 61337 68101 68251 73633 79381 79501 83333 84151 96727/ ],
+ 61 => [ qw/217 341 1261 2701 3661 6541 6697 7613 13213 16213 22177 23653 23959 31417 50117 61777 63139 67721 76301 77421 79381 80041/ ],
+ 73 => [ qw/205 259 533 1441 1921 2665 3439 5257 15457 23281 24617 26797 27787 28939 34219 39481 44671 45629 64681 67069 76429 79501 93521/ ],
+ 325 => [ qw/341 343 697 1141 2059 2149 3097 3537 4033 4681 4941 5833 6517 7987 8911 12403 12913 15043 16021 20017 22261 23221 24649 24929 31841 35371 38503 43213 44173 47197 50041 55909 56033 58969 59089 61337 65441 68823 72641 76793 78409 85879/ ],
+ 9375 => [ qw/11521 14689 17893 18361 20591 28093 32809 37969 44287 60701 70801 79957 88357 88831 94249 96247 99547/ ],
+ 28178 => [ qw/28179 29381 30353 34441 35371 37051 38503 43387 50557 51491 57553 79003 82801 83333 87249 88507 97921 99811/ ],
+ 75088 => [ qw/75089 79381 81317 91001 100101 111361 114211 136927 148289 169641 176661 191407 195649/ ],
+ 450775 => [ qw/465991 468931 485357 505441 536851 556421 578771 585631 586249 606361 631651 638731 641683 645679/ ],
+ 642735 => [ qw/653251 653333 663181 676651 714653 759277 794683 805141 844097 872191 874171 894671/ ],
+ 9780504 => [ qw/9780505 9784915 9826489 9882457 9974791 10017517 10018081 10084177 10188481 10247357 10267951 10392241 10427209 10511201/ ],
+ 203659041 => [ qw/204172939 204456793 206407057 206976001 207373483 209301121 210339397 211867969 212146507 212337217 212355793 214400629 214539841 215161459/ ],
+ 553174392 => [ qw/553174393 553945231 554494951 554892787 555429169 557058133 557163157 557165209 558966793 559407061 560291719 561008251 563947141/ ],
+ 1005905886 => [ qw/1005905887 1007713171 1008793699 1010415421 1010487061 1010836369 1012732873 1015269391 1016250247 1018405741 1020182041/ ],
+ 1340600841 => [ qw/1345289261 1345582981 1347743101 1348964401 1350371821 1353332417 1355646961 1357500901 1361675929 1364378203 1366346521 1367104639/ ],
+ 1795265022 => [ qw/1795265023 1797174457 1797741901 1804469753 1807751977 1808043283 1808205701 1813675681 1816462201 1817936371 1819050257/ ],
+ 3046413974 => [ qw/3046413975 3048698683 3051199817 3068572849 3069705673 3070556233 3079010071 3089940811 3090723901 3109299161 3110951251 3113625601/ ],
+ 3613982119 => [ qw/3626488471 3630467017 3643480501 3651840727 3653628247 3654142177 3672033223 3672036061 3675774019 3687246109 3690036017 3720856369/ ],
+);
+my $num_pseudoprimes = 0;
+foreach my $ppref (values %pseudoprimes) {
+ $num_pseudoprimes += scalar @$ppref;
+}
+
+
+plan tests => 0 + 3
+ + 4
+ + $num_pseudoprimes
+ + scalar @phis
+ + 1
+ + 1*$extra;
ok(!eval { miller_rabin(2047); }, "MR with no base fails");
ok(!eval { miller_rabin(2047,0); }, "MR base 0 fails");
@@ -19,39 +80,17 @@ is( miller_rabin(1, 2), 0, "MR with 0 shortcut composite");
is( miller_rabin(2, 2), 2, "MR with 2 shortcut prime");
is( miller_rabin(3, 2), 2, "MR with 3 shortcut prime");
-# small primes
-my @sp = qw/2 3 5 7 11 13 17 19 23 29 31 37/;
-# strong pseudoprimes for all prime bases 2 .. pn
-my @phis = qw/2047 1373653 25326001 3215031751 2152302898747 3474749660383 341550071728321 341550071728321/;
-
-# pseudoprimes from 2-100k for each prime base
-# perl -E 'use Math::Primality ":all"; for (2 .. 100000) { print "$_ " if is_strong_pseudoprime($_,17) && !is_prime($_); } print "\n"'
-my @psrp = (
- [ qw/2047 3277 4033 4681 8321 15841 29341 42799 49141 52633 65281 74665 80581 85489 88357 90751/ ],
- [ qw/121 703 1891 3281 8401 8911 10585 12403 16531 18721 19345 23521 31621 44287 47197 55969 63139 74593 79003 82513 87913 88573 97567/ ],
- [ qw/781 1541 5461 5611 7813 13021 14981 15751 24211 25351 29539 38081 40501 44801 53971 79381/ ],
- [ qw/25 325 703 2101 2353 4525 11041 14089 20197 29857 29891 39331 49241 58825 64681 76627 78937 79381 87673 88399 88831/ ],
- [ qw/133 793 2047 4577 5041 12403 13333 14521 17711 23377 43213 43739 47611 48283 49601 50737 50997 56057 58969 68137 74089 85879 86347 87913 88831/ ],
- [ qw/85 1099 5149 7107 8911 9637 13019 14491 17803 19757 20881 22177 23521 26521 35371 44173 45629 54097 56033 57205 75241 83333 85285 86347/ ],
- [ qw/9 91 145 781 1111 2821 4033 4187 5365 5833 6697 7171 15805 19729 21781 22791 24211 26245 31621 33001 33227 34441 35371 38081 42127 49771 71071 74665 77293 78881 88831 96433 97921 98671/ ],
- [ qw/9 49 169 343 1849 2353 2701 4033 4681 6541 6697 7957 9997 12403 13213 13747 15251 16531 18769 19729 24761 30589 31621 31861 32477 41003 49771 63139 64681 65161 66421 68257 73555 96049/ ],
- [ qw/169 265 553 1271 2701 4033 4371 4681 6533 6541 7957 8321 8651 8911 9805 14981 18721 25201 31861 34133 44173 47611 47783 50737 57401 62849 82513 96049/ ],
- [ qw/15 91 341 469 871 2257 4371 4411 5149 6097 8401 11581 12431 15577 16471 19093 25681 28009 29539 31417 33001 48133 49141 54913 79003/ ],
- [ qw/15 49 133 481 931 6241 8911 9131 10963 11041 14191 17767 29341 56033 58969 68251 79003 83333 87061 88183/ ],
- [ qw/9 451 469 589 685 817 1333 3781 8905 9271 18631 19517 20591 25327 34237 45551 46981 47587 48133 59563 61337 68101 68251 73633 79381 79501 83333 84151 96727/ ],
-);
# Check that each strong pseudoprime base b makes it through MR with that base
-my $bindex = 0;
-foreach my $base (@sp) {
- foreach my $p ( @{ $psrp[$bindex++] } ) {
+while (my($base, $ppref) = each (%pseudoprimes)) {
+ foreach my $p (@$ppref) {
ok(miller_rabin($p, $base), "Pseudoprime (base $base) $p passes MR");
}
}
# Check that phi_n makes passes MR with all prime bases < pn
-for my $phi (1 .. 8) {
- next if ($phi > 4) && (!$use64);
+for my $phi (1 .. scalar @phis) {
+ #next if ($phi > 4) && (!$use64);
ok( miller_rabin($phis[$phi-1], @sp[0 .. $phi-1]), "phi_$phi passes MR with first $phi primes");
}
@@ -66,7 +105,7 @@ for my $phi (1 .. 8) {
if (miller_rabin($_,2)) { $mr2fail = $_; last; }
}
}
- is($mr2fail, 0, "miller_rabin base 2 matches is_prime for 2-2046,2048-3276,3278-4032");
+ is($mr2fail, 0, "MR base 2 matches is_prime for 2-4032 (ex 2047,3277)");
}
# Verify MR base 2-3 for many small numbers (up to phi2)
@@ -81,29 +120,3 @@ if ($extra) {
}
is($mr2fail, 0, "miller_rabin bases 2,3 matches is_prime to 1,373,652");
}
-
-# More bases
-my @ebases = qw/61 73 325 9375 28178 75088 450775 642735/;
-my @epsrp = (
- [ qw/217 341 1261 2701 3661 6541 6697 7613 13213 16213 22177 23653 23959 31417 50117 61777 63139 67721 76301 77421 79381 80041/ ],
- [ qw/205 259 533 1441 1921 2665 3439 5257 15457 23281 24617 26797 27787 28939 34219 39481 44671 45629 64681 67069 76429 79501 93521/ ],
- [ qw/341 343 697 1141 2059 2149 3097 3537 4033 4681 4941 5833 6517 7987 8911 12403 12913 15043 16021 20017 22261 23221 24649 24929 31841 35371 38503 43213 44173 47197 50041 55909 56033 58969 59089 61337 65441 68823 72641 76793 78409 85879/ ],
- [ qw/11521 14689 17893 18361 20591 28093 32809 37969 44287 60701 70801 79957 88357 88831 94249 96247 99547/ ],
- [ qw/28179 29381 30353 34441 35371 37051 38503 43387 50557 51491 57553 79003 82801 83333 87249 88507 97921 99811/ ],
- [ qw/75089 79381 81317 91001 100101 111361 114211 136927 148289 169641 176661 191407 195649/ ],
- [ qw/465991 468931 485357 505441 536851 556421 578771 585631 586249 606361 631651 638731 641683 645679/ ],
- [ qw/653251 653333 663181 676651 714653 759277 794683 805141 844097 872191 874171 894671/ ],
-);
-
-# Check some of the extra bases we use
-$bindex = 0;
-foreach my $base (@ebases) {
- foreach my $p ( @{ $epsrp[$bindex++] } ) {
- ok(miller_rabin($p, $base), "Pseudoprime (base $base) $p passes MR");
- }
-}
-# TODO:
-# add tests for bases:
-# 1005905886, 1340600841, 553174392, 3046413974,
-# 203659041, 3613982119,
-# 9780504, 1795265022
diff --git a/t/80-pp.t b/t/80-pp.t
new file mode 100644
index 0000000..87c1a07
--- /dev/null
+++ b/t/80-pp.t
@@ -0,0 +1,438 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+# This is a subset of our tests. You really should run the whole test suite
+# on the PP code. What this will do is basic regression testing.
+
+use Test::More;
+my @small_primes = qw/
+2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71
+73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173
+179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281
+283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409
+419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541
+547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659
+661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809
+811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941
+947 953 967 971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049 1051
+1061 1063 1069
+/; # next prime is 1087
+my %small_primes;
+map { $small_primes{$_} = 1; } @small_primes;
+
+my @primes = qw/
+1129 1327 9551 15683 19609 31397 155921
+5 11 29 97 127 541 907 1151 1361 9587 15727 19661 31469 156007 360749
+370373 492227 1349651 1357333 2010881 4652507 17051887 20831533 47326913
+122164969 189695893 191913031
+/;
+
+my @composites = qw/
+0 1 4 6 8 9 10 12 14 15 16 18 20 21 22
+9 2047 1373653 25326001 3215031751
+561 1105 1729 2465 2821 6601 8911 10585 15841 29341 41041 46657 52633
+62745 63973 75361 101101 340561 488881 852841 1857241 6733693
+9439201 17236801 23382529 34657141 56052361 146843929
+341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371
+4681 5461 6601 7957 8321 52633 88357
+/;
+
+# pseudoprimes to various small prime bases
+my %pseudoprimes = (
+ 2 => [ qw/2047 3277 4033 4681 8321 15841 29341 42799 49141 52633 65281 74665 80581 85489 88357 90751/ ],
+ 3 => [ qw/121 703 1891 3281 8401 8911 10585 12403 16531 18721 19345 23521 31621 44287 47197 55969 63139 74593 79003 82513 87913 88573 97567/ ],
+ 5 => [ qw/781 1541 5461 5611 7813 13021 14981 15751 24211 25351 29539 38081 40501 44801 53971 79381/ ],
+ 7 => [ qw/25 325 703 2101 2353 4525 11041 14089 20197 29857 29891 39331 49241 58825 64681 76627 78937 79381 87673 88399 88831/ ],
+ 11 => [ qw/133 793 2047 4577 5041 12403 13333 14521 17711 23377 43213 43739 47611 48283 49601 50737 50997 56057 58969 68137 74089 85879 86347 87913 88831/ ],
+ 13 => [ qw/85 1099 5149 7107 8911 9637 13019 14491 17803 19757 20881 22177 23521 26521 35371 44173 45629 54097 56033 57205 75241 83333 85285 86347/ ],
+ 17 => [ qw/9 91 145 781 1111 2821 4033 4187 5365 5833 6697 7171 15805 19729 21781 22791 24211 26245 31621 33001 33227 34441 35371 38081 42127 49771 71071 74665 77293 78881 88831 96433 97921 98671/ ],
+ 19 => [ qw/9 49 169 343 1849 2353 2701 4033 4681 6541 6697 7957 9997 12403 13213 13747 15251 16531 18769 19729 24761 30589 31621 31861 32477 41003 49771 63139 64681 65161 66421 68257 73555 96049/ ],
+ 23 => [ qw/169 265 553 1271 2701 4033 4371 4681 6533 6541 7957 8321 8651 8911 9805 14981 18721 25201 31861 34133 44173 47611 47783 50737 57401 62849 82513 96049/ ],
+ 29 => [ qw/15 91 341 469 871 2257 4371 4411 5149 6097 8401 11581 12431 15577 16471 19093 25681 28009 29539 31417 33001 48133 49141 54913 79003/ ],
+ 31 => [ qw/15 49 133 481 931 6241 8911 9131 10963 11041 14191 17767 29341 56033 58969 68251 79003 83333 87061 88183/ ],
+ 37 => [ qw/9 451 469 589 685 817 1333 3781 8905 9271 18631 19517 20591 25327 34237 45551 46981 47587 48133 59563 61337 68101 68251 73633 79381 79501 83333 84151 96727/ ],
+ 61 => [ qw/217 341 1261 2701 3661 6541 6697 7613 13213 16213 22177 23653 23959 31417 50117 61777 63139 67721 76301 77421 79381 80041/ ],
+ 73 => [ qw/205 259 533 1441 1921 2665 3439 5257 15457 23281 24617 26797 27787 28939 34219 39481 44671 45629 64681 67069 76429 79501 93521/ ],
+);
+my $num_pseudoprimes = 0;
+foreach my $ppref (values %pseudoprimes) {
+ push @composites, @$ppref;
+ $num_pseudoprimes += scalar @$ppref;
+}
+{
+ my %uniq;
+ $uniq{$_}++ for (@composites);
+ @composites = sort {$a<=>$b} keys %uniq;
+}
+
+my %small_single = (
+ 0 => [],
+ 1 => [],
+ 2 => [2],
+ 3 => [2, 3],
+ 4 => [2, 3],
+ 5 => [2, 3, 5],
+ 6 => [2, 3, 5],
+ 7 => [2, 3, 5, 7],
+ 11 => [2, 3, 5, 7, 11],
+ 18 => [2, 3, 5, 7, 11, 13, 17],
+ 19 => [2, 3, 5, 7, 11, 13, 17, 19],
+ 20 => [2, 3, 5, 7, 11, 13, 17, 19],
+);
+
+my %small_range = (
+ "3 to 9" => [3,5,7],
+ "2 to 20" => [2,3,5,7,11,13,17,19],
+ "30 to 70" => [31,37,41,43,47,53,59,61,67],
+ "70 to 30" => [],
+ "20 to 2" => [],
+ "1 to 1" => [],
+ "2 to 2" => [2],
+ "3 to 3" => [3],
+ "2 to 3" => [2,3],
+ "2 to 5" => [2,3,5],
+ "3 to 6" => [3,5],
+ "3 to 7" => [3,5,7],
+ "4 to 8" => [5,7],
+ "2010733 to 2010881" => [2010733,2010881],
+ "2010734 to 2010880" => [],
+ "3088 to 3164" => [3089,3109,3119,3121,3137,3163],
+ "3089 to 3163" => [3089,3109,3119,3121,3137,3163],
+ "3090 to 3162" => [3109,3119,3121,3137],
+ "3842610773 to 3842611109" => [3842610773,3842611109],
+ "3842610774 to 3842611108" => [],
+);
+
+my %primegaps = (
+ 19609 => 52,
+ 360653 => 96,
+ 2010733 => 148,
+);
+
+my %pivals32 = (
+ 1 => 0,
+ 10 => 4,
+ 100 => 25,
+ 1000 => 168,
+ 10000 => 1229,
+ 100000 => 9592,
+ 1000000 => 78498,
+ 10000000 => 664579,
+ 100000000 => 5761455,
+ 1000000000 => 50847534,
+ 60067 => 6062,
+ 65535 => 6542,
+ 16777215 => 1077871,
+ 2147483647 => 105097565,
+ 4294967295 => 203280221,
+);
+my %pivals_small = map { $_ => $pivals32{$_} }
+ grep {$_ <= 2000000}
+ keys %pivals32;
+
+my %pi_intervals = (
+ "868396 to 9478505" => 563275,
+ "1118105 to 9961674" => 575195,
+ "24689 to 7973249" => 535368,
+ "1e10 +2**16" => 2821,
+ "17 to 13" => 0,
+ "3 to 17" => 6,
+ "4 to 17" => 5,
+ "4 to 16" => 4,
+ "191912783 +248" => 2,
+ "191912784 +247" => 1,
+ "191912783 +247" => 1,
+ "191912784 +246" => 0,
+);
+# Remove any entries where the high value is too large for us
+# ikegami++ for the delete from a hash slice idea
+delete @pi_intervals{ grep { (parse_range($_))[1] > ~0 } keys %pi_intervals };
+
+my %nthprimes32 = (
+ 1 => 2,
+ 10 => 29,
+ 100 => 541,
+ 1000 => 7919,
+ 10000 => 104729,
+ 100000 => 1299709,
+ 1000000 => 15485863,
+ 10000000 => 179424673,
+ 100000000 => 2038074743,
+);
+my %nthprimes_small = map { $_ => $nthprimes32{$_} }
+ grep {$_ <= 2000000}
+ keys %nthprimes32;
+
+my %eivals = (
+ -10 => -0.00000415696892968532438,
+ -0.5 => -0.55977359477616,
+ -0.1 => -1.8229239584193906660809,
+ -0.001 => -6.33153936413615,
+ -0.00001 => -10.9357198000436956,
+ -0.00000001 => -17.843465089050832587,
+ 0.693147180559945 => 1.0451637801174927848446, # log2
+ 1 => 1.8951178163559367554665,
+ 1.5 => 3.3012854491297978379574,
+ 2 => 4.9542343560018901633795,
+ 5 => 40.185275355803177455091,
+ 10 => 2492.2289762418777591384,
+ 12 => 14959.532666397528852292,
+ 20 => 25615652.664056588820481,
+ 40 => 6039718263611241.5783592,
+ 41 => 16006649143245041.110700,
+);
+my %livals = (
+ 0 => 0,
+ 1.01 => -4.0229586739299358695031,
+ 2 => 1.0451637801174927848446,
+ 10 => 6.1655995047872979375230,
+ 24 => 11.200315795232698830550,
+ 1000 => 177.60965799015222668764,
+ 100000 => 9629.8090010507982050343,
+ 100000000 => 5762209.3754480314675691,
+ 4294967295 => 203284081.95454158906409,
+ 10000000000 => 455055614.58662307560953,
+ 100000000000 => 4118066400.6216115150394,
+);
+my %rvals = (
+ 1.01 => 1.0060697180622924796117,
+ 2 => 1.5410090161871318832885,
+ 10 => 4.5645831410050902398658,
+ 1000 => 168.35944628116734806491,
+ 1000000 => 78527.399429127704858870,
+ 10000000 => 664667.44756474776798535,
+ 4294967295 => 203280697.51326064541983,
+ 10000000000 => 455050683.30684692446315,
+18446744073709551615 => 4.25656284014012122706963685602e17,
+);
+
+
+plan tests => 1 +
+ 2*(1087 + @primes + @composites) +
+ 3 + scalar(keys %small_single) + scalar(keys %small_range) +
+ 2*scalar(keys %primegaps) + 8 + 148 + 148 + 1 +
+ 3*scalar(keys %pivals32) + scalar(keys %pivals_small) + scalar(keys %pi_intervals) +
+ 2*scalar(keys %pivals_small) + 3*scalar(keys %nthprimes32) + scalar(keys %nthprimes_small) +
+ 4 + $num_pseudoprimes +
+ scalar(keys %eivals) + scalar(keys %livals) + scalar(keys %rvals) +
+ scalar @primes + 3*scalar @composites +
+ 0;
+
+use Math::Prime::Util qw/primes/;
+require_ok 'Math::Prime::Util::PP';
+ # This function skips some setup
+ undef *primes;
+ *primes = \&Math::Prime::Util::PP::primes;
+
+ *prime_count = \&Math::Prime::Util::PP::prime_count;
+ *prime_count_upper = \&Math::Prime::Util::PP::prime_count_upper;
+ *prime_count_lower = \&Math::Prime::Util::PP::prime_count_lower;
+ *prime_count_approx = \&Math::Prime::Util::PP::prime_count_approx;
+ *nth_prime = \&Math::Prime::Util::PP::nth_prime;
+ *nth_prime_upper = \&Math::Prime::Util::PP::nth_prime_upper;
+ *nth_prime_lower = \&Math::Prime::Util::PP::nth_prime_lower;
+ *nth_prime_approx = \&Math::Prime::Util::PP::nth_prime_approx;
+
+ *is_prime = \&Math::Prime::Util::PP::is_prime;
+ *next_prime = \&Math::Prime::Util::PP::next_prime;
+ *prev_prime = \&Math::Prime::Util::PP::prev_prime;
+
+ *miller_rabin = \&Math::Prime::Util::PP::miller_rabin;
+ *is_prob_prime = \&Math::Prime::Util::PP::is_prob_prime;
+
+ *factor = \&Math::Prime::Util::PP::factor;
+
+ *RiemannR = \&Math::Prime::Util::PP::RiemannR;
+ *LogarithmicIntegral = \&Math::Prime::Util::PP::LogarithmicIntegral;
+ *ExponentialIntegral = \&Math::Prime::Util::PP::ExponentialIntegral;
+
+###############################################################################
+
+foreach my $n (0 .. 1086) {
+ if (defined $small_primes{$n}) {
+ is( is_prime($n), 2, "$n is prime");
+ ok( is_prob_prime($n), "$n is probably prime");
+ } else {
+ ok(!is_prime($n), "$n is not prime");
+ ok(!is_prob_prime($n), "$n is not probably prime");
+ }
+}
+
+foreach my $n (@primes) {
+ is( is_prime($n), 2, "$n is prime" );
+ ok( is_prob_prime($n), "$n is probably prime");
+}
+foreach my $n (@composites) {
+ is( is_prime($n), 0, "$n is not prime" );
+ is( is_prob_prime($n), 0, "$n is not probably prime");
+}
+
+###############################################################################
+
+is_deeply( primes(1069), \@small_primes, "Primes between 0 and 1069" );
+is_deeply( primes(1070), \@small_primes, "Primes between 0 and 1070" );
+is_deeply( primes(1086), \@small_primes, "Primes between 0 and 1086" );
+
+while (my($high, $expect) = each (%small_single)) {
+ is_deeply( primes($high), $expect, "primes($high) should return [@{$expect}]")
+;
+}
+
+while (my($range, $expect) = each (%small_range)) {
+ my($low,$high) = $range =~ /(\d+) to (\d+)/;
+ is_deeply( primes($low, $high), $expect, "primes($low,$high) should return [@{$expect}]");
+}
+
+###############################################################################
+
+while (my($base, $range) = each (%primegaps)) {
+ is( next_prime($base), $base+$range, "next prime of $base is $base+$range" );
+ is( prev_prime($base+$range), $base, "prev prime of $base+$range is $base" );
+}
+
+is( next_prime(19608), 19609, "next prime of 19608 is 19609" );
+is( next_prime(19610), 19661, "next prime of 19610 is 19661" );
+is( next_prime(19660), 19661, "next prime of 19660 is 19661" );
+is( prev_prime(19662), 19661, "prev prime of 19662 is 19661" );
+is( prev_prime(19660), 19609, "prev prime of 19660 is 19609" );
+is( prev_prime(19610), 19609, "prev prime of 19610 is 19609" );
+
+is( prev_prime(2), 0, "Previous prime of 2 returns 0" );
+is( next_prime(~0-4), 0, "Next prime of ~0-4 returns 0" );
+
+foreach my $n (2010733 .. 2010880) {
+ is(next_prime($n), 2010881, "next_prime($n) == 2010881");
+}
+foreach my $n (2010734 .. 2010881) {
+ is(prev_prime($n), 2010733, "prev_prime($n) == 2010733");
+}
+# Similar test case to 2010870, where m=0 and next_prime is at m=1
+is(next_prime(1234567890), 1234567891, "next_prime(1234567890) == 1234567891)");
+
+###############################################################################
+
+while (my($n, $pin) = each (%pivals32)) {
+ cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
+ cmp_ok( prime_count_lower($n), '<=', $pin, "Pi($n) >= lower estimate" );
+ my $approx_range = abs($pin - prime_count_approx($n));
+ my $range_limit = ($n <= 100000000) ? 100 : 500;
+ cmp_ok( $approx_range, '<=', $range_limit, "prime_count_approx($n) within $range_limit");
+}
+while (my($n, $pin) = each (%pivals_small)) {
+ is( prime_count($n), $pin, "Pi($n) = $pin" );
+}
+
+while (my($range, $expect) = each (%pi_intervals)) {
+ my($low,$high) = parse_range($range);
+ is( prime_count($low,$high), $expect, "prime_count($range) = $expect");
+}
+
+###############################################################################
+
+while (my($n, $pin) = each (%pivals_small)) {
+ my $next = $pin+1;
+ cmp_ok( nth_prime($pin), '<=', $n, "nth_prime($pin) <= $n");
+ cmp_ok( nth_prime($next), '>=', $n, "nth_prime($next) >= $n");
+}
+
+while (my($n, $nth) = each (%nthprimes32)) {
+ cmp_ok( nth_prime_upper($n), '>=', $nth, "nth_prime($n) <= upper estimate" );
+ cmp_ok( nth_prime_lower($n), '<=', $nth, "nth_prime($n) >= lower estimate" );
+
+ my $approx = nth_prime_approx($n);
+ my $percent_limit = ($n >= 775) ? 1 : 2;
+ cmp_ok( abs($nth - $approx) / $nth, '<=', $percent_limit/100.0, "nth_prime_approx($n) = $approx within $percent_limit\% of $nth");
+}
+while (my($n, $nth) = each (%nthprimes_small)) {
+ is( nth_prime($n), $nth, "nth_prime($n) = $nth" );
+}
+
+###############################################################################
+
+is( miller_rabin(0, 2), 0, "MR with 0 shortcut composite");
+is( miller_rabin(1, 2), 0, "MR with 0 shortcut composite");
+is( miller_rabin(2, 2), 2, "MR with 2 shortcut prime");
+is( miller_rabin(3, 2), 2, "MR with 3 shortcut prime");
+
+while (my($base, $ppref) = each (%pseudoprimes)) {
+ foreach my $p (@$ppref) {
+ ok(miller_rabin($p, $base), "Pseudoprime (base $base) $p passes MR");
+ }
+}
+
+###############################################################################
+
+while (my($n, $ein) = each (%eivals)) {
+ cmp_closeto( ExponentialIntegral($n), $ein, 0.00000001 * abs($ein), "Ei($n) ~= $ein");
+}
+while (my($n, $lin) = each (%livals)) {
+ cmp_closeto( LogarithmicIntegral($n), $lin, 0.00000001 * abs($lin), "li($n) ~= $lin");
+}
+while (my($n, $rin) = each (%rvals)) {
+ cmp_closeto( RiemannR($n), $rin, 0.00000001 * abs($rin), "R($n) ~= $rin");
+}
+
+###############################################################################
+
+foreach my $n (@primes) {
+ my @f = factor($n);
+ is_deeply( \@f, [$n], "factor prime $n yields $n" );
+}
+foreach my $n (@composites) {
+ my @f = factor($n);
+ my $facstring = join(' * ', @f);
+
+ # Special case for 0 and 1
+ if ($n < 2) {
+ cmp_ok( scalar @f, '==', 1, "Factored small $n into itself" );
+ is( $f[0], $n, "$n = [ $facstring ]" );
+ ok( !is_prime($f[0]), "All factors [ $facstring ] of small $n are not prime" );
+ next;
+ }
+
+ # These are composites, so they should give us more than one factor
+ cmp_ok( scalar @f, '>=', 2, "Factored $n into multiple factors" );
+
+ # Do they multiply to the number?
+ my $product = 1; $product = int($product * $_) for @f;
+ is( $product, $n, "$n = [ $facstring ]" );
+
+ # Are they all prime?
+ my $isprime = 1; $isprime *= is_prime($_) for @f;
+ ok( $isprime, "All factors [ $facstring ] of $n are prime" );
+}
+
+###############################################################################
+
+sub parse_range {
+ my($range) = @_;
+ my($low,$high);
+ my $fixnum = sub {
+ my $nstr = shift;
+ $nstr =~ s/^(\d+)e(\d+)$/$1*(10**$2)/e;
+ $nstr =~ s/^(\d+)\*\*(\d+)$/$1**$2/e;
+ die "Unknown string in test" unless $nstr =~ /^\d+$/;
+ $nstr;
+ };
+ if ($range =~ /(\S+)\s+to\s+(\S+)/) {
+ $low = $fixnum->($1);
+ $high = $fixnum->($2);
+ } elsif ($range =~ /(\S+)\s*\+\s*(\S+)/) {
+ $low = $fixnum->($1);
+ $high = $low + $fixnum->($2);
+ } else {
+ die "Can't parse test data";
+ }
+ ($low,$high);
+}
+
+sub cmp_closeto {
+ my $got = shift;
+ my $expect = shift;
+ my $tolerance = shift;
+ my $message = shift;
+ cmp_ok( abs($got - $expect), '<=', $tolerance, $message );
+}
+
+done_testing();
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
More information about the Pkg-perl-cvs-commits
mailing list