#!/usr/bin/env perl
use 5.010;
use warnings;
use strict;
use Carp qw/croak/;
use Time::HiRes qw/time/;
use Benchmark qw/cmpthese/;
use Math::Prime::Util qw/is_prime forprimes forcomposites/;
use Math::Primality;
use Math::Prime::XS;
use Math::Pari;
# Other options include:
# Math::Prime::FastSieve
# Math::GMPz (some undocumented functions)
# Toby's example
sub isPrime1 (_)
{
my $num = shift;
return if $num == 1; # 1 is not prime
croak "usage: isPrime(NATURAL NUMBER)"
unless $num =~ /^[1-9][0-9]*$/;
for my $div (2 .. sqrt $num)
{
return if $num % $div == 0;
}
return 1;
}
# Johngg's conversion of Toby's example to a mod-2 wheel
sub tobyinkPlus (_)
{
my $toTest = shift;
return 0 if $toTest == 1;
return 1 if $toTest == 2;
return 0 unless $toTest % 2;
my $sqrtLimit = sqrt $toTest;
for ( my $div = 3; $div <= $sqrtLimit; $div += 2 )
{
return 0 unless $toTest % $div;
}
return 1;
}
# Slightly faster version of the trivial method
sub isPrime1a (_)
{
my($n) = @_;
croak "isPrime must be a positive integer" unless $n !~ tr/0123456789//c;
return (0,0,1,1,0,1,0,1,0,0)[$n] if $n < 10;
for my $div (2..sqrt $n) {
return 0 unless $n % $div;
}
1;
}
# A Mod-2 (i.e. odds-only) trial division
sub isPrime1b (_)
{
my($n) = @_;
croak "isPrime must be a positive integer" unless $n !~ tr/0123456789//c;
return (0,0,1,1,0,1,0,1,0,0)[$n] if $n < 10;
return unless $n & 1;
my $limit = int(sqrt($n));
for (my $div = 3; $div <= $limit; $div += 2) {
return 0 unless $n % $div;
}
1;
}
# Trial division on a mod-6 wheel, so we skip multiples of 2 and 3
sub isPrime1c (_)
{
my($n) = @_;
croak "isPrime must be a positive integer" unless $n !~ tr/0123456789//c;
return (0,0,1,1,0,1,0,1,0,0)[$n] if $n < 10;
return 0 unless $n & 1;
return 0 unless $n % 3;
my $limit = int(sqrt($n));
for (my $div = 5; $div <= $limit; $div += 6) {
return 0 unless $n % $div;
return 0 unless $n % ($div+2);
}
1;
}
# Johngg's SoE method
sub isPrime2
{
my $toTest = shift;
my $sqrtLimit = sqrt $toTest;
my $sieve = q{};
vec( $sieve, 0 , 1 ) = 1;
vec( $sieve, 1 , 1 ) = 1;
vec( $sieve, $toTest, 1 ) = 0;
my $marker = 1;
while ( $marker < $sqrtLimit )
{
my $possPrime = $marker + 1;
$possPrime ++ while vec( $sieve, $possPrime, 1 );
my $fill = 2 * $possPrime;
while ( $fill <= $toTest )
{
vec( $sieve, $fill, 1 ) = 1;
$fill += $possPrime;
}
last if vec( $sieve, $toTest, 1 );
$marker = $possPrime;
}
return not vec($sieve, $toTest, 1);
}
# Using better sieve. Still a bad way to do primality.
sub isPrime2a
{
my($n) = @_;
croak "isPrime must be a positive integer" unless $n !~ tr/0123456789//c;
return (0,0,1,1,0,1,0,1,0,0)[$n] if $n < 10;
return 0 unless $n & 1;
my ($sieve, $i, $limit, $s_end) = ( '', 3, int(sqrt($n)), $n >> 1 );
while ( $i <= $limit ) {
for (my $s = ($i*$i) >> 1; $s <= $s_end; $s += $i) {
vec($sieve, $s, 1) = 1;
}
do { $i += 2 } while vec($sieve, $i >> 1, 1) != 0;
}
vec($sieve, $n >> 1, 1) ? 0 : 1;
}
warn "Validating tests on small primes\n";
forprimes {
die "isPrime1 fail $_" unless isPrime1($_);
die "isPrime1a fail $_" unless isPrime1a($_);
die "isPrime1b fail $_" unless isPrime1b($_);
die "isPrime1c fail $_" unless isPrime1c($_);
die "is_prime fail $_" unless is_prime($_);
die "isPrime2 fail $_" unless $_ > 1000 || isPrime2($_);
die "isPrime2a fail $_" unless $_ > 1000 || isPrime2a($_);
} 100_000;
warn "Validating tests on small composites\n";
forcomposites {
die "isPrime1 fail $_" if isPrime1($_);
die "isPrime1a fail $_" if isPrime1a($_);
die "isPrime1b fail $_" if isPrime1b($_);
die "isPrime1c fail $_" if isPrime1c($_);
die "is_prime fail $_" if $_ < 1000 && is_prime($_);
die "isPrime2a fail $_" if $_ < 1000 && isPrime2a($_);
} 100_000;
my @small = qw/75 169 229 367 369 8794 9227 10807 11939 14803 19937 19938 39783 47083/;
my @medium = qw/526619 737159 872353 3022457 3932021 10180609 46051073 286914127 305669779 334473319 471623353 806431723 1615845311 1817131787/;
my @large = qw/15242882507 15442725479 33130600963 83197255529 86460333013 136073779057 626171233423 891666642109 1669672223953 1754764041599 3493560177931 3794553266983 43249947787433 112325626292951/;
my @ttest = qw/75 169 229 367 369 8794 9227 10807 11939 14803 19937 19938 39783 47083 199933 75640351 760149189769/;
say q[ TOBYINK DIV2 DIV6 MPU];
say q[ NUMBER : RES TIME RES TIME RES TIME RES TIME];
for my $num (qw(75 169 229 367 369 8794 9227 10807 11939 14803 19937 19938 39783 47083 199933 75640351 760149189769 635921898906263))
{
printf('%15d :', $num);
timeit($num, $_) for \&isPrime1, \&isPrime1b, \&isPrime1c, \&is_prime;
print "\n";
}
sub timeit {
my ($n, $f) = @_;
my $start = time();
print($f->($n) ? " Y " : " N ");
my $end = time();
printf '%.06f', ($end - $start);
}
say "\nSmall inputs (2-5 digits)\n";
cmpthese(-1, {
"Toby Div1" => sub { isPrime1($_) for @small },
"JohnGG Div2" => sub { tobyinkPlus($_) for @small },
"DJ Div1" => sub { isPrime1a($_) for @small },
"DJ Div2" => sub { isPrime1b($_) for @small },
"DJ Div6" => sub { isPrime1c($_) for @small },
"JohnGG SoE" => sub { isPrime2($_) for @small },
"DJ SoE" => sub { isPrime2a($_) for @small },
"Module:MPU" => sub { is_prime($_) for @small },
"Module:MP" => sub { Math::Primality::is_prime($_) for @small },
"Module:MPXS" => sub { Math::Prime::XS::is_prime($_) for @small },
"Module:Pari" => sub { Math::Pari::isprime($_) for @small },
});
say "\nMedium inputs (6-10 digits)\n";
cmpthese(-2, {
"Toby Div1" => sub { isPrime1($_) for @medium },
"JohnGG Div2" => sub { tobyinkPlus($_) for @medium },
"DJ Div1" => sub { isPrime1a($_) for @medium },
"DJ Div2" => sub { isPrime1b($_) for @medium },
"DJ Div6" => sub { isPrime1c($_) for @medium },
"JohnGG SoE" => sub { isPrime2($_) for @small },
"DJ SoE" => sub { isPrime2a($_) for @small },
"Module:MPU" => sub { is_prime($_) for @medium },
"Module:MP" => sub { Math::Primality::is_prime($_) for @medium },
"Module:MPXS" => sub { Math::Prime::XS::is_prime($_) for @medium },
"Module:Pari" => sub { Math::Pari::isprime($_) for @medium },
});
say "\nMedium inputs (10-15 digits)\n";
cmpthese(-8, {
"Toby Div1" => sub { isPrime1($_) for @large },
"JohnGG Div2" => sub { tobyinkPlus($_) for @large },
"DJ Div1" => sub { isPrime1a($_) for @large },
"DJ Div2" => sub { isPrime1b($_) for @large },
"DJ Div6" => sub { isPrime1c($_) for @large },
"Module:MPU" => sub { is_prime($_) for @large },
"Module:MP" => sub { Math::Primality::is_prime($_) for @large },
"Module:MPXS" => sub { Math::Prime::XS::is_prime($_) for @large },
"Module:Pari" => sub { Math::Pari::isprime($_) for @large },
});
####
Validating tests on small primes
Validating tests on small composites
TOBYINK DIV2 DIV6 MPU
NUMBER : RES TIME RES TIME RES TIME RES TIME
75 : N 0.000004 N 0.000002 N 0.000001 N 0.000002
169 : N 0.000002 N 0.000002 N 0.000001 N 0.000000
229 : Y 0.000004 Y 0.000001 Y 0.000001 Y 0.000000
367 : Y 0.000004 Y 0.000002 Y 0.000002 Y 0.000000
369 : N 0.000002 N 0.000001 N 0.000001 N 0.000000
8794 : N 0.000002 N 0.000001 N 0.000001 N 0.000000
9227 : Y 0.000008 Y 0.000005 Y 0.000004 Y 0.000000
10807 : N 0.000008 N 0.000005 N 0.000004 N 0.000000
11939 : Y 0.000009 Y 0.000005 Y 0.000005 Y 0.000000
14803 : N 0.000009 N 0.000005 N 0.000004 N 0.000000
19937 : Y 0.000011 Y 0.000006 Y 0.000005 Y 0.000001
19938 : N 0.000002 N 0.000001 N 0.000001 N 0.000001
39783 : N 0.000002 N 0.000001 N 0.000001 N 0.000000
47083 : N 0.000014 N 0.000008 N 0.000006 N 0.000000
199933 : Y 0.000030 Y 0.000017 Y 0.000014 Y 0.000008
75640351 : Y 0.000554 Y 0.000311 Y 0.000246 Y 0.000001
760149189769 : Y 0.056347 Y 0.031018 Y 0.024714 Y 0.000011
635921898906263 : Y 1.605829 Y 0.904163 Y 0.721130 Y 0.000004
####
Small inputs (2-5 digits)
Rate JohnGG SoE DJ SoE Module:MP Toby Div1 DJ Div1 JohnGG Div2 Module:Pari DJ Div2 DJ Div6 Module:MPXS Module:MPU
JohnGG SoE 16.7/s -- -62% -100% -100% -100% -100% -100% -100% -100% -100% -100%
DJ SoE 43.8/s 163% -- -99% -100% -100% -100% -100% -100% -100% -100% -100%
Module:MP 4567/s 27304% 10325% -- -72% -81% -82% -83% -84% -86% -99% -100%
Toby Div1 16439/s 98537% 37425% 260% -- -32% -37% -39% -42% -51% -97% -98%
DJ Div1 24093/s 144460% 54896% 428% 47% -- -8% -11% -15% -29% -95% -97%
JohnGG Div2 26065/s 156287% 59395% 471% 59% 8% -- -4% -8% -23% -94% -97%
Module:Pari 27049/s 162194% 61642% 492% 65% 12% 4% -- -4% -20% -94% -97%
DJ Div2 28183/s 168995% 64230% 517% 71% 17% 8% 4% -- -17% -94% -97%
DJ Div6 33810/s 202762% 77076% 640% 106% 40% 30% 25% 20% -- -93% -96%
Module:MPXS 472615/s 2835592% 1078696% 10248% 2775% 1862% 1713% 1647% 1577% 1298% -- -50%
Module:MPU 953746/s 5722376% 2176929% 20782% 5702% 3859% 3559% 3426% 3284% 2721% 102% --
Medium inputs (6-10 digits)
Rate JohnGG SoE DJ SoE Toby Div1 DJ Div1 JohnGG Div2 DJ Div2 DJ Div6 Module:MP Module:MPXS Module:Pari Module:MPU
JohnGG SoE 16.7/s -- -62% -78% -86% -87% -88% -91% -99% -99% -100% -100%
DJ SoE 44.0/s 164% -- -43% -62% -65% -69% -75% -97% -99% -99% -100%
Toby Div1 77.0/s 362% 75% -- -34% -39% -45% -57% -95% -98% -99% -100%
DJ Div1 116/s 598% 165% 51% -- -8% -18% -35% -93% -96% -98% -100%
JohnGG Div2 126/s 657% 187% 64% 8% -- -10% -29% -92% -96% -98% -100%
DJ Div2 141/s 746% 221% 83% 21% 12% -- -21% -91% -96% -98% -100%
DJ Div6 179/s 972% 306% 132% 54% 41% 27% -- -89% -95% -97% -100%
Module:MP 1619/s 9617% 3582% 2004% 1292% 1183% 1048% 807% -- -50% -72% -99%
Module:MPXS 3270/s 19519% 7335% 4149% 2710% 2490% 2218% 1731% 102% -- -44% -98%
Module:Pari 5798/s 34686% 13082% 7434% 4883% 4492% 4011% 3146% 258% 77% -- -97%
Module:MPU 189149/s 1134794% 429965% 245681% 162462% 149728% 134012% 105802% 11580% 5685% 3162% --
Medium inputs (10-15 digits)
Rate Toby Div1 DJ Div1 JohnGG Div2 DJ Div2 DJ Div6 Module:MPXS Module:MP Module:Pari Module:MPU
Toby Div1 0.586/s -- -32% -39% -45% -57% -98% -100% -100% -100%
DJ Div1 0.858/s 46% -- -10% -19% -36% -97% -100% -100% -100%
JohnGG Div2 0.957/s 63% 12% -- -10% -29% -96% -100% -100% -100%
DJ Div2 1.07/s 82% 24% 11% -- -21% -96% -100% -100% -100%
DJ Div6 1.35/s 130% 57% 41% 27% -- -95% -100% -100% -100%
Module:MPXS 24.9/s 4140% 2797% 2497% 2233% 1741% -- -94% -99% -100%
Module:MP 389/s 66256% 45241% 40546% 36418% 28718% 1465% -- -84% -100%
Module:Pari 2387/s 407067% 278118% 249307% 223982% 176731% 9504% 514% -- -98%
Module:MPU 114853/s 19593759% 13388435% 12002000% 10783286% 8509436% 462045% 29429% 4712% --