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 badbbc5ff71dd6af55544efec84b3130ff1801c7 Author: Dana Jacobsen <[email protected]> Date: Sat Oct 20 03:37:51 2012 -0600 Add Lucky primes, make Cuban primes via a generator instead of filter (hugely faster) --- bin/primes.pl | 150 ++++++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 121 insertions(+), 29 deletions(-) diff --git a/bin/primes.pl b/bin/primes.pl index 5d7a427..a0c11c7 100755 --- a/bin/primes.pl +++ b/bin/primes.pl @@ -11,6 +11,7 @@ $| = 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 @@ -20,6 +21,7 @@ $| = 1; # 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: @@ -38,7 +40,35 @@ $| = 1; # 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 + +# There are a few of these prime filters that Math::NumSeq supports, and in +# theory it will add them eventually since they are OEIS sequences. Many are +# of the form "primes from ####" so aren't hard to work up. Math::NumSeq is +# a really neat module for playing with OEIS sequences. +# +# Example: All Sophie Germain primes under 1M +# primes.pl --sophie 1 1000000 +# perl -MMath::NumSeq::SophieGermainPrimes=:all -E 'my $seq = Math::NumSeq::SophieGermainPrimes->new; my $v = 0; while (1) { $v = ($seq->next)[1]; last if $v > $end; say $v; } BEGIN {our $end = 1000000}' +# +# Timing from 1 .. N for small N is going to be similar. As N increases, the +# time difference grows rapidly. +# +# primes.pl Math::NumSeq::SophieGermainPrimes +# 1M 0.14s 0.18s +# 10M 0.82s 3.89s +# 100M 7.09s 793s +# 1000M 64.2s ? estimated >3 days +# +# If given a non-zero start value it spreads even more, as for most sequences +# primes.pl doesn't have to generate preceeding values, while NumSeq has to +# start at the beginning. Additionally, Math::NumSeq may or may not deal with +# numbers larger than 2^32 (many sequences do, but it uses Math::Factor::XS +# for factoring and primality, which is limited to 32-bit). +# +# Here's an example of a combination. Palindromic primes: +# primes.pl --palin 1 1000000000 +# perl -MMath::Prime::Util=is_prime -MMath::NumSeq::Palindromes=:all -E 'my $seq = Math::NumSeq::Palindromes->new; my $v = 0; while (1) { $v = ($seq->next)[1]; last if $v > $end; say $v if is_prime($v); } BEGIN {our $end = 1000000000}' + my %opts; GetOptions(\%opts, 'safe|A005385', @@ -46,6 +76,7 @@ GetOptions(\%opts, 'twin|A001359', 'lucas|A005479', 'fibonacci|A005478', + 'lucky|A031157', 'triplet|A007529', 'quadruplet|A007530', 'cousin|A023200', @@ -74,12 +105,15 @@ $end =~ s/^(\d+)\*\*(\d+)$/Math::BigInt->new($1)->bpow($2)/e; die "$start isn't a positive integer" if $start =~ tr/0123456789//c; die "$end isn't a positive integer" if $end =~ tr/0123456789//c; -# Turn start and end into bigints if they're very large -if ( ($start >= ~0 && $start ne ~0) || ($end >= ~0 && $end ne ~0) ) { - $start = Math::BigInt->new($start); - $end = Math::BigInt->new($end); +# Turn start and end into bigints if they're very large. +# Fun fact: Math::BigInt->new("1") <= 10000000000000000000 is false. Sigh. +if ( ($start >= 2**63) || ($end >= 2**63) ) { + $start = Math::BigInt->new("$start") unless ref($start) eq 'Math::BigInt'; + $end = Math::BigInt->new("$end") unless ref($end) eq 'Math::BigInt'; } +my $segment_size = $start - $start + 30 * 128_000; # 128kB + # Calculate the mod 210 pre-test. This helps with the individual filters, # but the real benefit is that it convolves the pretests, which can speed # up even more. @@ -91,10 +125,13 @@ if (scalar keys %mod_pass == 0) { if ($start > $end) { # Do nothing -} elsif (exists $opts{'lucas'} || - exists $opts{'fibonacci'} || - exists $opts{'euclid'} || - exists $opts{'mersenne'} +} elsif ( exists $opts{'lucas'} + || exists $opts{'fibonacci'} + || exists $opts{'euclid'} + || exists $opts{'lucky'} + || exists $opts{'mersenne'} + || exists $opts{'cuban1'} + || exists $opts{'cuban2'} ) { my @p = gen_and_filter($start, $end); print join("\n", @p), "\n" if scalar @p > 0; @@ -112,7 +149,7 @@ if ($start > $end) { } my $seg_start = $start; - my $seg_end = $start + $segment_size; + my $seg_end = int($start + $segment_size); $seg_end = $end if $end < $seg_end; $start = $seg_end+1; @@ -142,21 +179,21 @@ if ($start > $end) { } } -# Return all Lucas primes between start and end, using identity: -# L_n = F_n-1 + F_n+1 +# This is OEIS A000032, Lucas numbers beginning at 2. +# Use identity: L_n = F_n-1 + F_n+1. Would be faster if calculated directly. sub lucas_primes { my ($start, $end) = @_; my @lprimes; - my $prime = 2; my $k = 0; - while ($prime < $start) { + my $Lk = 2; + while ($Lk < $start) { $k++; - $prime = fib($k) + fib($k+2); + $Lk = fib($k+1) + fib($k-1); } - while ($prime <= $end) { - push @lprimes, $prime if is_prime($prime); + while ($Lk <= $end) { + push @lprimes, $Lk if is_prime($Lk); $k++; - $prime = fib($k) + fib($k+2); + $Lk = fib($k+1) + fib($k-1); } @lprimes; } @@ -164,11 +201,10 @@ sub lucas_primes { sub fibonacci_primes { my ($start, $end) = @_; my @fprimes; - my $Fk = 2; my $k = 3; + my $Fk = fib($k); while ($Fk < $start) { - $k++; - $Fk = fib($k); + $Fk = fib(++$k); } while ($Fk <= $end) { push @fprimes, $Fk if is_prime($Fk); @@ -207,6 +243,53 @@ sub euclid_primes { @eprimes; } +sub cuban_primes { + my ($start, $end, $add) = @_; + my @cprimes; + my $psub = ($add == 1) ? sub { 3*$_[0]*$_[0] + 3*$_[0] + 1 } + : sub { 3*$_[0]*$_[0] + 6*$_[0] + 4 }; + # Determine first y via quadratic equation (under-estimate) + my $y = ($start <= 2) ? 0 : + ($add == 1) + ? int((-3 + sqrt(3*3 - 4*3*(1-$start))) / (2*3)) + : int((-6 + sqrt(6*6 - 4*3*(4-$start))) / (2*3)); + die "Incorrect start calculation" if $y > 0 && $psub->($y - 1) >= $start; + + # skip forward until p >= start + $y++ while $psub->($y) < $start; + + my $p = $psub->($y); + while ($p <= $end) { + push @cprimes, $p if is_prime($p); + $p = $psub->(++$y); + } + @cprimes; +} + +sub lucky_primes { + my ($start, $end) = @_; + # First do a (very basic) lucky number sieve to generate A000959. + # Evens removed for k=1: + # my @lucky = map { $_*2+1 } (0 .. int(($end-1)/2)); + # Remove the 3rd elements for k=2: + # my @lucky = grep { my $m = $_ % 6; $m == 1 || $m == 3 } (0 .. $end); + # Remove the 4th elements for k=3: + my @lucky = grep { my $m = $_ % 21; $m != 18 && $m != 19 } + grep { my $m = $_ % 6; $m == 1 || $m == 3 } + map { $_*2+1 } (0 .. int(($end-1)/2)); + for (my $k = 3; $k < scalar @lucky; $k++) { + my $skip = $lucky[$k]; + my $index = $skip-1; + last if $index > $#lucky; + do { + splice(@lucky, $index, 1); + $index += $skip-1; + } while ($index <= $#lucky); + } + # Then restrict to primes to get A031157. + grep { $_ >= $start && is_prime($_) } @lucky; +} + # This is not a general palindromic digit function! sub ndig_palindromes { my $digits = shift; @@ -244,7 +327,7 @@ sub is_pillai { # 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) { + for (my $n = 2; $n < $p; $n++) { $n_factorial_mod_p->bmul($n)->bmod($p); next if $p % $n == 1; return 1 if $n_factorial_mod_p == ($p-1); @@ -297,6 +380,15 @@ sub gen_and_filter { if (exists $opts{'euclid'}) { merge_primes(\$gen, \@p, 'euclid', euclid_primes($start, $end, 1)); } + if (exists $opts{'lucky'}) { + merge_primes(\$gen, \@p, 'lucky', lucky_primes($start, $end)); + } + if (exists $opts{'cuban1'}) { + merge_primes(\$gen, \@p, 'cuban1', cuban_primes($start, $end, 1)); + } + if (exists $opts{'cuban2'}) { + merge_primes(\$gen, \@p, 'cuban2', cuban_primes($start, $end, 2)); + } if (exists $opts{'palindromic'}) { if (!defined $gen) { foreach my $d (length($start) .. length($end)) { @@ -360,13 +452,12 @@ sub gen_and_filter { if (exists $opts{'sophie'}) { @p = grep { is_prime( 2*$_+1 ); } @p; } - if (exists $opts{'cuban1'}) { - #@p = grep { my $n = sqrt((4*$_-1)/3); $n == int($n); } @p; - @p = grep { my $n = sqrt((4*$_-1)/3); 4*$_ == int($n)*int($n)*3+1; } @p; - } - if (exists $opts{'cuban2'}) { - @p = grep { my $n = sqrt(($_-1)/3); $_ == int($n)*int($n)*3+1; } @p; - } + #if (exists $opts{'cuban1'}) { + # @p = grep { my $n = sqrt((4*$_-1)/3); 4*$_ == int($n)*int($n)*3+1; } @p; + #} + #if (exists $opts{'cuban2'}) { + # @p = grep { my $n = sqrt(($_-1)/3); $_ == int($n)*int($n)*3+1; } @p; + #} if (exists $opts{'pnm1'}) { @p = grep { is_prime( primorial(Math::BigInt->new($_))-1 ) } @p; } @@ -438,6 +529,7 @@ to only those primes additionally meeting these conditions: --lucas Lucas L_p is prime --fibonacci Fibonacci F_p is prime --mersenne Mersenne M_p = 2^p-1 is prime + --lucky Lucky p is a lucky number --palindr Palindromic p is equal to p with its base-10 digits reversed --pillai Pillai n! % p = p-1 and p % n != 1 for some n --good Good p_n^2 > p_{n-i}*p_{n+i} for all i in (1..n-1) -- 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
