1#!/usr/bin/env perl
2use strict;
3use warnings;
4use Math::Prime::Util qw/:all/;
5use Math::Prime::Util::PrimeArray;
6use Math::NumSeq::Primes;
7use Math::Prime::TiedArray;
8use Benchmark qw/:all/;
9use List::Util qw/min max/;
10my $count = shift || -2;
11
12my ($s, $nlimit, $ilimit, $expect);
13
14if (1) {
15print '-' x 79, "\n";
16print "summation to 100k, looking for best methods (typically slice)\n";
17$nlimit = 100000;
18$ilimit = prime_count($nlimit)-1;
19$expect = 0; forprimes { $expect += $_ } $nlimit;
20
21cmpthese($count,{
22  'pa fetch'  => sub { $s=0; my $o = tie my @p, "Math::Prime::Util::PrimeArray";
23                       $s += $o->FETCH($_) for 0..$ilimit;
24                       die unless $s == $expect; },
25  'pa index'  => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
26                       $s += $primes[$_] for 0..$ilimit;
27                       die unless $s == $expect; },
28  'pa loop'   => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
29                       for (@primes) { last if $_ > $nlimit; $s += $_; }
30                       die $s unless $s == $expect; },
31  'pa slice'  => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
32                       $s += $_ for @primes[0..$ilimit];
33                       die unless $s == $expect; },
34  'pa each'   => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
35                       # Note: using last inside each is Very Bad Stuff.
36                       while(my(undef,$v) = each @primes) { last if $v > $nlimit; $s += $v; }
37                       die $s unless $s == $expect; },
38  'pa shift'  => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
39                       while ((my $p = shift @primes) <= $nlimit) { $s += $p; }
40                       die unless $s == $expect; },
41});
42}
43
44if (1) {
45print '-' x 79, "\n";
46print "summation to 100k, looking for best MPTA extension (typically ~1000)\n";
47$nlimit = 100000;
48$ilimit = prime_count($nlimit)-1;
49$expect = 0; forprimes { $expect += $_ } $nlimit;
50
51cmpthese($count,{
52  'MPTA'       => sub { $s=0; tie my @primes, "Math::Prime::TiedArray";
53                       $s += $primes[$_] for 0..$ilimit;
54                       die unless $s == $expect; },
55  'MPTA 400'   => sub { $s=0; tie my @primes, "Math::Prime::TiedArray", extend_step => 400;
56                       $s += $primes[$_] for 0..$ilimit;
57                       die unless $s == $expect; },
58  'MPTA 1000'  => sub { $s=0; tie my @primes, "Math::Prime::TiedArray", extend_step => 1000;
59                       $s += $primes[$_] for 0..$ilimit;
60                       die unless $s == $expect; },
61  'MPTA 4000'  => sub { $s=0; tie my @primes, "Math::Prime::TiedArray", extend_step => 4000;
62                       $s += $primes[$_] for 0..$ilimit;
63                       die unless $s == $expect; },
64});
65}
66
67if (1) {
68print '-' x 79, "\n";
69print "summation to 100k\n";
70print "Note: MPU::PrimeArray is about 30x faster than MPTA here.\n";
71print "      Math::NumSeq::Primes is reasonable fast (not random access)\n";
72print "      MPU's forprimes smashes everything else (not random access)\n";
73$nlimit = 100000;
74$ilimit = prime_count($nlimit)-1;
75$expect = 0; forprimes { $expect += $_ } $nlimit;
76
77cmpthese($count,{
78  'primes'    => sub { $s=0; $s += $_ for @{primes($nlimit)}; die unless $s == $expect; },
79  'forprimes' => sub { $s=0; forprimes { $s += $_ } $nlimit;  die unless $s == $expect; },
80  'iterator'  => sub { $s=0; my $it = prime_iterator();
81                       $s += $it->() for 0..$ilimit;
82                       die unless $s == $expect; },
83  'OO iter'   => sub { $s=0; my $it = prime_iterator_object();
84                       $s += $it->iterate() for 0..$ilimit;
85                       die unless $s == $expect; },
86  'pa slice'  => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
87                       $s += $_ for @primes[0..$ilimit];
88                       die unless $s == $expect; },
89  'NumSeq'    => sub { $s=0; my $seq = Math::NumSeq::Primes->new;
90                       while (1) { my($undev,$v) = $seq->next; last if $v > $nlimit; $s += $v; }
91                       die $s unless $s == $expect; },
92  # This was slightly faster than slice or shift
93  'MPTA'      => sub { $s=0; tie my @primes, "Math::Prime::TiedArray", extend_step => 1000;
94                       $s += $primes[$_] for 0..$ilimit;
95                       die unless $s == $expect; },
96});
97}
98
99if (0) {
100print '-' x 79, "\n";
101print "summation to 10M\n";
102print "Note: Math::Prime::TiedArray takes too long\n";
103print "      Math::NumSeq::Primes is now ~2x slower than PrimeArray\n";
104print "      forprimes is still the fastest solution for sequential access\n";
105$nlimit = 10_000_000;
106$ilimit = prime_count($nlimit)-1;
107$expect = 0; forprimes { $expect += $_ } $nlimit;
108
109cmpthese($count,{
110  'primes'    => sub { $s=0; $s += $_ for @{primes($nlimit)}; die unless $s == $expect; },
111  'forprimes' => sub { $s=0; forprimes { $s += $_ } $nlimit;  die unless $s == $expect; },
112  'pa index'  => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
113                       $s += $primes[$_] for 0..$ilimit;
114                       die unless $s == $expect; },
115  'pa loop'   => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
116                       for (@primes) { last if $_ > $nlimit; $s += $_; }
117                       die $s unless $s == $expect; },
118  'pa slice'  => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
119                       $s += $_ for @primes[0..$ilimit];
120                       die unless $s == $expect; },
121  'pa each'   => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
122                       while(my(undef,$v) = each @primes) { last if $v > $nlimit; $s += $v; }
123                       die $s unless $s == $expect; },
124  'pa shift'  => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
125                       while ((my $p = shift @primes) <= $nlimit) { $s += $p; }
126                       die unless $s == $expect; },
127  'numseq'    => sub { $s=0; my $seq = Math::NumSeq::Primes->new;
128                       while (1) { my($undev,$v) = $seq->next; last if $v > $nlimit; $s += $v; }
129                       die $s unless $s == $expect; },
130});
131}
132
133if (1) {
134print '-' x 79, "\n";
135print "Walk primes backwards from 1M\n";
136print "Note: MPTA takes 4x longer than just calling MPU's nth_prime!\n";
137$nlimit = 1_000_000;
138$ilimit = prime_count($nlimit)-1;
139$expect = 0; forprimes { $expect += $_ } $nlimit;
140
141cmpthese($count,{
142  'rev primes'=> sub { $s=0; $s += $_ for reverse @{primes($nlimit)}; die unless $s == $expect; },
143  'nthprime'  => sub { $s=0; $s += nth_prime($_) for reverse 1..$ilimit+1; die unless $s == $expect; },
144  'pa index'  => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
145                       $s += $primes[$_] for reverse 0..$ilimit;
146                       die unless $s == $expect; },
147  'OO iter'   => sub { $s=0; my $it = prime_iterator_object($nlimit);
148                       $s += $it->prev->value() for 0..$ilimit;
149                       die unless $s == $expect; },
150  'tiedarray' => sub { $s=0; tie my @primes, "Math::Prime::TiedArray", extend_step => 1000;
151                       $s += $primes[$_] for reverse 0..$ilimit;
152                       die unless $s == $expect; },
153});
154}
155
156if (1) {
157print '-' x 79, "\n";
158print "Random walk in 1M\n";
159print "MPTA takes about 2 minutes and lots of RAM per iteration.\n";
160srand(29);
161my @rindex;
162do { push @rindex, int(rand(1000000)) } for 1..10000;
163$expect = 0; $expect += nth_prime($_+1) for @rindex;
164
165cmpthese($count,{
166  'nthprime'  => sub { $s=0; $s += nth_prime($_+1) for @rindex; },
167  'pa index'  => sub { $s=0; tie my @primes, "Math::Prime::Util::PrimeArray";
168                       $s += $primes[$_] for @rindex;
169                       die unless $s == $expect; },
170   # Argh!  Is it possible to write a slower sieve than the one MPTA uses?
171  #'tiedarray' => sub { $s=0; tie my @primes, "Math::Prime::TiedArray", extend_step => 10000;
172  #                     $s += $primes[$_] for @rindex;
173  #                     die unless $s == $expect; },
174});
175}
176
177print '-' x 79, "\n";
178