[libmath-prime-util-perl] 05/59: Add some PP benchmarks
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:51 UTC 2015
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 <dana at 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";
+
+ at 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
More information about the Pkg-perl-cvs-commits
mailing list