| Subject: | primes incorrect for 664580, and slow |
Using 'say primes->take(1e6)->sum' takes almost 2 minutes and gives the wrong answer.
The correctness problem is in primes_gen where n is less than 1e7 but the next prime is higher (e.g. 1e7+19). It runs off the prime array. This could be fixed by making $build sieve enough to accommodate any next prime request in the range, or the simple workaround:
if ($n < 1e7-100) {
...
}
which makes us switch to trial division earlier. Test with 664580 to test faster.
The inefficiencies are more complicated.
(1) use correct SoE:
while (($n += 2) <= int(sqrt($max))) {
if (substr $prime, $n, 1) { init: $i = $n;
substr $prime, $i, 1, 0 while ($i += 2*$n) < $max}}
is a tiny change that gives a > 2x speedup.
(2) Use a faster SoE. See http://rosettacode.org/wiki/Sieve_of_Eratosthenes#Perl where the vector sieve should be at least as fast while using 1/16 the memory. The string sieves shown there are somewhat faster. The second one will cut memory use in half.
(3) we're rebuilding the entire prime list every time we expand, tossing out all the work we did earlier. At 10x growth this isn't terrible (Math::NumSeq::Primes is far worse here), so can probably be ignored.
(4) We could use $build for a aux sieve, then sieve a range. This ought to speed things up a couple orders of magnitude, while keeping memory use very low.
(5) Use faster trial division. A 2x speedup.
my $isprime = sub {
my $n = shift;
return ($n >= 2) if $n < 4;
return unless $n % 2 && $n % 3;
my $sqrtn = int(sqrt($n));
for (my $i = 5; $i <= $sqrtn; $i += 6) {
return unless $n % $i && $n % ($i+2);
}
1;
};
....
else {trial_division:
$n++ until $isprime->($n);
return $n++;
}
(6) use a module. A trivial example:
use ntheory "next_prime";
sub primes () {
$primes_gen ||= do {
my $n = 0;
&iterate(sub { $n = next_prime($n); return $n; })
}
}
removes the need to do any of this ourselves, and as a bonus runs faster. We could also use one of the iterators from that module or a tied array. Math::Prime::FastSieve, Math::Prime::XS, and Math::NumSeq::Primes are other alternatives.