Wow. Thank you.
Nice article, but one minor quibble.
> An attempt to allocate workspace of 4294967295 bytes quickly crashes the program
That is not necessarily true. Unless malloc() has a bug it will return NULL. If the code is properly checking returns, it will not crash (it may exit if it's policy is to give up on out of memory or it may abandon that unit of work and recover).
Samuel, thanks for the correction. It looks as if I can't even code up a bug correctly 
I'm not sure I was the first to use the phrase "ship happens", but I do recall using in 2002 with colleagues and then in my book "Practical Development Environments", p.240, as a subchapter heading.
p225 to be correct.
What I really meant to say was "thanks for an interesting article" and "please don't let this be my 15 minutes of fame!"
Thanks for stepping forward, Matt, and I'm glad you found the article of interest. 2002 is the earliest so far. Unless anyone else knows otherwise ... ?
Interesting article.
A python coding question. Why this:
for i in islice(count(2), sqrt_hi):
if sieve[i] != 0:
remove = slice(i * i, hi, i)
sieve[remove] = zeros[remove]
instead of this:
for i in xrange(2, sqrt_hi+1):
if sieve[i] != 0:
for remove in xrange(i * i, hi, i) :
sieve[remove] = 0
Is it faster to precalculate the indexes and then copy the zeroes from a list of zeroes? The itertools variant seems harder to grok.
And another thing. Implicitly, your code for is_fprime() is:
def is_fprime(n, tries=3):
'''Use Fermat little theorem to guess if n is prime.
'''
xs = (randrange(1, n) for _ in range(tries))
return all(prime_mod_test(x, n) for x in xs)
where prime_mod_test() is:
def prime_mod_test(test, base) :
return (test ** base) % base == test
This code is brutally slow. It takes, well I don't know how long it takes.it ran for 25 minutes and I killed it. However, with this as prime_mod_test():
def prime_mod_test(test, base) :
a = 1
for i in xrange(1, base) :
a = (a * test) % base
if a == 1 :
return (base - 1) % i == 0
return False
the runtime is 83 seconds. The observation is that if you hit a remainder of 1 before doing all the powers, then you will cycle back to 1 eventually, and the only question is whether you will cycle on the the base-1 iteration. Well, you'll do that if the loop iterator is a multiple of base-1, so test for (base-1) % i == 0.
Hugh, thanks for your comments and well done on figuring out how to workaround the code-formatting issues! To answer the first point, about xrange, islice and slicing: I've simply never liked xrange, even though I recognise the need for it, and I'm pleased Python 3.0 removes this need. I'm starting to think about trying to make my Python 2.* code work unchanged with Python 3.0.
I certainly was not concerned about efficiency or indeed memory use when I used a slice assignment to sieve out multiples of i, but I do like the use of a single assignment to sieve out all these multiples.
I'll respond separately to your second point.
Hugh, to respond to the second point: yes, using is_fprime() repeatedly as I do in the article is indeed very slow, though it may have some use in spotting if a single large number is composite. Whereas the classic Eratosthenes sieve finds primes in order, with each result building on the ones before. Your speedy variant on the Fermat inspired test is ingenious but quite subtle.
Really the code here is orthogonal to my main point, that algorithms which are almost right can be hard to detect and hard to fix. However much we tweak the is_fprime() implementation, it will still return false positives for the Carmichael numbers. I wanted my code to express the ideas in the article as simply as possible, and let the computer feel the (considerable) strain.
Interestingly, the Rabin-Miller primality test uses a similar idea but can be made correct. Of course, finding large prime numbers turns out to be a very important and specialised problem, and I'm well out of my depth here already!
Thanks again for your comments.
Commenting by HaloScan
