# Sieve of Erathosthenes speedup

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]?

• What is your program supposed to compute? For N = 20 the final value of P is [2, 3, 5, 7, 23, 43] which is neither the list of primes up to 20 nor the list of the first 20 primes. – Martin R Nov 12 '14 at 15:07
• @MartinR if you put for M in range(1) then P should return the prime numbers of up to N. In this case, I was trying to compute primes up to 3*N, carrying over all the primes computed previously. – john mangual Nov 12 '14 at 15:14
• With for M in range(1) and N = 20 I get P = [2, 3, 5, 7]. – Martin R Nov 12 '14 at 15:16

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:

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... ;-)

• This is most interesting – but is it a review of the code in the question? – Martin R Nov 12 '14 at 15:27
• Apart from the P.S. about prime number sources it says where the code in question needs to go in order to achieve certain performance characteristics, without specifying the particalur means of achieving a certain effect (e.g. array of integers for simulating a bitmap and so on). In that sense it should be helpful in progressing from the current state of the reviewed code to where its author wants it to be. Mapping the solution space, so to speak. A 100% pure review would have to restrict itself to mentioning that the code does not in fact do the Sieve of Eratosthenes and that division is bad. – DarthGizka Nov 12 '14 at 15:35
• @DarthGizka please look at my answer, this version should print the primes now... I had posted it earlier but my own question got edited – john mangual Nov 13 '14 at 2:48
• why is_composite[pp::p] Can u explain this point ? This doesn't help filter in faster way as it produce the multipliers of p, you need another step like (set(N)-N[pp::p]) to filter multipliers of p from gloabl list N. This is to replace n[np.where(n% p != 0)] – user3378649 Nov 13 '14 at 7:06
• The logic of the sieve becomes a bit cleaner and more natural if a mark in the sieve array signifies compositeness; instead of is_prime = [False] * 2 + [True] * (N - 1) you use is_composite = [False] * (N + 1). Simpler, cleaner, and it goes well with implicit initialisation to 0/False/nil in many languages. – DarthGizka Nov 13 '14 at 7:06

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)).
• Revising the question — even if you just append a newer version of the code to the question — causes confusion since older answers will discuss the old code and newer answers will discuss the new code. You're welcome to post your new code as an answer, but your answer cannot be a code dump — you need to explain what you changed and why. See a deeper explanation here of our procedures. – 200_success Nov 13 '14 at 6:50
• @ john mangual: your code still uses trial division instead of additive sieving as in the Sieve of Eratosthenes. As long as you keep the division filter your code will remain horribly inefficient, regardless of any numpy/scipy firepower that you bring to bear. – DarthGizka Nov 13 '14 at 7:42
• @DarthGizka I don't understand these terms additive filter and division filter. Which lines are you referring to? – john mangual Nov 13 '14 at 11:01
• n[np.where(n% n[0] != 0)] asks python to perform a trial division n[i] % n[0] on each element of n and to keep only those where the remainder is nonzero. Divisions are extremely expensive. Eratosthenes devised a way of eliminating multiples of primes (a.k.a. 'composites') without any division whatsoever, which consists of writing down all the numbers in order - to make them indexable, in our parlance - and to eliminate multiples of a prime p by going to the (p*p)th number, crossing it out, striding forward by p (to the (p*p + p)th number now), crossing it out, striding forward by p again – DarthGizka Nov 13 '14 at 11:42
• The 'striding forward by p' is where the 'additive' thing comes in, as it is just updating an index i by adding p (or i += p). – DarthGizka Nov 13 '14 at 11:46