[libmath-prime-util-perl] 20/50: Next round of primes.pl mods
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:45:34 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.13
in repository libmath-prime-util-perl.
commit 6c850f8a310be40a48bf437d821b9904fc1c41d1
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Oct 18 22:56:42 2012 -0600
Next round of primes.pl mods
---
bin/primes.pl | 50 +++++++++++--
examples/make-script-test-data.pl | 153 ++++++++++++++++++++++++++++++++++++++
examples/test-primes-script.pl | 133 ++++++++++++++++++++++++---------
lib/Math/Prime/Util.pm | 2 +-
4 files changed, 293 insertions(+), 45 deletions(-)
diff --git a/bin/primes.pl b/bin/primes.pl
index f44420d..80bfb7d 100755
--- a/bin/primes.pl
+++ b/bin/primes.pl
@@ -11,6 +11,33 @@ $| = 1;
# http://en.wikipedia.org/wiki/List_of_prime_numbers
# http://mathworld.wolfram.com/IntegerSequencePrimes.html
+# This program shouldn't contain any special knowledge about the series
+# members other than perhaps the start. It can know patterns, but don't
+# include a static list of the members, for instance. It should actually
+# compute the entries in a range (though go ahead and be clever about it).
+# Example:
+# DO use knowledge that F_k is prime only if k <= 4 or k is prime.
+# DO use knowledge that safe primes are <= 7 or congruent to 11 mod 12.
+# DO NOT use knowledge that fibprime(14) = 19134702400093278081449423917
+
+# The various primorial primes are confusing. Some things to consider:
+# 1) there are two definitions of primorial: p# and p_n#
+# 2) three sequences:
+# p where 1+p# is prime
+# n where 1+p_n# is prime
+# p_n#+1 where 1+p_n# is prime
+# 3) intersections of sequences (e.g. p_n#+1 and p_n#-1)
+# 4) other sequences like A057705: p where p+1 is an A002110 primorial
+# plus all the crazy primorial sequences (unlikely to be confused)
+#
+# A005234 p where p#+1 prime
+# A136351 p# where p#+1 prime 2,6,30,210,2310,200560490130
+# A014545 n where p_n#+1 prime 1,2,3,4,5,11,75,171,172
+# A018239 p_n#+1 where p_n#+1 prime
+#
+# A006794 p where p#-1 prime 3,5,11,13,41,89,317,337
+# A057704 n where p_n#-1 prime 2,3,5,6,13,24,66,68,167
+
my $segment_size = 30 * 128_000; # 128kB
my %opts;
GetOptions(\%opts,
@@ -32,9 +59,11 @@ GetOptions(\%opts,
'pnp1|A005234',
'pnm1|A006794',
'euclid|A018239',
+ 'nompugmp', # turn off MPU::GMP for debugging
'help',
) || die_usage();
die_usage() if exists $opts{'help'};
+Math::Prime::Util::prime_set_config(gmp=>0) if exists $opts{'nompugmp'};
# Get the start and end values. Verify they're positive integers.
@@ -91,17 +120,17 @@ sub lucas_primes {
sub fibonacci_primes {
my ($start, $end) = @_;
my @fprimes;
- my $prime = 2;
+ my $Fk = 2;
my $k = 3;
- while ($prime < $start) {
+ while ($Fk < $start) {
$k++;
- $prime = fib($k);
+ $Fk = fib($k);
}
- while ($prime <= $end) {
- push @fprimes, $prime if is_prime($prime);
+ while ($Fk <= $end) {
+ push @fprimes, $Fk if is_prime($Fk);
# For all but k=4, F_k is prime only when k is prime.
$k = ($k <= 4) ? $k+1 : next_prime($k);
- $prime = fib($k);
+ $Fk = fib($k);
}
@fprimes;
}
@@ -163,6 +192,13 @@ sub is_pillai {
# foreach my $n (grep { ($p % $_) != 1 } (2 .. $p-1)) {
# return 1 if Math::BigInt->new($n)->bfac()->bmod($p) == ($p-1);
# }
+ # About 20% faster but sucks memory:
+ # foreach my $n (grep { ($p % $_) != 1 } (2 .. $p-1)) {
+ # $fac_c[$n] = Math::BigInt->new($n)->bfac() if !defined $fac_c[$n];
+ # return 1 if ($fac_c[$n] % $p) == ($p-1);
+ # }
+ # Once p gets large (say 20,000+), then calculating and storing n! is
+ # unreasonable, and the code below will be much faster.
my $n_factorial_mod_p = Math::BigInt->bone();
for my $n (2 .. $p-1) {
$n_factorial_mod_p->bmul($n)->bmod($p);
@@ -172,7 +208,7 @@ sub is_pillai {
0;
}
-# Also pretty slow.
+# Not nearly as slow as Pillai, but not fast.
sub is_good_prime {
my $p = shift;
return 0 if $p <= 2; # 2 isn't a good prime
diff --git a/examples/make-script-test-data.pl b/examples/make-script-test-data.pl
new file mode 100755
index 0000000..bfee664
--- /dev/null
+++ b/examples/make-script-test-data.pl
@@ -0,0 +1,153 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use File::Spec::Functions;
+use FindBin;
+use bigint;
+use Data::BitStream::XS;
+use Math::Prime::Util qw/is_prime/;
+$|++; #flush the output buffer after every write() or print() function
+
+
+my @test_data = (
+ [ 668, "Mersenne", "--mersenne", 10**100],
+ [ 7529, "Triplet", "--triplet", 0],
+ [ 7530, "Quadruplet", "--quadruplet", 0],
+ [23200, "Cousin", "--cousin", 0],
+ [23201, "Sexy", "--sexy", 0],
+ [ 1359, "Twin", "--twin", 0],
+ [ 5385, "Safe", "--safe", 0],
+ [ 5384, "SG", "--sophie", 0],
+ [ 2385, "Palindromic", "--palin", 32_965_656_923],
+ [ 5479, "Lucas", "--lucas", 0],
+ [ 5478, "Fibonacci", "--fibonacci", 0],
+ [63980, "Pillai", "--pillai", 2000],
+ [28388, "Good", "--good", 20000],
+ [ 2407, "Cuban y+1", "--cuban1", 0],
+ [ 2648, "Cuban y+2", "--cuban2", 100_000_000],
+ [ 5234, "Primorial+1", "--pnp1", 2500],
+ [ 6794, "Primorial-1", "--pnm1", 2500],
+ [18239, "Euclid", "--euclid", 0],
+);
+
+foreach my $test (@test_data) {
+ my $oeis_no = $test->[0];
+ my $filename = sprintf("b%06d.txt", $oeis_no);
+ my $link = sprintf("http://oeis.org/A%06d/b%06d.txt", $oeis_no, $oeis_no);
+ if (!-r $filename) {
+ warn "Getting $filename from $link\n";
+ qx/wget $link/;
+ die "Could not retrieve. Bailing\n" unless -r $filename;
+ }
+ my $ref_data = read_oeis(@$test);
+ push @$test, $ref_data;
+}
+
+my $stream = Data::BitStream::XS->new( file => 'script-test-data.bs', mode => 'w' );
+foreach my $test (@test_data) {
+ test_oeis(@$test);
+}
+$stream->write_close();
+
+sub read_oeis {
+ my($oeis_no, $name, $script_arg, $restrict) = @_;
+ die "Restrict isn't defined for $oeis_no : $name" unless defined $restrict;
+
+ my $filename = sprintf("b%06d.txt", $oeis_no);
+ my $link = sprintf("http://oeis.org/A%06d/b%06d.txt", $oeis_no, $oeis_no);
+
+ my @ref;
+ {
+ open my $fh, '<', $filename
+ or die "Can't read $filename.\nYou should run:\n wget $link\n";
+ printf "%12s primes: reading %12s...", $name, $filename;
+ my $char = " ";
+ while (<$fh>) {
+ next unless /^(\d+)\s+(\d+)/;
+ my $v = (length($2) < 20) ? $2 : Math::BigInt->new("$2");
+ if ($restrict > 0 && $v > $restrict) {
+ $char = '*';
+ last;
+ }
+ push @ref, $v;
+ }
+ close $fh;
+ print "$char";
+ }
+ printf " %7d.", scalar @ref;
+ print " Testing..";
+ if ($ref[-1] > 18446744073709551615) {
+ print ",";
+ # Check for monotonic and primeness
+ foreach my $i (0 .. $#ref) {
+ die "non-prime in $oeis_no $name\n" unless is_prime($ref[$i]);
+ if ($i > 0) {
+ die "non-monotonic sequence in $oeis_no $name ($i $ref[$i-1] $ref[$i])\n" if $ref[$i] <= $ref[$i-1];
+ die "even number in $oeis_no $name\n" if ($ref[$i] % 2) == 0;
+ }
+ }
+ } else {
+ no bigint;
+ print ".";
+ # Check for monotonic and primeness
+ foreach my $i (0 .. $#ref) {
+ die "non-prime in $oeis_no $name\n" unless is_prime($ref[$i]);
+ if ($i > 0) {
+ die "non-monotonic sequence in $oeis_no $name\n" if $ref[$i] <= $ref[$i-1];
+ die "even number in $oeis_no $name\n" if ($ref[$i] % 2) == 0;
+ }
+ }
+ }
+ print "done\n";
+ return \@ref;
+}
+
+sub test_oeis {
+ my($oeis_no, $name, $script_arg, $restrict, $ref_data) = @_;
+
+ my @ref = @$ref_data;
+ printf "%12s primes: stream..", $name;
+
+ if ($ref[-1] > 18446744073709551615) {
+ print ",";
+ # Store the first two values, then a list of deltas
+ $stream->put_gamma($oeis_no, 1, scalar @ref, $ref[0], $ref[1]);
+ print ".";
+ my @deltas = map { ($ref[$_] - $ref[$_-1] - 2)/2 } (2..$#ref);
+ print ".";
+ # Ugly... Check for anything really big;
+ my @giant;
+ foreach my $d (@deltas) {
+ if ($d >= 18446744073709551614) {
+ push @giant, $d;
+ $d = 18446744073709551614;
+ }
+ }
+ print ".";
+ my $k = 2;
+ $stream->put_arice($k, @deltas);
+ print ".";
+ # Store giant deltas raw
+ foreach my $d (@giant) {
+ if (ref($d) ne 'Math::BigInt') {
+ warn "big delta $d isn't a bigint.\n";
+ $d = Math::BigInt->new(0);
+ }
+ my $binstr = substr($d->as_bin, 2);
+ $stream->put_gamma(length($binstr));
+ $stream->put_string($binstr);
+ }
+ } else {
+ no bigint;
+ print ".";
+ # Store the first two values, then a list of deltas
+ $stream->put_gamma($oeis_no, 0, scalar @ref, $ref[0], $ref[1]);
+ print ".";
+ my @deltas = map { ($ref[$_] - $ref[$_-1] - 2)/2 } (2..$#ref);
+ print ".";
+ my $k = 2;
+ $stream->put_arice($k, @deltas);
+ }
+
+ print "done\n";
+}
diff --git a/examples/test-primes-script.pl b/examples/test-primes-script.pl
index 12d2090..e91ccdb 100755
--- a/examples/test-primes-script.pl
+++ b/examples/test-primes-script.pl
@@ -3,51 +3,110 @@ use strict;
use warnings;
use File::Spec::Functions;
use FindBin;
+use Time::HiRes qw(gettimeofday tv_interval);
+use bigint;
+use Data::BitStream::XS;
$|++; #flush the output buffer after every write() or print() function
-test_oeis(668, "Mersenne", "--mersenne", '1' . '0' x 100);
-#test_oeis(2407, "Cuban y+1", "--cuban1");
-#test_oeis(2648, "Cuban y+2", "--cuban2", 100_000_000);
-test_oeis(5234, "Primorial+1", "--pnp1", 2500);
-test_oeis(6794, "Primorial-1", "--pnm1", 2500);
-test_oeis(18239, "Euclid", "--euclid");
-test_oeis(7529, "Triplet", "--triplet");
-test_oeis(7530, "Quadruplet", "--quadruplet");
-test_oeis(23200, "Cousin", "--cousin");
-test_oeis(23201, "Sexy", "--sexy");
-test_oeis(1359, "Twin", "--twin");
-test_oeis(5385, "Safe", "--safe");
-test_oeis(5384, "SG", "--sophie");
-test_oeis(2385, "Palindromic", "--palin");
-test_oeis(5479, "Lucas", "--lucas");
-test_oeis(5478, "Fibonacci", "--fibonacci");
-test_oeis(63980, "Pillai", "--pillai", 2000);
-test_oeis(28388, "Good", "--good", 20000);
+# Should use a "put_test_string" sort of call on the test data, so we
+# wouldn't need this.
+my @test_data = (
+ [ 668, "Mersenne", "--mersenne", 10**100],
+ [ 7529, "Triplet", "--triplet", 0],
+ [ 7530, "Quadruplet", "--quadruplet", 0],
+ [23200, "Cousin", "--cousin", 0],
+ [23201, "Sexy", "--sexy", 0],
+ [ 1359, "Twin", "--twin", 0],
+ [ 5385, "Safe", "--safe", 0],
+ [ 5384, "SG", "--sophie", 0],
+ [ 2385, "Palindromic", "--palin", 32_965_656_923],
+ [ 5479, "Lucas", "--lucas", 0],
+ [ 5478, "Fibonacci", "--fibonacci", 0],
+ [63980, "Pillai", "--pillai", 2000],
+ [28388, "Good", "--good", 20000],
+ [ 2407, "Cuban y+1", "--cuban1", 0],
+ [ 2648, "Cuban y+2", "--cuban2", 100_000_000],
+ [ 5234, "Primorial+1", "--pnp1", 2500],
+ [ 6794, "Primorial-1", "--pnm1", 2500],
+ [18239, "Euclid", "--euclid", 0],
+);
+my %oeis_name = map { $_->[0] => $_->[1] } @test_data;
-sub test_oeis {
- my($oeis_no, $name, $script_arg, $restrict) = @_;
+my $test_data_hash = read_script_data('script-test-data.bs');
- my $filename = sprintf("b%06d.txt", $oeis_no);
- my $link = sprintf("http://oeis.org/A%06d/b%06d.txt", $oeis_no, $oeis_no);
+foreach my $test (@test_data) {
+ my $oeis_no = $test->[0];
+ my $test_data = $test_data_hash->{$oeis_no};
+ if (!defined $test_data) {
+ die "No test data found for OEIS $oeis_no : $test->[1] primes\n";
+ }
+ test_oeis(@$test, $test_data);
+}
- my @ref;
- my @scr;
- {
- open my $fh, '<', $filename
- or die "Can't read $filename.\nYou should run:\n wget $link\n";
- printf "%12s primes: reading %15s...", $name, $filename;
- while (<$fh>) {
- next unless /^(\d+)\s+(\d+)/;
- last if defined $restrict && $2 > $restrict;
- push @ref, "$2";
+
+
+
+sub read_script_data {
+ my ($filename) = @_;
+
+ my $stream = Data::BitStream::XS->new( file => $filename, mode => 'ro' );
+ my %hash;
+
+ while (!$stream->exhausted) {
+ my ($oeis_no, $is_bigint, $num_entries, @ref) = $stream->get_gamma(5);
+ printf "%12s primes (OEIS A%06d): reading %7d entries..", $oeis_name{$oeis_no}, $oeis_no, $num_entries;
+ if ($is_bigint) {
+ print ",";
+ my $k = 2;
+ my @deltas = $stream->get_arice($k, $num_entries-2);
+ print ".";
+ # Check to see if we have any giant deltas
+ foreach my $d (@deltas) {
+ if ($d >= 18446744073709551614) {
+ my $len = $stream->get_gamma;
+ my $binstr = $stream->read_string($len);
+ $d = Math::BigInt->new('0b' . $binstr);
+ }
+ }
+ print ".";
+ my $prev = $ref[1];
+ push @ref, map { $prev = $_*2+$prev+2; } @deltas;
+ $hash{$oeis_no} = \@ref;
+ print ".\n";
+ } else {
+ no bigint;
+ print ".";
+ my $k = 2;
+ my @deltas = $stream->get_arice($k, $num_entries-2);
+ print ".";
+ my $prev = $ref[1];
+ push @ref, map { $prev = $_*2+$prev+2; } @deltas;
+ $hash{$oeis_no} = \@ref;
+ print ".\n";
}
- close $fh;
}
- printf " %7d.", scalar @ref;
+ \%hash;
+}
+
+sub test_oeis {
+ my($oeis_no, $name, $script_arg, $restrict, $ref_data) = @_;
- printf " Generating...";
- @scr = split /\s+/, qx+$FindBin::Bin/../bin/primes.pl $script_arg 1 $ref[-1]+;
- printf " %7d.\n", scalar @scr;
+ my @ref = @$ref_data;
+ printf "%12s primes (OEIS A%06d): generating..", $name, $oeis_no;
+
+ #print "\n";
+ #print "reference data:\n";
+ #print " $_\n" for @ref;
+ #print "primes.pl $script_arg 1 $ref[-1]\n";
+ my $start = [gettimeofday];
+ my @scr = split /\s+/, qx+$FindBin::Bin/../bin/primes.pl $script_arg 1 $ref[-1]+;
+ {
+ no bigint;
+ my $num_generated = scalar @scr;
+ my $seconds = tv_interval($start);
+ my $msperprime = ($seconds * 1000.0) / $num_generated;
+ printf " %7d. %7.2f ms/prime\n", $num_generated, $msperprime;
+ }
die "Not equal numbers: ", scalar @ref, " - ", scalar @scr, "\n"
unless @ref == @scr;
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 4b5a207..b858b03 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -663,7 +663,7 @@ sub primorial {
}
my $pn = 1;
- $pn = Math::BigInt->new->bone if defined $bigint::VERSION &&
+ $pn = Math::BigInt->new->bone if defined $Math::BigInt::VERSION &&
$n >= (($_Config{'maxbits'} == 32) ? 29 : 53);
foreach my $p ( @{ primes($n) } ) {
--
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