| Subject: | prime 1928099 not returned [0.22] |
| Date: | Tue, 02 Nov 2010 11:38:12 +1100 |
| To: | bug-Math-Prime-XS [...] rt.cpan.org |
| From: | Kevin Ryde <user42 [...] zip.com.au> |
In Math::Prime::XS 0.22 on recent i386 debian perl 5.10.1,
sieve_primes() doesn't return the prime 1928099. Eg. foo.pl below
prints "not found".
I think the line
inc = (n * n) - i;
overflows 32-bits at n=802649, causing the loop to mark 1928099 as
composite when it's not.
I think the loop ends up going i=322801, i=1125450, i=1928099. The
first 322801 is before n so no harm, the even number 1125450 is excluded
by the i%2 check, but I believe i=1928099 is mis-marked.
Perhaps the main loop could stop when n > sqrt(hi) and at that point
just push the remaining values from composite[] onto the stack.
I struck this when exercising some code generating sophie germain
primes. I don't know if this particular prime is the only one affected.
Incidentally, I wonder if the test
if (i % 2 == 0)
continue;
could be eliminated by taking inc=2*n. If n and i are both odd then i+n
is even and i+2*n would be the next odd to mark, not needing a check.
(Or is there a risk of overflow near 2^32? I suppose if hi=2^32-1 then
some care would be needed either way. That hi would take what 250Mb of
ram, and wouldn't be very fast, but might be a conceivable size ...)
use strict;
use warnings;
use Math::Prime::XS;
foreach my $p (Math::Prime::XS::sieve_primes (2000000)) {
if ($p == 1928099) {
print "found, ok\n";
exit 0;
}
}
print "not found, bad\n";
exit 1;