``````#!/usr/bin/env perl
use strict;
use warnings;

use Math::Prime::Util qw/prime_iterator is_prime
next_prime nth_prime_upper prime_precalc forprimes/;

my \$count = shift || 20;
my \$method = shift || 'forprimes';
my \$precalc = 0;   # If set, precalc all the values we'll call is_prime on

# Find Sophie Germain primes (numbers where p and 2p+1 are both prime).

# Four methods are shown: forprimes, iter, iter2, and MNS.

# Times for 300k:
#
#                300k               1M
# precalc:
#   forprimes    1.3s   9.0MB       7.1s   21.6MB
#   iter         2.8s   8.7MB      12.6s   21.4MB
#   iter2        1.9s   8.7MB       9.4s   21.4MB
# no precalc:
#   forprimes    1.5s   4.5MB       5.6s   4.5MB
#   iter         9.5s   4.3MB      37.5s   4.3MB
#   iter2        8.5s   4.3MB      33.9s   4.3MB
#   MNS        254.3s  11.3MB   >1500s   >15 MB

if (\$precalc) {
prime_precalc(2 * sg_upper_bound(\$count));
}

if (\$method eq 'forprimes') {

my \$estimate = sg_upper_bound(\$count);
my \$numfound = 0;
forprimes {
if (\$numfound < \$count && is_prime(2*\$_+1)) {
print "\$_\n";
\$numfound++;
}
} \$estimate;
die "Estimate too low" unless \$numfound >= \$count;

} elsif (\$method eq 'iter') {

# Wrap the standard iterator
sub get_sophie_germain_iterator {
my \$p = shift || 2;
my \$it = prime_iterator(\$p);
return sub {
do { \$p = \$it->() } while !is_prime(2*\$p+1);
\$p;
};
}
my \$sgit = get_sophie_germain_iterator();
print \$sgit->(), "\n" for 1 .. \$count;

} elsif (\$method eq 'iter2') {

# Iterate directly using next_prime
my \$prime = 2;
for (1 .. \$count) {
\$prime = next_prime(\$prime) while !is_prime(2*\$prime+1);
print "\$prime\n";
\$prime = next_prime(\$prime);
}

} elsif (\$method eq 'MNS') {

# Use Math::NumSeq
require Math::NumSeq::SophieGermainPrimes;
my \$seq = Math::NumSeq::SophieGermainPrimes->new;
for (1 .. \$count) {
print 0+(\$seq->next)[1];
}

}

# Used for precalc and the forprimes example
sub sg_upper_bound {
my \$count = shift;
my \$nth = nth_prime_upper(\$count);
# For lack of a better formula, do this step-wise estimate.
my \$estimate = (\$count <   5000) ? 150 + int( \$nth * log(\$nth) * 1.2 )
: (\$count <  19000) ? int( \$nth * log(\$nth) * 1.135 )
: (\$count <  45000) ? int( \$nth * log(\$nth) * 1.10 )
: (\$count < 100000) ? int( \$nth * log(\$nth) * 1.08 )
: (\$count < 165000) ? int( \$nth * log(\$nth) * 1.06 )
: (\$count < 360000) ? int( \$nth * log(\$nth) * 1.05 )
: (\$count < 750000) ? int( \$nth * log(\$nth) * 1.04 )
: (\$count <1700000) ? int( \$nth * log(\$nth) * 1.03 )
:                     int( \$nth * log(\$nth) * 1.02 );

return \$estimate;
}
``````