1#!/usr/bin/env perl 2use strict; 3use warnings; 4use Math::Prime::Util qw/:all/; 5use Time::HiRes qw(gettimeofday tv_interval); 6$| = 1; # fast pipes 7 8my $nprimes = shift || 50_000_000; 9 10# 1. forprimes does a segmented sieve and calls us for each prime. This is 11# independent of is_prime and the main sieve. So for each entry let's 12# compare next_prime and prev_prime. 13{ 14 print "Using MPU forprimes to $nprimes\n"; 15 my $start_time = [gettimeofday]; 16 my $nextprint = 5000000; 17 my $n = 0; 18 forprimes { 19 die "next $n not $_" unless next_prime($n) == $_; 20 die "prev $n" unless $n == 0 || prev_prime($_) == $n; 21 $n = $_; 22 if ($n > $nextprint) { print "$n.."; $nextprint += 5000000; } 23 } $nprimes; 24 my $seconds = tv_interval($start_time); 25 my $micro_per_call = ($seconds * 1000000) / (2*prime_count($nprimes)); 26 printf "Success using forprimes to $nprimes. %6.2f uSec/call\n", $micro_per_call; 27} 28 29print "\n"; 30 31# 2. Just like before, but now we'll call prime_precalc first. This makes the 32# prev_prime and next_prime functions really fast since they just look in 33# the cached sieve. 34{ 35 print "Using MPU forprimes to $nprimes with prime_precalc\n"; 36 my $start_time = [gettimeofday]; 37 prime_precalc($nprimes); 38 my $nextprint = 5000000; 39 my $n = 0; 40 forprimes { 41 die "next $n not $_" unless next_prime($n) == $_; 42 die "prev $n" unless $n==0 || prev_prime($_) == $n; 43 $n = $_; 44 if ($n > $nextprint) { print "$n.."; $nextprint += 5000000; } 45 } $nprimes; 46 my $seconds = tv_interval($start_time); 47 my $micro_per_call = ($seconds * 1000000) / (2*prime_count($nprimes)); 48 printf "Success using forprimes/precalc to $nprimes. %6.2f uSec/call\n", $micro_per_call; 49} 50 51print "\n\n"; 52 53# Now do some more comparative timing. 54my @pr = @{primes($nprimes)}; 55my $numpr = scalar @pr; 56prime_memfree(); 57 58{ 59 print "MPU forprimes..."; 60 my $start_time = [gettimeofday]; 61 my $i = 0; 62 forprimes { 63 die "next $_ not ", $pr[$i-1] unless $pr[$i++] == $_; 64 } $nprimes; 65 my $seconds = tv_interval($start_time); 66 my $micro_per_call = ($seconds * 1000000) / (1*prime_count($nprimes)); 67 printf "%8.2f uSec/call\n", $micro_per_call; 68 prime_memfree(); 69} 70{ 71 print "MPU prev/next..."; 72 my $start_time = [gettimeofday]; 73 my $n = 0; 74 foreach my $p (@pr) { 75 my $next = next_prime($n); 76 my $prev = prev_prime($p); 77 die "MPU next($n) is not $p\n" unless $next == $p; 78 die "MPU prev($p) is not $n\n" unless $n==0 || $prev == $n; 79 $n = $next; 80 } 81 my $seconds = tv_interval($start_time); 82 my $micro_per_call = ($seconds * 1000000) / (2*$numpr); 83 printf "%8.2f uSec/call\n", $micro_per_call; 84} 85{ 86 print "MPU precalc prev/next..."; 87 my $start_time = [gettimeofday]; 88 prime_precalc($pr[-1]+1000); 89 my $n = 0; 90 foreach my $p (@pr) { 91 my $next = next_prime($n); 92 my $prev = prev_prime($p); 93 die "MPU next($n) is not $p\n" unless $next == $p; 94 die "MPU prev($p) is not $n\n" unless $n==0 || $prev == $n; 95 $n = $next; 96 } 97 my $seconds = tv_interval($start_time); 98 my $micro_per_call = ($seconds * 1000000) / (2*$numpr); 99 printf "%8.2f uSec/call\n", $micro_per_call; 100 prime_memfree(); 101} 102 103# Math::Prime::FastSieve 104if (eval { require Math::Prime::FastSieve; Math::Prime::FastSieve->import(); Inline->init(); 1; }) { 105 print "Math::Prime::FastSieve......"; 106 my $start_time = [gettimeofday]; 107 my $sieve = Math::Prime::FastSieve::Sieve->new( $pr[-1]+1000 ); 108 my $n = 0; 109 foreach my $p (@pr) { 110 my $next = $sieve->nearest_ge($n+1); 111 my $prev = $sieve->nearest_le($p-1); 112 die "MPFS next($n) is not $p\n" unless $next == $p; 113 die "MPFS prev($p) is not $n\n" unless $n==0 || $prev == $n; 114 $n = $next; 115 } 116 my $seconds = tv_interval($start_time); 117 my $micro_per_call = ($seconds * 1000000) / (2*$numpr); 118 printf "%8.2f uSec/call\n", $micro_per_call; 119} else { 120 print "Math::Prime::FastSieve not installed. Skipping\n"; 121} 122 123# Math::Pari. 124if (eval { require Math::Pari; 1; }) { 125 print "Math::Pari prec/next..."; 126 my @pari_pr = grep { $_ < 5_000_000 } @pr; 127 my $pari_numpr = scalar @pari_pr; 128 my $start_time = [gettimeofday]; 129 my $n = 0; 130 foreach my $p (@pari_pr) { 131 my $next = Math::Pari::nextprime($n+1); 132 my $prev = Math::Pari::precprime($p-1); 133 die "Pari next($n) is not $p\n" unless $next == $p; 134 die "Pari prec($p) is not $n\n" unless $n==0 || $prev == $n; 135 $n = $next; 136 } 137 my $seconds = tv_interval($start_time); 138 my $micro_per_call = ($seconds * 1000000) / (2*$pari_numpr); 139 printf "%8.2f uSec/call\n", $micro_per_call; 140} else { 141 print "Math::Pari not installed. Skipping\n"; 142} 143 144# Math::NumSeq::Primes 145if (eval { require Math::NumSeq::Primes; 1; }) { 146 print "Math::NumSeq::Primes next..."; 147 my $start_time = [gettimeofday]; 148 my $seq = Math::NumSeq::Primes->new(); 149 my $n = 0; 150 foreach my $p (@pr) { 151 my $next = ($seq->next)[1]; 152 die "MNP next($n) is not $p\n" unless $next == $p; 153 $n = $next; 154 } 155 my $seconds = tv_interval($start_time); 156 my $micro_per_call = ($seconds * 1000000) / (1*$numpr); 157 printf "%8.2f uSec/call\n", $micro_per_call; 158} else { 159 print "Math::NumSeq::Primes not installed. Skipping\n"; 160} 161 162# Math::Primality 163if (eval { require Math::Primality; 1; }) { 164 print "Math::Primality prev/next..."; 165 my @mp_pr = grep { $_ < 100_000 } @pr; 166 my $mp_numpr = scalar @mp_pr; 167 my $start_time = [gettimeofday]; 168 my $n = 0; 169 foreach my $p (@mp_pr) { 170 my $next = Math::Primality::next_prime($n); 171 my $prev = ($p == 2) ? 0 : Math::Primality::prev_prime($p); 172 die "MP next($n) is not $p\n" unless $next == $p; 173 die "MP prev($p) is not $n\n" unless $n==0 || $prev == $n; 174 $n = $next; 175 } 176 my $seconds = tv_interval($start_time); 177 my $micro_per_call = ($seconds * 1000000) / (2*$mp_numpr); 178 printf "%8.2f uSec/call\n", $micro_per_call; 179} else { 180 print "Math::Primality not installed. Skipping\n"; 181} 182