Sieve of Erathosthenes speedup

Question

While reading HackerNews, I learned about a data set of the first 50 million primes. This would be a lot, but his website has disclaimer:

n this directory I have the first fifty million primes in blocks of one million. Usually it is faster to run a program on your own computer than to download them, but by popular demand, here they are!

N = 10000

n = np.arange(N)+2
P = []

for M in range(3):

    n = np.arange(N)+2 + M*N
    for p in P:
        n =  n[np.where(n% p != 0)]

    while(len(n) > np.sqrt(N)):
        P += [n[0]]
        n = n[np.where(n% P[-1] != 0)]
        #print len(n),

print len(P), M*N/np.log(M*N), M

When I posted this on HackerNews, someone mentioned this was sub-optimal.

I like this way of using NumPy. This resembles the sieve of Eratosthenes, so you'd think it was fast. However n%p[-1] isn't fast to do, your inner loop ends up taking nearly linear time. At least linear in the final number of primes, which is larger than \$N \log N\$. Hence your algorithm runs in \$N^{3/2}\$ whereas Eratosthenes is \$N \log \log N\$.

How do I fix it? For large values, say N=10**6 it's a bit slow. What is wrong with n%p[-1]?

Answer 1

Sieving in the sense of the old Greek is about eliminating composites by crossing out all multiples of primes (striding additively), and the original SoE does not divide, ever. Trial division is something completely different, and for numbers approaching 2^64 it takes several million times as long as a segmented Sieve of Eratosthenes, even if you trial-divide only by primes.

The memory footprint can be reduced by sieving only the odd numbers, since the only even prime can be pulled out of thin air if need be.

Sieves are usually represented as packed bitmaps rather than one number per array cell, because otherwise memory consumption would be huge and things would be exceedingly slow. More compact forms of representations exist, which extend the 'odds-only' idea to dropping multiples of more small primes. They are called 'wheels', and one of the regulars here has written a brilliant explanation of them.

The famous mod 30 wheel drops 2, 3, and 5, and thus effectively stuffs 30 integers into a single byte. Packed odds-only bitmaps are much less complicated than such wheels though and good enough for many practical purposes.

The biggest speedup - almost an order of magnitude - comes from sieving in small blocks that fit into the L1 cache of the CPU, 32 KBytes or thereabouts. However, this makes division creep back into the code, because it is necessary for computing the start offset for each prime when processing a segment.

The division for recomputing the striding position for a prime for every new segment can be eliminated by remembering the last offsets for each prime from segment to segment (or rather the next offsets). That results in a further speedup.

The last optimisation is called 'presieving', which in practical terms means blasting the bit pattern for the composites of a bunch of small primes all over a bitmap before commencing normal operation at the first prime not in that bunch.

There are two variations: the general one blasts a segment immediately before commencing work on it, thus also warming the caches. The special one blasts the whole sieve in one go, and it is slightly faster overall given certain conditions. Perfect for initialising an auxiliary small factors sieve (up to 2^32) for a sieve that operates up to 2^64.

My pastebin has several standalone .cpp programs demonstrating these optimisations and their effects, so that you can gauge what each optimisation can buy you in terms of speed. Presieving might be slightly more difficult to do in python but all the other optimisations apply directly.

As a rough guide, sieving the full 32-bit range using the plain odds-only Eratosthenes takes 20 to 30 seconds on modern CPUs (single-threaded). If you apply all the optimisations then you can cut that to 2 seconds. The segmented sieve logic is also a precondition for parallel sieving.

EDIT: if you are interested in primes in general and large lists - or fast bulk generation - in particular then have a peek at my topic Checksumming large swathes of prime numbers? (for verification). It mentions quite a few sources, up to 2^64 - 10 * 2^32 in fast bulk, from then on at a rate of about 1 million primes per minute (by having gp/PARI dump forprime()). primos.mat.br has lists up to 10^12, by the way.

EDIT: some resources on how to do Sieve of Eratosthenes in python and who to apply some of the optimisations:

• Gareth Rees covers the basics in his excellent post under Sieve of Eratosthenes - Python
• the python SoE examples on RosettaCode show many variants and some optimisations in action
• the numpy/scipy docs have the lowdown on fancy indexing and dtype (to force arrays to specific data types like unsigned 32-bit integers a.k.a. uint32)
• LiterateProgramming explains the basics of SoE in python at a gentler pace

The keys to speedy sieving in languages like python are two:

• efficient, compact representation of the sieve (sets, bool arrays, uint32/uintp arrays)
• pushing processing into the engine/runtime without incurring gratuitous inefficiencies (e.g. instead of division filter where(n % p != 0), use additive striding like is_composite[p*p::p] = True for a bool array)

Things like presieving pretty much require using uintp arrays to be really effective but sets and boolean arrays are more interesting for exploring the logic of the Sieve vs. the capabilities of python, numpy and scipy.

Note: your use of the division filter turns your code into 100% bona fide trial division, and it would make the old Greek turn in his grave... ;-)

Answer 2

I don't know why my question was reverted... I guess I will post my improvement as an answer:

N = 100000

n = np.arange(N)+2
P = []
M = 2

for m in range(M):
    print m

    n = np.arange(N)+2 + m*N
    for p in P:
        n =  n[np.where(n% p != 0)]

    while( n[0]  < np.sqrt(M*N)):
        P += [n[0]]
        n = n[np.where(n% n[0] != 0)]
        #print len(n),
    P = np.hstack((P,n))

print len(P), M*N/np.log(M*N), M

I tried address some of the problems with my original code:

• I introduce the variable m to iterate over range(M).
• Instead of n[np.where(n% P[-1] != 0)] I now say n[np.where(n% n[0] != 0)]. Comments say there is still problems that I use the % operator.
• The termination condition len(n) > np.sqrt(N) has been replaced with n[0] < np.sqrt(M*N)
• Once I get the primes up to \$\sqrt{M \times N}\$ the rest of the numbers in my list are prime. Hence P = np.hstack((P,n)).