This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.10 in repository libmath-prime-util-perl.
commit 1362a6b36a6ec352c3be287961ba0bbc90b2a78b Author: Dana Jacobsen <d...@acm.org> Date: Thu Jun 28 10:24:44 2012 -0600 Add some PP benchmarks --- examples/bench-pp-count.pl | 499 +++++++++++++++++++++++++++++++++++++++++++ examples/bench-pp-isprime.pl | 215 +++++++++++++++++++ examples/bench-pp-sieve.pl | 479 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 1193 insertions(+) diff --git a/examples/bench-pp-count.pl b/examples/bench-pp-count.pl new file mode 100644 index 0000000..45e3151 --- /dev/null +++ b/examples/bench-pp-count.pl @@ -0,0 +1,499 @@ +#!/usr/bin/env perl +use strict; +use warnings; + +use Benchmark qw/:all/; +#use Devel::Size qw/total_size/; +#use Math::Prime::Util; +#use Math::Prime::FastSieve; +#*mpu_erat = \&Math::Prime::Util::erat_primes; +#*fs_erat = \&Math::Prime::FastSieve::primes; + +my $upper = shift || 8192; +my $count = shift || -1; +my $countarg; + +#atkin2(100); exit(0); + +# Shows sizes for sieving to 100k, and rate/second for sieving to 16k +my $pc_subs = { + "Rosetta 4" => sub { rosetta4($countarg) }, # 25k 60/s + "Atkin MPTA" => sub { atkin($countarg) }, # 3430k 90/s + "Merlyn" => sub { merlyn($countarg)}, # 13k 96/s + "Rosetta 2" => sub { rosetta2($countarg) }, # 13k 109/s + "Atkin 2" => sub { atkin2($countarg) }, # 1669k 110/s + "DO Vec" => sub {daoswald_vec($countarg)}, # 13k 112/s + "Rosetta 3" => sub { rosetta3($countarg) }, # 4496k 165/s + "Rosetta 1" => sub { rosetta1($countarg) }, # 3449k 187/s + "Shootout" => sub { shootout($countarg) }, # 3200k 231/s + "DJ Vec" => sub { dj1($countarg) }, # 7k 245/s + "Scriptol" => sub { scriptol($countarg) }, # 3200k 290/s + "DO Array" => sub {daoswald_array($countarg)},# 3200k 306/s + "DJ Array" => sub { dj2($countarg) }, # 1494k 475/s + "In Many" => sub { inmany($countarg) }, # 2018k 666/s + "DJ String1" => sub { dj3($countarg) }, # 50k 981/s + "DJ String2" => sub { dj4($countarg) }, # 50k 1682/s +# "MPU Sieve" => sub { +# scalar @{mpu_erat(2,$countarg)}; }, # 3k 14325/s +# "MPFS Sieve" => sub { +# scalar @{fs_erat($countarg)}; }, # 7k 14325/s +}; + +my %verify = ( + 10 => 4, + 11 => 5, + 100 => 25, + 112 => 29, + 113 => 30, + 114 => 30, + 1000 => 168, + 10000 => 1229, + 100000 => 9592, +); + +# Verify +while (my($name, $sub) = each (%$pc_subs)) { + while (my($n, $pin) = each (%verify)) { + $countarg = $n; + my $picount = $sub->(); + die "$name ($n) = $picount, should be $pin" unless $picount == $pin; + } +} +print "Done with verification, starting benchmark\n"; + +$countarg = $upper; +cmpthese($count, $pc_subs); + + + +# www.scriptol.com/programming/sieve.php +sub scriptol { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my @flags = (0 .. $max); + for my $i (2 .. int(sqrt($max)) + 1) + { + next unless defined $flags[$i]; + for (my $k=$i+$i; $k <= $max; $k+=$i) + { + undef $flags[$k]; + } + } + #print "scriptol size: ", total_size(\@flags), "\n" if $max > 90000; + my $count = 0; + for my $j (2 .. $max) { + $count++ if defined $flags[$j]; + } + $count; +} + +# http://dada.perl.it/shootout/sieve.perl.html +sub shootout { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my $count = 0; + my @flags = (0 .. $max); + for my $i (2 .. $max) { + next unless defined $flags[$i]; + for (my $k=$i+$i; $k <= $max; $k+=$i) { + undef $flags[$k]; + } + $count++; + } + #print "shootout size: ", total_size(\@flags), "\n" if $max > 90000; + $count; +} + +# http://c2.com/cgi/wiki?SieveOfEratosthenesInManyProgrammingLanguages +sub inmany { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + $max++; + + my @c; + for(my $t=3; $t*$t<$max; $t+=2) { + if (!$c[$t]) { + for(my $s=$t*$t; $s<$max; $s+=$t*2) { $c[$s]++ } + } + } + #print "inmany size: ", total_size(\@c), "\n" if $max > 90000; + my $count = 1; + for(my $t=3; $t<$max; $t+=2) { + $c[$t] || $count++; + } + $count; +} + +# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl +sub rosetta1 { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my $count = 0; #my @primes; + my @tested = (1); + my $j = 1; + while ($j < $max) { + next if $tested[$j++]; + $count++; #push @primes, $j; + for (my $k= $j; $k <= $max; $k+=$j) { + $tested[$k-1]= 1; + } + } + #print "R1 size: ", total_size(\@tested), "\n" if $max > 90000; + $count; #scalar @primes; +} + +# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl +sub rosetta2 { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my $count = 0; #my @primes; + my $nonPrimes = ''; + foreach my $p (2 .. $max) { + unless (vec($nonPrimes, $p, 1)) { + for (my $i = $p * $p; $i <= $max; $i += $p) { + vec($nonPrimes, $i, 1) = 1; + } + $count++; #push @primes, $p; + } + } + #print "R2 size: ", total_size(\$nonPrimes), "\n" if $max > 90000; + $count; #scalar @primes; +} + +# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl +sub rosetta3 { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my $i; + my @s; + my $count = scalar + grep { not $s[ $i = $_ ] and do + { $s[ $i += $_ ]++ while $i <= $max; 1 } + } 2 .. $max; + #print "R3 size: ", total_size(\@s), "\n" if $max > 90000; + $count; #scalar @primes; +} + +# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl +sub rosetta4 { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my $i; + my $s = ''; + my $count = scalar + grep { not vec $s, $i = $_, 1 and do + { (vec $s, $i += $_, 1) = 1 while $i <= $max; 1 } + } 2 .. $max; + #print "R4 size: ", total_size(\$s), "\n" if $max > 90000; + $count; #scalar @primes; +} + +# From Math::Primes::TiedArray +sub atkin { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + return 2 if $max < 5; + + my $sqrt = sqrt($max); + my %sieve; + foreach my $x ( 1 .. $sqrt ) { + + foreach my $y ( 1 .. $sqrt ) { + + my $n = 3 * $x**2 - $y**2; + if ( $x > $y + and $n <= $max + and $n % 12 == 11 ) + { + $sieve{$n} = not $sieve{$n}; + } + + $n = 3 * $x**2 + $y**2; + if ( $n <= $max and $n % 12 == 7 ) { + $sieve{$n} = not $sieve{$n}; + } + + $n = 4 * $x**2 + $y**2; + if ( $n <= $max + and ( $n % 12 == 1 or $n % 12 == 5 ) ) + { + $sieve{$n} = not $sieve{$n}; + } + } + } + # eliminate composites by sieving + foreach my $n ( 5 .. $sqrt ) { + + next unless $sieve{$n}; + + my $k = int(1/$n**2) * $n**2; + while ( $k <= $max ) { + $sieve{$k} = 0; + $k += $n**2; + } + } + $sieve{2} = 1; + $sieve{3} = 1; + #print "Atkin size: ", total_size(\%sieve), "\n" if $max > 90000; + + # save the found primes in our cache + my $count = 0; + foreach my $n ( 2 .. $max ) { + next unless $sieve{$n}; + $count++; + } + $count; +} + +# Naive Sieve of Atkin, basically straight from Wikipedia. +# +# <rant> +# +# First thing to note about SoA, is that people love to quote things like +# "memory use is O(N^(1/2+o(1)))" then proceed to _clearly_ use N bytes in +# their implementation. If your data structures between SoA and SoE are the +# same, then all talk about comparative O(blah..blah) memory use is stupid. +# +# Secondly, assuming you're not Dan Bernstein, if your Sieve of Atkin is +# faster than your Sieve of Eratosthenes, then I strongly suggest you verify +# your code actually _works_, and secondly I would bet you made stupid mistakes +# in your SoE implementation. If your SoA code even remotely resembles the +# Wikipedia code and it comes out faster than SoE, then I _guarantee_ your +# SoE is borked. +# +# SoA does have a slightly better asymptotic operation count O(N/loglogN) vs. +# O(N) for SoE. The Wikipedia-like code that most people use is O(N) so it +# isn't even theoretically better unless you pull lots of stunts like primegen +# does. Even if you do, loglogN is essentially a small constant for most uses +# (it's under 4 for all 64-bit values), so you need to make sure all the rest +# of your overhead is controlled. +# +# Sumarizing, in practice the SoE is faster, and often a LOT faster. +# +# </rant> +# +sub atkin2 { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my @sieve; + + my $sqrt = int(sqrt($max)); + for my $x (1 .. $sqrt) { + for my $y (1 .. $sqrt) { + my $n; + + $n = 4*$x*$x + $y*$y; + if ( ($n <= $max) && ( (($n%12) == 1) || (($n%12) == 5) ) ) { + $sieve[$n] ^= 1; + } + $n = 3*$x*$x + $y*$y; + if ( ($n <= $max) && (($n%12) == 7) ) { + $sieve[$n] ^= 1; + } + $n = 3*$x*$x - $y*$y; + if ( ($x > $y) && ($n <= $max) && (($n%12) == 11) ) { + $sieve[$n] ^= 1; + } + } + } + + for my $n (5 .. $sqrt) { + if ($sieve[$n]) { + my $k = $n*$n; + my $z = $k; + while ($z <= $max) { + $sieve[$z] = 0; + $z += $k; + } + } + } + $sieve[2] = 1; + $sieve[3] = 1; + #print "Atkin size: ", total_size(\@sieve), "\n" if $max > 90000; + + my $count = scalar grep { $sieve[$_] } 2 .. $#sieve; + $count; +} + +# https://github.com/daoswald/Inline-C-Perl-Mongers-Talk/blob/master/primesbench.pl +sub daoswald_array { + my($top) = @_; + return 0 if $top < 2; + return 1 if $top < 3; + $top++; + + my @primes = (1) x $top; + my $i_times_j; + for my $i ( 2 .. sqrt $top ) { + if ( $primes[$i] ) { + for ( my $j = $i; ( $i_times_j = $i * $j ) < $top; $j++ ) { + undef $primes[$i_times_j]; + } + } + } + #print "do_array size: ", total_size(\@primes), "\n" if $top > 90000; + my $count = scalar grep { $primes[$_] } 2 .. $#primes; + $count; +} + +sub daoswald_vec { + my($top) = @_; + return 0 if $top < 2; + return 1 if $top < 3; + + my $primes = ''; + vec( $primes, $top, 1 ) = 0; + my $i_times_j; + for my $i ( 2 .. sqrt $top ) { + if ( !vec( $primes, $i, 1 ) ) { + for ( my $j = $i; ( $i_times_j = $i * $j ) <= $top; $j++ ) { + vec( $primes, $i_times_j, 1 ) = 1; + } + } + } + #print "do_vec size: ", total_size(\$primes), "\n" if $top > 90000; + my $count = scalar grep { !vec( $primes, $_, 1 ) } 2 .. $top ; + $count; +} + +# Merlyn's Unix Review Column 26, June 1999 +# http://www.stonehenge.com/merlyn/UnixReview/col26.html +sub merlyn { + my($UPPER) = @_; + return 0 if $UPPER < 2; + return 1 if $UPPER < 3; + + my $count = 0; + my $sieve = ""; + GUESS: for (my $guess = 2; $guess <= $UPPER; $guess++) { + next GUESS if vec($sieve,$guess,1); + $count++; + for (my $mults = $guess * $guess; $mults <= $UPPER; $mults += $guess) { + vec($sieve,$mults,1) = 1; + } + } + #print "Merlyn size: ", total_size(\$sieve), "\n" if $UPPER > 90000; + $count; +} + + +sub dj1 { + my($end) = @_; + return 0 if $end < 2; + return 1 if $end < 3; + + # vector + my $sieve = ''; + my $n = 3; + while ( ($n*$n) <= $end ) { + my $s = $n*$n; + while ($s <= $end) { + vec($sieve, $s >> 1, 1) = 1; + $s += 2*$n; + } + do { $n += 2 } while vec($sieve, $n >> 1, 1) != 0; + } + #print "DJ1 size: ", total_size(\$sieve), "\n" if $end > 90000; + my $count = 1; + $n = 3; + while ($n <= $end) { + $count++ if !vec($sieve, $n >> 1, 1); + $n += 2; + } + $count; +} + +sub dj2 { + my($end) = @_; + return 0 if $end < 2; + return 1 if $end < 3; + + # array + my @sieve; + my $n = 3; + while ( ($n*$n) <= $end ) { + my $s = $n*$n; + while ($s <= $end) { + $sieve[$s>>1] = 1; + $s += 2*$n; + } + do { $n += 2 } while $sieve[$n>>1]; + } + #print "DJ2 size: ", total_size(\@sieve), "\n" if $end > 90000; + my $count = 1; + $n = 3; + while ($n <= $end) { + $count++ if !$sieve[$n>>1]; + $n += 2; + } + $count; +} + +# ~2x faster than inmany, lots faster than the others. Only loses to dj4, +# which is just this code with a presieve added. +sub dj3 { + my($end) = @_; + return 0 if $end < 2; + return 1 if $end < 3; + $end-- if ($end & 1) == 0; + + # string + my $sieve = '1' . '0' x ($end>>1); + my $n = 3; + while ( ($n*$n) <= $end ) { + my $s = $n*$n; + my $filter_s = $s >> 1; + my $filter_end = $end >> 1; + while ($filter_s <= $filter_end) { + substr($sieve, $filter_s, 1) = '1'; + $filter_s += $n; + } + do { $n += 2 } while substr($sieve, $n>>1, 1); + } + #print "DJ3 size: ", total_size(\$sieve), "\n" if $end > 90000; + my $count = 1 + $sieve =~ tr/0//; + $count; +} + +# 2-3x faster than inmany, 6-7x faster than any of the other non-DJ methods. +sub dj4 { + my($end) = @_; + return 0 if $end < 2; + return 1 if $end < 3; + $end-- if ($end & 1) == 0; + + # string with prefill + my $whole = int( ($end>>1) / 15); + my $sieve = '100010010010110' . '011010010010110' x $whole; + substr($sieve, ($end>>1)+1) = ''; + my $n = 7; + while ( ($n*$n) <= $end ) { + my $s = $n*$n; + my $filter_s = $s >> 1; + my $filter_end = $end >> 1; + while ($filter_s <= $filter_end) { + substr($sieve, $filter_s, 1) = '1'; + $filter_s += $n; + } + do { $n += 2 } while substr($sieve, $n>>1, 1); + } + #print "DJ4 size: ", total_size(\$sieve), "\n" if $end > 90000; + my $count = 1 + $sieve =~ tr/0//; + $count; +} diff --git a/examples/bench-pp-isprime.pl b/examples/bench-pp-isprime.pl new file mode 100644 index 0000000..6d739cf --- /dev/null +++ b/examples/bench-pp-isprime.pl @@ -0,0 +1,215 @@ +#!/usr/bin/env perl +use strict; +use warnings; + +use Benchmark qw/:all/; +use Devel::Size qw/total_size/; +use Math::Prime::Util; +*mpu_isprime = \&Math::Prime::Util::is_prime; + +my $count = shift || -1; + +my @numlist; +my @testnums = (0..1000, 5_000_000 .. 5_001_000, 30037, 20359*41117, 92987*65171, 27361*31249, 70790191, 3211717*9673231); + +my $ip_subs = { + #"Abigail" => sub { my$r;$r=abigail($_) for @numlist; $r;}, + "Rosetta" => sub { my$r;$r=rosetta($_) for @numlist; $r;}, + "Rosetta2"=> sub { my$r;$r=rosetta2($_) for @numlist; $r;}, + "DJ" => sub { my$r;$r=dj($_) for @numlist; $r;}, + "DJ2" => sub { my$r;$r=dj2($_) for @numlist; $r;}, + "DJ3" => sub { my$r;$r=dj3($_) for @numlist; $r;}, + "DJ4" => sub { my$r;$r=dj4($_) for @numlist; $r;}, + "MPU" => sub { my$r;$r=mpu_isprime($_) for @numlist; $r;}, +}; + +my %verify = ( + 0 => 0, + 1 => 0, + 2 => 1, + 3 => 1, + 4 => 0, + 5 => 1, + 6 => 0, + 7 => 1, + 13 => 1, + 20 => 0, + 377 => 0, +70790191 => 1, +); + +# Verify +while (my($name, $sub) = each (%$ip_subs)) { + while (my($n, $v_ip) = each (%verify)) { + @numlist = ($n); +#print "$name($n): ", $sub->(), "\n"; + my $isprime = ($sub->() ? 1 : 0); + die "$name($n) = $isprime, should be $v_ip\n" unless $isprime == $v_ip; + } +} +for my $n (0 .. 50000) { + die "dj($n) != mpu($n)" unless dj($n) == mpu_isprime($n); + die "dj2($n) != mpu($n)" unless dj2($n) == mpu_isprime($n); + die "dj3($n) != mpu($n)" unless dj3($n) == mpu_isprime($n); + die "dj4($n) != mpu($n)" unless dj4($n) == mpu_isprime($n); + die "rosetta($n) != mpu($n)" unless rosetta($n) == mpu_isprime($n)/2; + die "rosetta2($n) != mpu($n)" unless rosetta2($n) == mpu_isprime($n)/2; +} +print "Done with verification, starting benchmark\n"; + +@numlist = @testnums; +cmpthese($count, $ip_subs); + + +sub rosetta { + my $n = shift; + $n % $_ or return 0 for 2 .. sqrt $n; + $n > 1; +} + +sub rosetta2 { + my $p = shift; + if ($p == 2) { + return 1; + } elsif ($p <= 1 || $p % 2 == 0) { + return 0; + } else { + my $limit = sqrt($p); + for (my $i = 3; $i <= $limit; $i += 2) { + return 0 if $p % $i == 0; + } + return 1; + } +} + +# Terrifically clever, but useless for large numbers +sub abigail { + ('1' x shift) !~ /^1?$|^(11+?)\1+$/ +} + +sub dj { + my($n) = @_; + return 0 if $n < 2; # 0 and 1 are composite + return 2 if ($n == 2) || ($n == 3) || ($n == 5); # 2, 3, 5 are prime + # multiples of 2,3,5 are composite + return 0 if (($n % 2) == 0) || (($n % 3) == 0) || (($n % 5) == 0); + + my $q; + foreach my $i (qw/7 11 13 17 19 23 29 31 37 41 43 47 53 59/) { + $q = int($n/$i); return 2 if $q < $i; return 0 if $n == ($q*$i); + } + + my $i = 61; # mod-30 loop + while (1) { + $q = int($n/$i); last if $q < $i; return 0 if $n == ($q*$i); $i += 6; + $q = int($n/$i); last if $q < $i; return 0 if $n == ($q*$i); $i += 4; + $q = int($n/$i); last if $q < $i; return 0 if $n == ($q*$i); $i += 2; + $q = int($n/$i); last if $q < $i; return 0 if $n == ($q*$i); $i += 4; + $q = int($n/$i); last if $q < $i; return 0 if $n == ($q*$i); $i += 2; + $q = int($n/$i); last if $q < $i; return 0 if $n == ($q*$i); $i += 4; + $q = int($n/$i); last if $q < $i; return 0 if $n == ($q*$i); $i += 6; + $q = int($n/$i); last if $q < $i; return 0 if $n == ($q*$i); $i += 2; + } + 2; +} + +sub dj2 { + my($n) = @_; + return 2 if ($n == 2) || ($n == 3) || ($n == 5); # 2, 3, 5 are prime + return 0 if $n < 7; # everything else below 7 is composite + # multiples of 2,3,5 are composite + return 0 if (($n % 2) == 0) || (($n % 3) == 0) || (($n % 5) == 0); + + foreach my $i (qw/7 11 13 17 19 23 29 31 37 41 43 47 53 59/) { + return 2 if $i*$i > $n; + return 0 if ($n % $i) == 0; + } + my $limit = int(sqrt($n)); + + my $i = 61; # mod-30 loop + while (1) { + return 0 if ($n % $i) == 0; $i += 6; last if $i > $limit; + return 0 if ($n % $i) == 0; $i += 4; last if $i > $limit; + return 0 if ($n % $i) == 0; $i += 2; last if $i > $limit; + return 0 if ($n % $i) == 0; $i += 4; last if $i > $limit; + return 0 if ($n % $i) == 0; $i += 2; last if $i > $limit; + return 0 if ($n % $i) == 0; $i += 4; last if $i > $limit; + return 0 if ($n % $i) == 0; $i += 6; last if $i > $limit; + return 0 if ($n % $i) == 0; $i += 2; last if $i > $limit; + } + 2; +} + +sub dj3 { + my($n) = @_; + return 2 if ($n == 2) || ($n == 3) || ($n == 5); # 2, 3, 5 are prime + return 0 if $n < 7; # everything else below 7 is composite + # multiples of 2,3,5 are composite + return 0 if (($n % 2) == 0) || (($n % 3) == 0) || (($n % 5) == 0); + + foreach my $i (qw/7 11 13 17 19 23 29 31 37 41 43 47 53 59/) { + return 2 if $i*$i > $n; + return 0 if ($n % $i) == 0; + } + my $limit = int(sqrt($n)); + + my $i = 61; # mod-30 loop + while (($i+30) <= $limit) { + return 0 if ($n % $i) == 0; $i += 6; + return 0 if ($n % $i) == 0; $i += 4; + return 0 if ($n % $i) == 0; $i += 2; + return 0 if ($n % $i) == 0; $i += 4; + return 0 if ($n % $i) == 0; $i += 2; + return 0 if ($n % $i) == 0; $i += 4; + return 0 if ($n % $i) == 0; $i += 6; + return 0 if ($n % $i) == 0; $i += 2; + } + while (1) { + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 6; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 4; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 2; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 4; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 2; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 4; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 6; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 2; + } + 2; +} + +sub dj4 { + my($n) = @_; + return 2 if ($n == 2) || ($n == 3) || ($n == 5); # 2, 3, 5 are prime + return 0 if $n < 7; # everything else below 7 is composite + # multiples of 2,3,5 are composite + return 0 if (($n % 2) == 0) || (($n % 3) == 0) || (($n % 5) == 0); + + foreach my $i (qw/7 11 13 17 19 23 29/) { + return 2 if $i*$i > $n; + return 0 if ($n % $i) == 0; + } + my $limit = int(sqrt($n)); + + my $i = 31; + while (($i+30) <= $limit) { + return 0 if ($n % $i) == 0; $i += 6; + return 0 if ($n % $i) == 0; $i += 4; + return 0 if ($n % $i) == 0; $i += 2; + return 0 if ($n % $i) == 0; $i += 4; + return 0 if ($n % $i) == 0; $i += 2; + return 0 if ($n % $i) == 0; $i += 4; + return 0 if ($n % $i) == 0; $i += 6; + return 0 if ($n % $i) == 0; $i += 2; + } + while (1) { + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 6; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 4; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 2; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 4; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 2; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 4; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 6; + last if $i > $limit; return 0 if ($n % $i) == 0; $i += 2; + } + 2; +} diff --git a/examples/bench-pp-sieve.pl b/examples/bench-pp-sieve.pl new file mode 100644 index 0000000..5464eb3 --- /dev/null +++ b/examples/bench-pp-sieve.pl @@ -0,0 +1,479 @@ +#!/usr/bin/env perl +use strict; +use warnings; + +use Benchmark qw/:all/; +use Devel::Size qw/total_size/; +use Math::Prime::Util; +use Math::Prime::FastSieve; +*mpu_erat = \&Math::Prime::Util::erat_primes; +*fs_erat = \&Math::Prime::FastSieve::primes; + +my $upper = shift || 8192; +my $count = shift || -1; +my $countarg; +my $sum; + +# This is like counting, but we want an array returned. +# The subs will compute a sum on the results. + +# Times for 100k. +# Vs. MPU sieve, as we move from 8k to 10M: +# Atkin MPTA, Rosetta 3 & 1, Shootout, Scriptol, DO Array, DJ Array, and +# InMany all slow down. Atkin 2 speeds up (from 65x slower to 54x slower). +# The DJ string methods have almost no relative slowdown, so stretch out their +# advantage over the others (In Many, DJ Array, DJ Vec, and DO Array). +my $pc_subs = { + "Rosetta 4" => sub {$sum=0; $sum+=$_ for rosetta4($countarg);$sum;}, # 9/s + "Atkin MPTA"=> sub {$sum=0; $sum+=$_ for atkin($countarg);$sum;}, # 11/s + "Merlyn" => sub {$sum=0; $sum+=$_ for merlyn($countarg);$sum;}, # 15/s + "Rosetta 2" => sub {$sum=0; $sum+=$_ for rosetta2($countarg);$sum; }, # 16/s + "DO Vec" => sub {$sum=0; $sum+=$_ for daos_vec($countarg);$sum;}, # 16/s + "Atkin 2" => sub {$sum=0; $sum+=$_ for atkin2($countarg);$sum; }, # 17/s + "Rosetta 3" => sub {$sum=0; $sum+=$_ for rosetta3($countarg);$sum; }, # 23/s + "Rosetta 1" => sub {$sum=0; $sum+=$_ for rosetta1($countarg);$sum; }, # 26/s + "Shootout" => sub {$sum=0; $sum+=$_ for shootout($countarg);$sum; }, # 30/s + "Scriptol" => sub {$sum=0; $sum+=$_ for scriptol($countarg);$sum; }, # 33/s + "DJ Vec" => sub {$sum=0; $sum+=$_ for dj1($countarg);$sum; }, # 34/s + "DO Array" => sub {$sum=0; $sum+=$_ for daos_array($countarg);$sum;},# 41/s + "DJ Array" => sub {$sum=0; $sum+=$_ for dj2($countarg);$sum; }, # 63/s + "In Many" => sub {$sum=0; $sum+=$_ for inmany($countarg);$sum; }, # 86/s + "DJ String1"=> sub {$sum=0; $sum+=$_ for dj3($countarg);$sum; }, # 99/s + "DJ String2"=> sub {$sum=0; $sum+=$_ for dj4($countarg);$sum; }, # 134/s + "MPFS Sieve"=> sub { # 1216/s + $sum=0; $sum+=$_ for @{fs_erat($countarg)};;$sum; }, + "MPU Sieve" => sub { # 1290/s + $sum=0; $sum+=$_ for @{mpu_erat(2,$countarg)};;$sum; }, +}; + +my %verify = ( + 10 => 17, + 11 => 28, + 100 => 1060, + 112 => 1480, + 113 => 1593, + 114 => 1593, + 1000 => 76127, + 10000 => 5736396, + 100000 => 454396537, +); + +# Verify +while (my($name, $sub) = each (%$pc_subs)) { + while (my($n, $v_pi_sum) = each (%verify)) { + $countarg = $n; + my $pi_sum = $sub->(); + die "$name ($n) = $pi_sum, should be $v_pi_sum" unless $pi_sum == $v_pi_sum; + } +} +print "Done with verification, starting benchmark\n"; + +$countarg = $upper; +cmpthese($count, $pc_subs); + + + +# www.scriptol.com/programming/sieve.php +sub scriptol { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my @flags = (0 .. $max); + for my $i (2 .. int(sqrt($max)) + 1) + { + next unless defined $flags[$i]; + for (my $k=$i+$i; $k <= $max; $k+=$i) + { + undef $flags[$k]; + } + } + return grep { defined $flags[$_] } 2 .. $max; +} + +# http://dada.perl.it/shootout/sieve.perl.html +sub shootout { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my @primes; + my @flags = (0 .. $max); + for my $i (2 .. $max) { + next unless defined $flags[$i]; + for (my $k=$i+$i; $k <= $max; $k+=$i) { + undef $flags[$k]; + } + push @primes, $i; + } + @primes; +} + +# http://c2.com/cgi/wiki?SieveOfEratosthenesInManyProgrammingLanguages +sub inmany { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my @c; + for(my $t=3; $t*$t<=$max; $t+=2) { + if (!$c[$t]) { + for(my $s=$t*$t; $s<=$max; $s+=$t*2) { $c[$s]++ } + } + } + my @primes = (2); + for(my $t=3; $t<=$max; $t+=2) { + $c[$t] || push @primes, $t; + } + @primes; + # grep { $c[$_] } 3 .. $max; +} + +# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl +sub rosetta1 { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my @primes; + my @tested = (1); + my $j = 1; + while ($j < $max) { + next if $tested[$j++]; + push @primes, $j; + for (my $k= $j; $k <= $max; $k+=$j) { + $tested[$k-1]= 1; + } + } + @primes; +} + +# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl +sub rosetta2 { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my @primes; + my $nonPrimes = ''; + foreach my $p (2 .. $max) { + unless (vec($nonPrimes, $p, 1)) { + for (my $i = $p * $p; $i <= $max; $i += $p) { + vec($nonPrimes, $i, 1) = 1; + } + push @primes, $p; + } + } + @primes; +} + +# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl +sub rosetta3 { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my(@s, $i); + grep { not $s[ $i = $_ ] and do + { $s[ $i += $_ ]++ while $i <= $max; 1 } + } 2 .. $max; +} + +# http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl +sub rosetta4 { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my $i; + my $s = ''; + grep { not vec $s, $i = $_, 1 and do + { (vec $s, $i += $_, 1) = 1 while $i <= $max; 1 } + } 2 .. $max; +} + +# From Math::Primes::TiedArray +sub atkin { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + return 2 if $max < 5; + + my $sqrt = sqrt($max); + my %sieve; + foreach my $x ( 1 .. $sqrt ) { + + foreach my $y ( 1 .. $sqrt ) { + + my $n = 3 * $x**2 - $y**2; + if ( $x > $y + and $n <= $max + and $n % 12 == 11 ) + { + $sieve{$n} = not $sieve{$n}; + } + + $n = 3 * $x**2 + $y**2; + if ( $n <= $max and $n % 12 == 7 ) { + $sieve{$n} = not $sieve{$n}; + } + + $n = 4 * $x**2 + $y**2; + if ( $n <= $max + and ( $n % 12 == 1 or $n % 12 == 5 ) ) + { + $sieve{$n} = not $sieve{$n}; + } + } + } + # eliminate composites by sieving + foreach my $n ( 5 .. $sqrt ) { + + next unless $sieve{$n}; + + my $k = int(1/$n**2) * $n**2; + while ( $k <= $max ) { + $sieve{$k} = 0; + $k += $n**2; + } + } + my @primes = (2, 3); + push @primes, grep { $sieve{$_} } 5 .. $max; + @primes; +} + +# Naive Sieve of Atkin, basically straight from Wikipedia. +# +# <rant> +# +# First thing to note about SoA, is that people love to quote things like +# "memory use is O(N^(1/2+o(1)))" then proceed to _clearly_ use N bytes in +# their implementation. If your data structures between SoA and SoE are the +# same, then all talk about comparative O(blah..blah) memory use is stupid. +# +# Secondly, assuming you're not Dan Bernstein, if your Sieve of Atkin is +# faster than your Sieve of Eratosthenes, then I strongly suggest you verify +# your code actually _works_, and secondly I would bet you made stupid mistakes +# in your SoE implementation. If your SoA code even remotely resembles the +# Wikipedia code and it comes out faster than SoE, then I _guarantee_ your +# SoE is borked. +# +# SoA does have a slightly better asymptotic operation count O(N/loglogN) vs. +# O(N) for SoE. The Wikipedia-like code that most people use is O(N) so it +# isn't even theoretically better unless you pull lots of stunts like primegen +# does. Even if you do, loglogN is essentially a small constant for most uses +# (it's under 4 for all 64-bit values), so you need to make sure all the rest +# of your overhead is controlled. +# +# Sumarizing, in practice the SoE is faster, and often a LOT faster. +# +# </rant> +# +sub atkin2 { + my($max) = @_; + return 0 if $max < 2; + return 1 if $max < 3; + + my @sieve; + + my $sqrt = int(sqrt($max)); + for my $x (1 .. $sqrt) { + for my $y (1 .. $sqrt) { + my $n; + + $n = 4*$x*$x + $y*$y; + if ( ($n <= $max) && ( (($n%12) == 1) || (($n%12) == 5) ) ) { + $sieve[$n] ^= 1; + } + $n = 3*$x*$x + $y*$y; + if ( ($n <= $max) && (($n%12) == 7) ) { + $sieve[$n] ^= 1; + } + $n = 3*$x*$x - $y*$y; + if ( ($x > $y) && ($n <= $max) && (($n%12) == 11) ) { + $sieve[$n] ^= 1; + } + } + } + + for my $n (5 .. $sqrt) { + if ($sieve[$n]) { + my $k = $n*$n; + my $z = $k; + while ($z <= $max) { + $sieve[$z] = 0; + $z += $k; + } + } + } + + $sieve[2] = 1; + $sieve[3] = 1; + grep { $sieve[$_] } 2 .. $max; +} + +# https://github.com/daoswald/Inline-C-Perl-Mongers-Talk/blob/master/primesbench.pl +sub daos_array { + my($top) = @_; + return 0 if $top < 2; + return 1 if $top < 3; + $top++; + + my @primes = (1) x $top; + my $i_times_j; + for my $i ( 2 .. sqrt $top ) { + if ( $primes[$i] ) { + for ( my $j = $i; ( $i_times_j = $i * $j ) < $top; $j++ ) { + undef $primes[$i_times_j]; + } + } + } + return grep { $primes[$_] } 2 .. $#primes; +} + +sub daos_vec { + my($top) = @_; + return 0 if $top < 2; + return 1 if $top < 3; + + my $primes = ''; + vec( $primes, $top, 1 ) = 0; + my $i_times_j; + for my $i ( 2 .. sqrt $top ) { + if ( !vec( $primes, $i, 1 ) ) { + for ( my $j = $i; ( $i_times_j = $i * $j ) <= $top; $j++ ) { + vec( $primes, $i_times_j, 1 ) = 1; + } + } + } + return grep { !vec( $primes, $_, 1 ) } 2 .. $top; +} + +# Merlyn's Unix Review Column 26, June 1999 +# http://www.stonehenge.com/merlyn/UnixReview/col26.html +sub merlyn { + my($UPPER) = @_; + return 0 if $UPPER < 2; + return 1 if $UPPER < 3; + + my @primes; + my $sieve = ""; + GUESS: for (my $guess = 2; $guess <= $UPPER; $guess++) { + next GUESS if vec($sieve,$guess,1); + push @primes, $guess; + for (my $mults = $guess * $guess; $mults <= $UPPER; $mults += $guess) { + vec($sieve,$mults,1) = 1; + } + } + @primes; +} + + +sub dj1 { + my($end) = @_; + return 0 if $end < 2; + return 1 if $end < 3; + + # vector + my $sieve = ''; + my $n = 3; + while ( ($n*$n) <= $end ) { + my $s = $n*$n; + while ($s <= $end) { + vec($sieve, $s >> 1, 1) = 1; + $s += 2*$n; + } + do { $n += 2 } while vec($sieve, $n >> 1, 1) != 0; + } + + my @primes = (2); + $n = 3; + while ($n <= $end) { + push @primes, $n if !vec($sieve, $n >> 1, 1); + $n += 2; + } + @primes; +} + +sub dj2 { + my($end) = @_; + return 0 if $end < 2; + return 1 if $end < 3; + + # array + my @sieve; + my $n = 3; + while ( ($n*$n) <= $end ) { + my $s = $n*$n; + while ($s <= $end) { + $sieve[$s>>1] = 1; + $s += 2*$n; + } + do { $n += 2 } while $sieve[$n>>1]; + } + my @primes = (2); + $n = 3; + while ($n <= $end) { + push @primes, $n if !$sieve[$n>>1]; + $n += 2; + } + @primes; +} + +sub dj3 { + my($end) = @_; + return 0 if $end < 2; + return 1 if $end < 3; + $end-- if ($end & 1) == 0; + + # string + my $sieve = '1' . '0' x ($end>>1); + my $n = 3; + while ( ($n*$n) <= $end ) { + my $s = $n*$n; + my $filter_s = $s >> 1; + my $filter_end = $end >> 1; + while ($filter_s <= $filter_end) { + substr($sieve, $filter_s, 1) = '1'; + $filter_s += $n; + } + do { $n += 2 } while substr($sieve, $n>>1, 1); + } + my @primes = (2); + $n = 3-2; + foreach my $s (split("0", substr($sieve, 1), -1)) { + $n += 2 + 2 * length($s); + push @primes, $n if $n <= $end; + } + @primes; +} + +sub dj4 { + my($end) = @_; + return 0 if $end < 2; + return 1 if $end < 3; + $end-- if ($end & 1) == 0; + + # string with prefill + my $whole = int( ($end>>1) / 15); + my $sieve = '100010010010110' . '011010010010110' x $whole; + substr($sieve, ($end>>1)+1) = ''; + my $n = 7; + while ( ($n*$n) <= $end ) { + my $s = $n*$n; + my $filter_s = $s >> 1; + my $filter_end = $end >> 1; + while ($filter_s <= $filter_end) { + substr($sieve, $filter_s, 1) = '1'; + $filter_s += $n; + } + do { $n += 2 } while substr($sieve, $n>>1, 1); + } + my @primes = (2, 3, 5); + $n = 7-2; + foreach my $s (split("0", substr($sieve, 3), -1)) { + $n += 2 + 2 * length($s); + push @primes, $n if $n <= $end; + } + @primes; +} -- 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 Pkg-perl-cvs-commits@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits