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