#!/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% -- ```