[libmath-prime-util-perl] 02/02: Update for 0.07
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:31 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.07
in repository libmath-prime-util-perl.
commit 3f327f75d1566a3488a5205caf241c7ef76f3b44
Author: Dana Jacobsen <dana at acm.org>
Date: Sun Jun 17 15:11:47 2012 -0600
Update for 0.07
---
Changes | 5 +
MANIFEST | 7 +-
README | 2 +-
TODO | 2 +
XS.xs | 1 +
examples/README | 108 +++++++++++++++++++++
examples/bench-factor-extra.pl | 2 +-
examples/bench-primecount.pl | 38 ++++++++
.../{test-factoring.pl => test-factor-mpxs.pl} | 2 +-
examples/test-factor-yafu.pl | 105 ++++++++++++++++++++
examples/test-nextprime-yafu.pl | 82 ++++++++++++++++
examples/test-pcapprox.pl | 12 ++-
examples/test-primes-yafu.pl | 107 ++++++++++++++++++++
lib/Math/Prime/Util.pm | 4 +-
t/12-nextprime.t | 30 ++++--
util.c | 7 +-
16 files changed, 498 insertions(+), 16 deletions(-)
diff --git a/Changes b/Changes
index 0ae8163..7e4d706 100644
--- a/Changes
+++ b/Changes
@@ -1,5 +1,10 @@
Revision history for Perl extension Math::Prime::Util.
+0.07 17 June 2012
+ - Fixed a bug in next_prime found by Lou Godio (thank you VERY much!).
+ Added more tests for this. This had been changed in another area but
+ hadn't been brought into next_prime.
+
0.06 14 June 2012
- Change to New/Safefree from malloc. Oops.
diff --git a/MANIFEST b/MANIFEST
index 7102ba0..6432974 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -16,6 +16,7 @@ sieve.h
sieve.c
util.h
util.c
+examples/README
examples/bench-factor.pl
examples/bench-factor-extra.pl
examples/bench-factor-semiprime.pl
@@ -23,8 +24,12 @@ examples/bench-is-prime.pl
examples/bench-miller-rabin.pl
examples/bench-nthprime.pl
examples/bench-pcapprox.pl
+examples/bench-primecount.pl
examples/bench-random-prime.pl
-examples/test-factoring.pl
+examples/test-factor-yafu.pl
+examples/test-factor-mpxs.pl
+examples/test-nextprime-yafu.pl
+examples/test-primes-yafu.pl
examples/test-holf.pl
examples/test-nthapprox.pl
examples/test-pcapprox.pl
diff --git a/README b/README
index ed662a0..7d7b7f3 100644
--- a/README
+++ b/README
@@ -1,4 +1,4 @@
-Math::Prime::Util version 0.06
+Math::Prime::Util version 0.07
A set of utilities related to prime numbers. These include multiple sieving
methods, is_prime, prime_count, nth_prime, approximations and bounds for
diff --git a/TODO b/TODO
index 93c3a6a..303f8c6 100644
--- a/TODO
+++ b/TODO
@@ -37,3 +37,5 @@
- Iterator or tie?
- Get rid of erat_simple and bitarray.h. They're only there for comparison.
+
+- Need to be more aware of precision for the math operations.
diff --git a/XS.xs b/XS.xs
index 4053120..3f460f4 100644
--- a/XS.xs
+++ b/XS.xs
@@ -283,6 +283,7 @@ factor(IN UV n)
if (!split_success) {
split_success = holf_factor(n, factor_stack+nstack, 2000)-1;
}
+ /* Very, very few numbers make it past here */
}
if (split_success) {
MPUassert( split_success == 1, "split factor returned more than 2 factors");
diff --git a/examples/README b/examples/README
new file mode 100644
index 0000000..ac51401
--- /dev/null
+++ b/examples/README
@@ -0,0 +1,108 @@
+
+There are two main types of scripts here: benchmarks and correctness tests.
+
+The test-* scripts are generally trying to test one part of the module
+against another part of the module, an external module, or an external program.
+These usually consist of a combination of fixed tests and a long sequence of
+testing with random numbers, trying to find things the standard testing might
+have missed.
+
+
+test-factor-yafu.pl
+
+ Tests factorization compared with YAFU (v1.31.1). No arguments.
+
+test-factor-mpxs.pl
+
+ Tests factorization compared with Math::Factor::XS (v0.26).
+ One argument gives the number of random tests to perform.
+
+test-nextprime-yafu.pl
+
+ Tests next_prime() compared with YAFU (v1.31.1). No arguments.
+
+test-primes-yafu.pl
+
+ Tests primes($a,$b+$interval) compared with YAFU (v1.31.1). No arguments.
+ The interval is currently 8000.
+
+test-holf.pl
+
+ Tests the holf_factor() function vs. the factor() function. Given enough
+ rounds, HOLF (like Fermat) should be able to factor a number. We keep
+ calling it on each non-prime return value until it's done.
+
+test-nthapprox.pl
+
+ Tests the nth_prime approximation and upper/lower bounds vs. known values
+ for the nth prime on large values.
+
+test-pcapprox.pl
+
+ Tests the prime_count approximation and upper/lower bounds vs. known values
+ for Pi(x) on large values. Also examines the Schoenfeld and Stoll
+ inequalities.
+
+bench-factor-extra.pl
+
+ Benchmarks the various factoring options (prho, pbrent, pminus1, fermat, holf,
+ squfof, trial) on random n-digit numbers. Also gives the percent of the time
+ solutions were found with the given number of rounds (256k for SQUFOF, 2000
+ for HOLF, and 400 for probabilistic methods).
+
+bench-factor.pl
+
+ Benchmarks factoring random and semiprime n-digit numbers using the factor
+ method, and compares vs. Math::Factor::XS (v0.26). The latter uses a trial
+ division algorithm. MPU 0.05 and later use a threshold of 10M (8 digits)
+ to switch between trial division and methods like SQUFOF, Pollard's Rho, and
+ HOLF.
+
+bench-factor-semiprime.pl
+
+ Benchmarks factoring semiprimes, and compares with Math::Prime::XS and
+ Math::Pari. Takes two optional arguments: the number of digits (default 15)
+ and the benchmark count (default -2, meaning 2 seconds).
+
+bench-is-prime.pl
+
+ Benchmarks is_prime on random n-digit numbers, n from 5 to 10/20. Compares
+ MPU's is_prime and is_prob_prime vs. Math::Prime::XS and optionally
+ Math::Primality, Math::Pari, and Math::Prime::FastSieve. The first two of
+ the optional modules use methods more appropriate for big numbers, so are up
+ to an order of magnitude slower for 64-bit numbers. The last module (MPFS)
+ is extremely fast, but requires presieving to at least the number to be
+ tested, which is great for small numbers, but not for large.
+ Also, no additional precalc is done for MPU. If you really want blazing
+ fast is_prime and don't care about the memory and time to sieve, run
+ prime_precalc to the limit of your numbers and is_prime will turn into a
+ sub-100 microsecond bit array lookup for any number in the range.
+ Takes one optional argument of the benchmark count (default -5).
+
+bench-miller-rabin.pl
+
+ Benchmarks the strong miller_rabin test using 7 bases at various digit counts.
+ Takes one optional argument of the benchmark count (default -5).
+
+bench-nthprime.pl
+
+ Benchmarks the speed of nth_prime with various digit sizes.
+
+bench-pcapprox.pl
+
+ Benchmarks the speed of prime_count related functions for random n-digit
+ numbers (n = 5 to 10/20). This includes lower bound, lower/upper average,
+ the prime_count_approx function, li(x), and R(x).
+ Takes one optional argument of the benchmark count (default -5).
+
+bench-primecount.pl
+
+ Benchmarks the speed of prime_count on random n-digit numbers (n = 2 .. 8).
+ Takes one optional argument of the benchmark count (default -5).
+ I'll note you can easily see the transition from where we're just counting
+ existing value to where we have to sieve. Adding a prime_precalc(10**9)
+ line will speed up the 5-,6-,7-, and 8-digit prime_counts greatly.
+
+bench-random-prime.pl
+
+ Benchmarks the speed of random_ndigit_prime for various digits.
diff --git a/examples/bench-factor-extra.pl b/examples/bench-factor-extra.pl
index 18962c4..160f759 100755
--- a/examples/bench-factor-extra.pl
+++ b/examples/bench-factor-extra.pl
@@ -11,7 +11,7 @@ my $maxdigits = ($is64bit) ? 20 : 10; # Noting the range is limited for max.
srand(29);
my $rounds = 400;
my $sqrounds = 256*1024;
-my $hrounds = 100000;
+my $hrounds = 2000;
test_at_digits($_) for ( 3 .. $maxdigits );
diff --git a/examples/bench-primecount.pl b/examples/bench-primecount.pl
new file mode 100755
index 0000000..32b2f0b
--- /dev/null
+++ b/examples/bench-primecount.pl
@@ -0,0 +1,38 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util ":all";
+use Benchmark qw/:all/;
+my $maxdigits = (~0 <= 4294967295) ? 10 : 20;
+my $nnums = 100;
+
+my $count = shift || -5;
+
+srand(29);
+my @darray;
+push @darray, [gendigits($_)] for (2 .. 10);
+
+ my $sum;
+ cmpthese($count,{
+ ' 2' => sub { $sum += prime_count($_) for @{$darray[2-2]} },
+ ' 3' => sub { $sum += prime_count($_) for @{$darray[3-2]} },
+ ' 4' => sub { $sum += prime_count($_) for @{$darray[4-2]} },
+ ' 5' => sub { $sum += prime_count($_) for @{$darray[5-2]} },
+ ' 6' => sub { $sum += prime_count($_) for @{$darray[6-2]} },
+ ' 7' => sub { $sum += prime_count($_) for @{$darray[7-2]} },
+ ' 8' => sub { $sum += prime_count($_) for @{$darray[8-2]} },
+ #' 9' => sub { $sum += prime_count($_) for @{$darray[9-2]} },
+ #'10' => sub { $sum += prime_count($_) for @{$darray[10-2]} },
+ });
+ print "\n";
+
+sub gendigits {
+ my $digits = shift;
+ die "Digits must be > 0" unless $digits > 0;
+
+ my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
+ my $max = int(10 ** $digits);
+ $max = ~0 if $max > ~0;
+ my @nums = map { $base+int(rand($max-$base)) } (1 .. $nnums);
+ return @nums;
+}
diff --git a/examples/test-factoring.pl b/examples/test-factor-mpxs.pl
similarity index 96%
rename from examples/test-factoring.pl
rename to examples/test-factor-mpxs.pl
index a559e60..da6e4ff 100755
--- a/examples/test-factoring.pl
+++ b/examples/test-factor-mpxs.pl
@@ -33,6 +33,6 @@ while ($nrandom-- > 0) {
my @mpu = sort { $a<=>$b } factor($n);
die "failure for $n" unless scalar @mfxs == scalar @mpu;
for (0 .. $#mfxs) { die "failure for $n" unless $mfxs[$_] == $mpu[$_]; }
- print "." if ($nrandom % 1024) == 0;
+ print "." if ($nrandom % 256) == 0;
}
print "\n";
diff --git a/examples/test-factor-yafu.pl b/examples/test-factor-yafu.pl
new file mode 100755
index 0000000..bec01f2
--- /dev/null
+++ b/examples/test-factor-yafu.pl
@@ -0,0 +1,105 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util qw/factor/;
+use File::Temp qw/tempfile/;
+use autodie;
+use Text::Diff;
+my $maxdigits = (~0 <= 4294967295) ? 10 : 20;
+$| = 1; # fast pipes
+my $num = 10000;
+my $yafu_fname = "yafu_batchfile_$$.txt";
+$SIG{'INT'} = \&gotsig;
+
+
+{ # Test from 2 to 10000
+ print " 2 - 1000"; test_array( 2 .. 1000);
+ print " 1001 - 5000"; test_array( 1001 .. 5000);
+ print " 5001 - 10000"; test_array( 5001 .. 10000);
+}
+
+foreach my $digits (5 .. $maxdigits) {
+ printf "%5d %2d-digit numbers", $num, $digits;
+ my @narray = gendigits($digits, $num);
+ test_array(@narray);
+ $num = int($num * 0.9); # reduce as we go
+}
+
+sub test_array {
+ my @narray = @_;
+ print ".";
+ my @mpuarray = mpu_factors(@narray);
+ print ".";
+ my @yafuarray = yafu_factors(@narray);
+ print ".";
+ if ($#mpuarray != $#yafuarray) {
+ die "MPU got $#mpuarray primes, YAFU got $#yafuarray\n";
+ }
+ foreach my $n (@narray) {
+ my @mpu = @{shift @mpuarray};
+ my @yafu = @{shift @yafuarray};
+ die "mpu array is for the wrong n?" unless $n == shift @mpu;
+ die "yafu array is for the wrong n?" unless $n == shift @yafu;
+ my $diff = diff \@mpu, \@yafu, { STYLE => 'Table' };
+ die "factor($n):\n$diff\n" if length($diff) > 0;
+ }
+ print ".";
+ print "OK\n";
+}
+
+sub gendigits {
+ my $digits = shift;
+ die "Digits must be > 0" unless $digits > 0;
+ my $howmany = shift;
+
+ my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
+ my $max = int(10 ** $digits);
+ $max = ~0 if $max > ~0;
+ my @nums = map { $base+int(rand($max-$base)) } (1 .. $howmany);
+ return @nums;
+}
+
+sub mpu_factors {
+ my @piarray;
+ push @piarray, [$_, sort {$a<=>$b} factor($_)] for @_;
+ @piarray;
+}
+
+sub yafu_factors {
+ my @ns = @_;
+ my @piarray;
+
+ #my $fh = File::Temp->new; # .... autodie
+ #print $fh, "$_\n" for @_;
+ #$fh->flush;
+
+ # Shudder. Yafu must have a file in the current directory.
+ open(my $fh, '>', $yafu_fname);
+ print $fh "$_\n" for @ns;
+ close $fh;
+
+ open my $yafu, "yafu \"factor(\@)\" -batchfile $yafu_fname |";
+ my @curfactors;
+ while (<$yafu>) {
+ chomp;
+ if (/^P(RP)?\d+ = (\d+)/) {
+ push @curfactors, $2;
+ } elsif (/^C\d+ = (\d+)/) {
+ # Yafu didn't factor this one completely. Sneakily do it ourselves.
+ push @curfactors, factor($1);
+ } elsif (/ans = (\d+)/) {
+ push @piarray, [shift @ns, sort {$a<=>$b} @curfactors];
+ @curfactors = ();
+ }
+ }
+ close($yafu);
+ @piarray;
+}
+sub gotsig { my $sig = shift; die "Die because SIG$sig\n"; }
+END {
+ unlink $yafu_fname if -e $yafu_fname;
+ # YAFU leaves stuff around
+ unlink "__tmpbatchfile" if -e "__tmpbatchfile";
+ unlink "session.log" if -e "session.log";
+ unlink "factor.log" if -e "factor.log";
+}
diff --git a/examples/test-nextprime-yafu.pl b/examples/test-nextprime-yafu.pl
new file mode 100755
index 0000000..e30c348
--- /dev/null
+++ b/examples/test-nextprime-yafu.pl
@@ -0,0 +1,82 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util qw/next_prime/;
+use File::Temp qw/tempfile/;
+use autodie;
+my $maxdigits = (~0 <= 4294967295) ? 10 : 20;
+$| = 1; # fast pipes
+my $num = 10000;
+my $yafu_fname = "yafu_batchfile_$$.txt";
+$SIG{'INT'} = \&gotsig;
+
+foreach my $digits (1 .. $maxdigits) {
+ printf "%2d-digit numbers", $digits;
+ my @narray = gendigits($digits, $num);
+ print ".";
+ my @mpuarray = mpu_next_primes(@narray);
+ print ".";
+ die "mpu_next_primes didn't get enough numbers" unless $#mpuarray == $#narray;
+ my @yafuarray = yafu_next_primes(@narray);
+ die "yafunext_primes didn't get enough numbers" unless $#yafuarray == $#narray;
+ print ".";
+ foreach my $n (@narray) {
+ my $mpu = shift @mpuarray;
+ my $yafu = shift @yafuarray;
+ die "next_prime($n): MPU: $mpu YAFU: $yafu\n" unless $mpu == $yafu;
+ }
+ print ".";
+ print "OK\n";
+}
+
+sub gendigits {
+ my $digits = shift;
+ die "Digits must be > 0" unless $digits > 0;
+ my $howmany = shift;
+
+ my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
+ my $max = int(10 ** $digits);
+ $max = ~0 if $max > ~0;
+ my @nums = map { $base+int(rand($max-$base)) } (1 .. $howmany);
+ return @nums;
+}
+
+sub mpu_next_primes {
+ my @nparray;
+ push @nparray, next_prime($_) for @_;
+ @nparray;
+}
+
+sub yafu_next_primes {
+ my @nparray;
+ # Yafu 1.31 seems to go out of its way to make it hard to process more than
+ # one number at a time. The batchfile system will infinite loop if the data
+ # file isn't in the current directory.
+ # It does its darndest to see if you're on a terminal or not, and if not it
+ # just cuts you off after one number. So any sort of tempfile or pipe stuff
+ # just plain doesn't work. Faking it using IO::*tty* would probably work.
+
+ #my $fh = File::Temp->new; # .... autodie
+ #print $fh, "$_\n" for @_;
+ #$fh->flush;
+
+ # Shudder. Read comments above about why I have to do this.
+ open(my $fh, '>', $yafu_fname);
+ print $fh "$_\n" for @_;
+ close $fh;
+
+ open my $yafu, "yafu \"nextprime(\@)\" -batchfile $yafu_fname |";
+ while (<$yafu>) {
+ next unless /ans = (\d+)/;
+ push @nparray, $1;
+ }
+ close($yafu);
+ @nparray;
+}
+sub gotsig { my $sig = shift; die "Die because SIG$sig\n"; }
+END {
+ unlink $yafu_fname if -e $yafu_fname;
+ # YAFU leaves stuff around
+ unlink "__tmpbatchfile" if -e "__tmpbatchfile";
+ unlink "session.log" if -e "session.log";
+}
diff --git a/examples/test-pcapprox.pl b/examples/test-pcapprox.pl
index 70b1041..cae95d2 100755
--- a/examples/test-pcapprox.pl
+++ b/examples/test-pcapprox.pl
@@ -55,7 +55,7 @@ foreach my $n (sort {$a<=>$b} keys %pivals) {
my $pin = $pivals{$n};
my $pcl = prime_count_lower($n);
my $pcu = prime_count_upper($n);
- my ($scl,$scu) = schoenfeld($n);
+ my ($scl,$scu) = stoll($n);
#printf "10^%2d %12d %12d\n", length($n)-1, $pin-$pcl, $pcu-$pin;
printf "10^%2d %12.7f %12.7f %12.7f %12.7f\n",
@@ -66,7 +66,15 @@ foreach my $n (sort {$a<=>$b} keys %pivals) {
sub schoenfeld {
my $x = shift;
my $lix = LogarithmicIntegral($x);
- my $bound = ((sqrt($x)*log($x)) / 8*3.1415926535);
+ my $bound = (sqrt($x)*log($x)) / 8*3.1415926535;
+ ($lix-$bound,$lix+$bound);
+}
+
+# http://www.ams.org/journals/mcom/2011-80-276/S0025-5718-2011-02477-4/home.html
+sub stoll {
+ my $x = shift;
+ my $lix = LogarithmicIntegral($x);
+ my $bound = sqrt($x) * (log(log(log($x))) + exp(1) + 1) / (exp(1)*log($x));
($lix-$bound,$lix+$bound);
}
diff --git a/examples/test-primes-yafu.pl b/examples/test-primes-yafu.pl
new file mode 100755
index 0000000..e22b0e9
--- /dev/null
+++ b/examples/test-primes-yafu.pl
@@ -0,0 +1,107 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util qw/primes/;
+use File::Temp qw/tempfile/;
+use autodie;
+use Text::Diff;
+my $maxdigits = (~0 <= 4294967295) ? 10 : 18;
+$| = 1; # fast pipes
+my $num = 5000;
+my $interval = 8000;
+my $yafu_fname = "yafu_batchfile_$$.txt";
+$SIG{'INT'} = \&gotsig;
+
+
+# Note -- yafu 1.31 will not sieve 19 digit numbers. E.g.:
+# primes(8631424695497106432,8631424695497114432,0)
+# gives:
+# input too high
+
+
+foreach my $digits (3 .. $maxdigits) {
+ printf "%5d %2d-digit numbers", $num, $digits;
+ my @narray = gendigits($digits, $num);
+ print ".";
+ my @mpuarray = mpu_primes(@narray);
+ print ".";
+ #die "mpu_next_primes didn't get enough numbers" unless $#mpuarray-1 == $#narray;
+ my @yafuarray = yafu_primes(@narray);
+ #die "yafunext_primes didn't get enough numbers" unless $#yafuarray-1 == $#narray;
+ print ".";
+ if ($#mpuarray != $#yafuarray) {
+ die "MPU got $#mpuarray primes, YAFU got $#yafuarray\n";
+ }
+ foreach my $n (@narray) {
+ my @mpu = @{shift @mpuarray};
+ my @yafu = @{shift @yafuarray};
+ die "mpu array is for the wrong n?" unless $n == shift @mpu;
+ die "yafu array is for the wrong n?" unless $n == shift @yafu;
+ my $diff = diff \@mpu, \@yafu, { STYLE => 'Table' };
+ die "primes($n,$n+$interval):\n$diff\n" if length($diff) > 0;
+ }
+ print ".";
+ print "OK\n";
+ $num = int($num * 0.75); # reduce as we go
+}
+
+sub gendigits {
+ my $digits = shift;
+ die "Digits must be > 0" unless $digits > 0;
+ my $howmany = shift;
+
+ my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
+ my $max = int(10 ** $digits);
+ #$max = ~0 if $max > ~0;
+ $max = ~0-$interval if $max > (~0-$interval); # special for us
+ my @nums = map { $base+int(rand($max-$base)) } (1 .. $howmany);
+ return @nums;
+}
+
+sub mpu_primes {
+ my @piarray;
+ push @piarray, [$_, @{primes($_, $_+$interval)}] for @_;
+ @piarray;
+}
+
+sub yafu_primes {
+ my @ns = @_;
+ my @piarray;
+ # Yafu 1.31 seems to go out of its way to make it hard to process more than
+ # one number at a time. The batchfile system will infinite loop if the data
+ # file isn't in the current directory.
+ # It does its darndest to see if you're on a terminal or not, and if not it
+ # just cuts you off after one number. So any sort of tempfile or pipe stuff
+ # just plain doesn't work. Faking it using IO::*tty* would probably work.
+
+ #my $fh = File::Temp->new; # .... autodie
+ #print $fh, "$_\n" for @_;
+ #$fh->flush;
+
+ # Shudder. Read comments above about why I have to do this.
+ open(my $fh, '>', $yafu_fname);
+ print $fh "$_,", $_+$interval, ",0\n" for @ns;
+ close $fh;
+
+ open my $yafu, "yafu \"primes(\@)\" -pscreen -batchfile $yafu_fname |";
+ my @curprimes;
+ while (<$yafu>) {
+ chomp;
+ if (/^\d+/) {
+ push @curprimes, split(/ /);
+ } elsif (/ans = (\d+)/) {
+ foreach my $p (@curprimes) { die "Entry is '$p'" unless $p =~ /^\d+$/; }
+ push @piarray, [shift @ns, @curprimes];
+ @curprimes = ();
+ }
+ }
+ close($yafu);
+ @piarray;
+}
+sub gotsig { my $sig = shift; die "Die because SIG$sig\n"; }
+END {
+ unlink $yafu_fname if -e $yafu_fname;
+ # YAFU leaves stuff around
+ unlink "__tmpbatchfile" if -e "__tmpbatchfile";
+ unlink "session.log" if -e "session.log";
+}
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index cc22f8e..e225335 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -5,7 +5,7 @@ use Carp qw/croak confess/;
BEGIN {
$Math::Prime::Util::AUTHORITY = 'cpan:DANAJ';
- $Math::Prime::Util::VERSION = '0.06';
+ $Math::Prime::Util::VERSION = '0.07';
}
# parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
@@ -213,7 +213,7 @@ Math::Prime::Util - Utilities related to prime numbers, including fast sieves an
=head1 VERSION
-Version 0.06
+Version 0.07
=head1 SYNOPSIS
diff --git a/t/12-nextprime.t b/t/12-nextprime.t
index 6cd225b..dd8a4eb 100644
--- a/t/12-nextprime.t
+++ b/t/12-nextprime.t
@@ -8,7 +8,7 @@ use Math::Prime::Util qw/next_prime prev_prime/;
my $use64 = Math::Prime::Util::_maxbits > 32;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
-plan tests => 499*2 + 3*2 + 6 + 2;
+plan tests => 3573*2 + 3*2 + 6 + 2 + 148 + 148 + 1;
my @small_primes = qw/
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71
@@ -38,11 +38,19 @@ my @small_primes = qw/
3433 3449 3457 3461 3463 3467 3469 3491 3499 3511 3517 3527 3529 3533 3539 3541 3547 3557 3559 3571
/;
-for (my $i = 0; $i < (scalar @small_primes) - 1; $i++) {
- my $n = next_prime($small_primes[$i]);
- is($n, $small_primes[$i+1], "the next prime after $small_primes[$i] is $small_primes[$i+1] ?= $n");
- my $p = prev_prime($small_primes[$i+1]);
- is($p, $small_primes[$i], "the prev prime before $small_primes[$i+1] is $small_primes[$i] ?= $n");
+{
+ # insert primes before and after
+ unshift @small_primes, 0;
+ push @small_primes, 3581;
+ # Now test next_prime and prev_prime for all numbers 0 to 3572
+ my $prev_index = 0;
+ my $next_index = 1;
+ foreach my $n (0 .. 3572) {
+ $next_index++ if $n >= $small_primes[$next_index];
+ $prev_index++ if $n > $small_primes[$prev_index+1];
+ is(next_prime($n), $small_primes[$next_index], "next_prime($n) == $small_primes[$next_index]");
+ is(prev_prime($n), $small_primes[$prev_index], "prev_prime($n) == $small_primes[$prev_index]");
+ }
}
my %primegaps = (
@@ -65,3 +73,13 @@ 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" );
+
+# Turns out the testing of prev/next from 0-3572 still misses some cases.
+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)");
diff --git a/util.c b/util.c
index 614136f..86983ca 100644
--- a/util.c
+++ b/util.c
@@ -154,7 +154,8 @@ UV next_trial_prime(UV n)
d = n/30;
m = n - d*30;
- m = nextwheel30[m]; if (m == 1) d++;
+ /* Move forward one, knowing we may not be on the wheel */
+ if (m == 29) { d++; m = 1; } else { m = nextwheel30[m]; }
while (!_is_trial_prime7(d*30+m)) {
m = nextwheel30[m]; if (m == 1) d++;
}
@@ -189,8 +190,10 @@ UV next_prime(UV n)
d = n/30;
m = n - d*30;
- m = nextwheel30[m]; if (m == 1) d++;
+ /* Move forward one, knowing we may not be on the wheel */
+ if (m == 29) { d++; m = 1; } else { m = nextwheel30[m]; }
while (!_is_prime7(d*30+m)) {
+ /* Move forward one, knowing we are on the wheel */
m = nextwheel30[m]; if (m == 1) d++;
}
return(d*30+m);
--
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