#!/usr/bin/env perl
use strict;
use warnings;
use Getopt::Long;
use Math::BigInt try => 'GMP';
use Math::Prime::Util qw/primes prime_count next_prime prev_prime
twin_primes sieve_prime_cluster mulmod is_pillai
is_prime is_provable_prime is_mersenne_prime
lucasu lucasv
nth_prime prime_count primorial pn_primorial/;
$| = 1;
# For many more types, see:
# 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
#
# As an aside, the 18th p#-1 is 15877, but the 19th is 843301.
# The p#+1's are a bit denser, with the 22nd at 392113.
# 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.06s 0.13s
# 10M 0.21 2.91
# 100M 1.52 396
# 1000M 13.7 > a day
#
# 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;
# Make Getopt not capture +
Getopt::Long::Configure(qw/no_getopt_compat/);
GetOptions(\%opts,
'safe|A005385',
'sophie|sg|A005384',
'twin|A001359',
'lucas|A005479',
'fibonacci|A005478',
'lucky|A031157',
'triplet|A007529',
'quadruplet|A007530',
'cousin|A023200',
'sexy|A023201',
'mersenne|A000668',
'palindromic|palindrome|palendrome|A002385',
'pillai|A063980',
'good|A028388',
'cuban1|A002407',
'cuban2|A002648',
'pnp1|A005234',
'pnm1|A006794',
'euclid|A018239',
'circular|A068652',
'panaitopol|A027862',
'provable',
'nompugmp', # turn off MPU::GMP for debugging
'version',
'help',
) || die_usage();
Math::Prime::Util::prime_set_config(gmp=>0) if exists $opts{'nompugmp'};
if (exists $opts{'version'}) {
my $version_str =
"primes.pl version 1.3 using Math::Prime::Util $Math::Prime::Util::VERSION";
$version_str .= " and MPU::GMP $Math::Prime::Util::GMP::VERSION"
if Math::Prime::Util::prime_get_config->{'gmp'};
$version_str .= "\nWritten by Dana Jacobsen.\n";
die "$version_str";
}
die_usage() if exists $opts{'help'};
# Get the start and end values. Verify they're positive integers.
@ARGV = (0,@ARGV) if @ARGV == 1;
die_usage() unless @ARGV == 2;
my ($start, $end) = @ARGV;
# Allow some expression evaluation on the input, but don't just eval it.
$end = "($start)$end" if $end =~ /^\+/;
$start =~ s/\s*$//; $start =~ s/^\s*//;
$end =~ s/\s*$//; $end =~ s/^\s*//;
$start = eval_expr($start) unless $start =~ /^\d+$/;
$end = eval_expr($end ) unless $end =~ /^\d+$/;
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.
# 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.
my ($min_pass, %mod_pass) = find_mod210_restriction();
# Find out if they've filtered so much nothing passes (e.g. cousin quad)
if (scalar keys %mod_pass == 0) {
$end = $min_pass if $end > $min_pass;
}
if ($start > $end) {
# Do nothing
} 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;
} else {
while ($start <= $end) {
# Adjust segment sizes for some cases
$segment_size = 10000 if $start > ~0; # small if doing bigints
if (exists $opts{'pillai'}) {
$segment_size = ($start < 10000) ? 100 : 1000; # very small for Pillai
}
if (exists $opts{'pnp1'} || exists $opts{'pnm1'}) {
$segment_size = 500;
}
if (exists $opts{'palindromic'}) {
$segment_size = 10**length($start) - $start - 1; # all n-digit numbers
}
if (exists $opts{'panaitopol'}) {
$segment_size = (~0 == 4294967295) ? 2147483648 : int(10**12);
}
my $seg_start = $start;
my $seg_end = int($start + $segment_size);
$seg_end = $end if $end < $seg_end;
$start = $seg_end+1;
my $p = gen_and_filter($seg_start, $seg_end);
# print this segment
print join("\n", @$p), "\n" if scalar @$p > 0;
}
}
# This is OEIS A000032, Lucas numbers beginning at 2.
sub lucas_primes {
my ($start, $end) = @_;
my ($k, $Lk, @lprimes) = (0);
do {
$Lk = lucasv(1,-1,$k);
push @lprimes, $Lk if $Lk >= $start && is_prime($Lk);
$k++;
} while $Lk < $end;
@lprimes;
}
sub fibonacci_primes {
my ($start, $end) = @_;
my ($k, $Fk, @fprimes) = (3);
do {
$Fk = lucasu(1,-1,$k);
push @fprimes, $Fk if $Fk >= $start && is_prime($Fk);
$k = ($k <= 4) ? $k+1 : next_prime($k);
} while $Fk < $end;
@fprimes;
}
sub mersenne_primes {
my ($start, $end) = @_;
my @mprimes;
my $p = 1;
while (1) {
$p = next_prime($p); # Mp is not prime if p is not prime
next if $p > 3 && ($p % 4) == 3 && is_prime(2*$p+1);
my $Mp = Math::BigInt->bone->blsft($p)->bdec;
last if $Mp > $end;
push @mprimes, $Mp if $Mp >= $start && is_mersenne_prime($p);
}
@mprimes;
}
sub euclid_primes {
my ($start, $end, $add) = @_;
my @eprimes;
my $k = 0;
while (1) {
my $primorial = pn_primorial(Math::BigInt->new($k)) + $add;
last if $primorial > $end;
push @eprimes, $primorial if $primorial >= $start && is_prime($primorial);
$k++;
}
@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 panaitopol_primes {
my ($start, $end) = @_;
my @init;
push @init, 5 if $start <= 5 && $end >= 5;
push @init, 13 if $start <= 13 && $end >= 13;
return @init if $end < 41;
my $nbeg = ($start <= 41) ? 4 : int( sqrt( ($start-1)/2) );
my $nend = int( sqrt(($end-1)/2) );
$nbeg++ while (2*$nbeg*($nbeg+1)+1) < $start;
$nend-- while (2*$nend*($nend+1)+1) > $end;
# TODO: BigInts
return @init,
grep { is_prime($_) }
grep { ($_%5) && ($_%13) && ($_%17) && ($_%29) && ($_%37) }
map { 2*$_*($_+1)+1 }
$nbeg .. $nend;
}
sub lucky_primes {
my ($start, $end) = @_;
# First do a lucky number sieve to generate A000959.
my @_lf63; # Lucky: 1,3,7,9,13,15,... 63=7*9.
$_lf63[$_] = 1 for (qw/2 5 8 11 14 17 18 19 20 23 26 27 28 29 32 35 38 39 40 41 44 47 50 53 56 57 58 59 60 61 62/);
my @lucky;
my $n = 1;
while ($n <= $end) {
my $m63 = $n % 63;
push @lucky, $n unless $_lf63[$m63];
push @lucky, $n+2 unless $_lf63[$m63+2];
$n += 6;
}
delete $lucky[-1] if $lucky[-1] > $end;
for (my $k = 4; $k < scalar @lucky && $lucky[$k]-1 <= $#lucky; $k++) {
my $skip = $lucky[$k]-1;
my $index = $skip;
while ($index <= $#lucky) {
splice(@lucky, $index, 1);
$index += $skip;
}
}
shift @lucky while $lucky[0] < $start;
# Then restrict to primes to get A031157.
grep { is_prime($_) } @lucky;
}
# This is not a general palindromic digit function!
sub ndig_palindromes {
my $digits = shift;
return (2,3,5,7) if $digits == 1;
return (11) if $digits == 2;
return () if ($digits % 2) == 0;
my $rhdig = int(($digits - 1) / 2);
return grep { is_prime($_) }
map { $_ . reverse substr($_,0,$rhdig) }
map { $_ * int(10**$rhdig) .. ($_+1) * int(10**$rhdig) - 1 }
1, 3, 7, 9;
}
# Not fast.
sub is_good_prime {
my $p = shift;
return 0 if $p <= 2; # 2 isn't a good prime
my $lower = $p;
my $upper = $p;
while ($lower > 2) {
$lower = prev_prime($lower);
$upper = next_prime($upper);
return 0 if ($p*$p) <= ($upper * $lower);
}
1;
}
# Assumes the input is prime. Returns 1 if all digit rotations are prime.
sub is_circular_prime {
my $p = shift;
return 1 if $p < 10;
return 0 if $p =~ tr/024568//;
# TODO: BigInts
foreach my $rot (1 .. length($p)-1) {
return 0 unless is_prime( substr($p, $rot) . substr($p, 0, $rot) );
}
1;
}
sub merge_primes {
my ($genref, $pref, $name, @primes) = @_;
if (!defined $$genref) {
@$pref = @primes;
$$genref = $name;
} else {
my %f;
undef @f{ @primes };
@$pref = grep { exists $f{$_} } @$pref;
}
}
# This is used for things that can generate a filtered list faster than
# searching through all primes in the range.
sub gen_and_filter {
my ($start, $end) = @_;
my $gen;
my $p = [];
$end-- if ($end % 2) == 0 && $end > 2;
if (exists $opts{'lucas'}) {
merge_primes(\$gen, $p, 'lucas', lucas_primes($start, $end));
}
if (exists $opts{'fibonacci'}) {
merge_primes(\$gen, $p, 'fibonacci', fibonacci_primes($start, $end));
}
if (exists $opts{'mersenne'}) {
merge_primes(\$gen, $p, 'mersenne', mersenne_primes($start, $end));
}
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{'panaitopol'}) {
merge_primes(\$gen, $p, 'panaitopol', panaitopol_primes($start, $end));
}
if (exists $opts{'palindromic'}) {
if (!defined $gen) {
foreach my $d (length($start) .. length($end)) {
push @$p, grep { $_ >= $start && $_ <= $end }
ndig_palindromes($d);
}
$gen = 'palindromic';
} else {
@$p = grep { $_ eq reverse $_; } @$p;
}
}
# Combine the cluster types and use an efficient cluster sieve if possible
if (!defined $gen) {
my @cluster;
if (defined $opts{'twin'}) { $cluster[2] = 1; }
if (defined $opts{'cousin'}) { $cluster[4] = 1; }
if (defined $opts{'sexy'}) { $cluster[6] = 1; }
if (defined $opts{'triplet'}) { $cluster[6] = 1; }
if (defined $opts{'quadruplet'}) { $cluster[$_] = 1 for (2,6,8); }
@cluster = grep { defined $cluster[$_] } 0 .. $#cluster;
if (scalar @cluster) {
if (scalar(@cluster) == 1 && $cluster[0] == 2) {
$p = twin_primes($start, $end);
} else {
$p = [sieve_prime_cluster($start, $end, @cluster)];
}
$gen = 'cluster';
}
}
if (!defined $gen) {
$p = primes($start, $end);
$gen = 'primes';
}
# Apply the mod 210 pretest
if ($min_pass > 0) {
@$p = grep { $_ <= $min_pass || exists $mod_pass{$_ % 210} } @$p;
}
# If we didn't generate the list with a cluster sieve, grep them out
if ($gen ne 'cluster') {
if (exists $opts{'twin'}) {
@$p = grep { is_prime( $_+2 ); } @$p;
}
if (exists $opts{'quadruplet'}) {
@$p = grep { is_prime($_+2) && is_prime($_+6) && is_prime($_+8); } @$p;
}
if (exists $opts{'triplet'}) {
@$p = grep { is_prime($_+6) && (is_prime($_+2) || is_prime($_+4)); } @$p;
}
if (exists $opts{'cousin'}) {
@$p = grep { is_prime($_+4); } @$p;
}
if (exists $opts{'sexy'}) {
@$p = grep { is_prime($_+6); } @$p;
}
} else {
# Cluster sieve for triplet gives us just p+6.
if (exists $opts{'triplet'} && !exists $opts{'twin'} && !exists $opts{'cousin'} && !exists $opts{'quadruplet'}) {
@$p = grep { is_prime($_+2) || is_prime($_+4); } @$p;
}
}
if (exists $opts{'safe'}) {
@$p = grep { is_prime( ($_-1) >> 1 ); }
grep { ($_ <= 7) || ($_ % 12) == 11; }
@$p;
}
if (exists $opts{'sophie'}) {
@$p = grep { is_prime( 2*$_+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;
}
if (exists $opts{'pnp1'}) {
@$p = grep { is_prime( primorial(Math::BigInt->new($_))+1 ) } @$p;
}
if (exists $opts{'circular'}) {
@$p = grep { is_circular_prime($_) } @$p;
}
if (exists $opts{'pillai'}) {
# See: http://en.wikipedia.org/wiki/Pillai_prime
@$p = grep { is_pillai($_); } @$p;
}
if (exists $opts{'good'}) {
@$p = grep { is_good_prime($_); } @$p;
}
if (exists $opts{'provable'}) {
@$p = grep { is_provable_prime($_) == 2; } @$p;
}
$p;
}
{
my %_mod210_restrict = (
cuban1 => {min=> 7, mod=>[1,19,37,61,79,121,127,169,187]},
cuban2 => {min=> 2, mod=>[1,13,43,109,139,151,169,181,193]},
twin => {min=> 5, mod=>[11,17,29,41,59,71,101,107,137,149,167,179,191,197,209]},
triplet => {min=> 7, mod=>[11,13,17,37,41,67,97,101,103,107,137,163,167,187,191,193]},
quadruplet => {min=> 5, mod=>[11,101,191]},
cousin => {min=> 7, mod=>[13,19,37,43,67,79,97,103,109,127,139,163,169,187,193]},
sexy => {min=> 7, mod=>[11,13,17,23,31,37,41,47,53,61,67,73,83,97,101,103,107,121,131,137,143,151,157,163,167,173,181,187,191,193]},
safe => {min=>11, mod=>[17,23,47,53,59,83,89,107,137,143,149,167,173,179,209]},
sophie => {min=> 5, mod=>[11,23,29,41,53,71,83,89,113,131,149,173,179,191,209]},
panaitopol => {min=> 5, mod=>[1,11,13,41,43,53,61,71,83,103,113,131,151,173,181,193]},
# Nothing for good, pillai, palindromic, fib, lucas, mersenne, primorials
);
sub find_mod210_restriction {
my %mods_left;
undef @mods_left{ grep { ($_%2) && ($_%3) && ($_%5) && ($_%7) } (0..209) };
my $min = 0;
while (my($filter,$data) = each %_mod210_restrict) {
next unless exists $opts{$filter};
$min = $data->{min} if $min < $data->{min};
my %thismod;
undef @thismod{ @{$data->{mod}} };
foreach my $m (keys %mods_left) {
delete $mods_left{$m} unless exists $thismod{$m};
}
}
return ($min, %mods_left);
}
}
# This is rather braindead. We're going to eval their input so they can give
# arbitrary expressions. But we only want to allow math-like strings.
sub eval_expr {
my $expr = shift;
die "$expr cannot be evaluated" if $expr =~ /:/; # Use : for escape
$expr =~ s/nth_prime\(/:1(/g;
$expr =~ s/log\(/:2(/g;
die "$expr cannot be evaluated" if $expr =~ tr|-0123456789+*/() :||c;
$expr =~ s/:1/nth_prime/g;
$expr =~ s/:2/log/g;
$expr =~ s/(\d+)/ Math::BigInt->new($1) /g;
my $res = eval $expr; ## no critic
die "Cannot eval: $expr\n" if !defined $res;
$res = int($res->bstr) if ref($res) eq 'Math::BigInt' && $res <= ~0;
$res;
}
sub die_usage {
die <<EOU;
Usage: $0 [options] [START] END
Displays all primes between the positive integers START and END, inclusive.
The START and END values must be integers or simple expressions. This allows
inputs like "10**500+100" or "2**64-1000" or "2 * nth_prime(560)".
Additionally, if END starts with '+' then it is assumed to add to START.
If only one number is given, primes up to that number are shown (START = 0).
General options:
--help displays this help message
Filter options, which will cause the list of primes to be further filtered
to only those primes additionally meeting these conditions:
--twin Twin p+2 is prime
--triplet Triplet p+6 and (p+2 or p+4) are prime
--quadruplet Quadruplet p+2, p+6, and p+8 are prime
--cousin Cousin p+4 is prime
--sexy Sexy p+6 is prime
--safe Safe (p-1)/2 is also prime
--sophie Sophie Germain 2p+1 is also prime
--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)
--cuban1 Cuban (y+1) p = (x^3 - y^3)/(x-y), x=y+1
--cuban2 Cuban (y+2) p = (x^3 - y^3)/(x-y), x=y+2
--pnp1 Primorial+1 p#+1 is prime
--pnm1 Primorial-1 p#-1 is prime
--euclid Euclid pn#+1 is prime
--circular Circular all digit rotations of p are prime
--panaitopol Panaitopol p = (x^4-y^4)/(x^3+y^3) for some x,y
--provable Ensure all primes are provably prime
Note that options can be combined, e.g. display only safe twin primes.
In all cases involving multiples (twin, triplet, etc.), the value returned
is p -- the least value of the set.
EOU
}