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 <[email protected]> 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 _______________________________________________ Pkg-perl-cvs-commits mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
