How to speed up a sieve from the eratosthenes pit list generator

Question

My problem comes directly from the CS circles site. This is the last issue at the bottom of this page , called 'Primed for Takeoff'. The main summary is that they need a list of length 1,000.001, where the index of each element is True if the index is simple, and the index of each element is False if it is not simple.

So, for example, isPrime [13] has the value True. isPrime [14] is false.

The first small part of the isPrime list will look like this:

isPrime = [False, False, True, True, False, True, False, True, False, False, ...] 

My problem is the 7 second time they have. I have a working code below with a lower number 101 for debugging purposes. When I hit it to the required list length of 1000.001, it takes too much time, I finished killing the program in a few minutes. This is my work (length 101), very slow code:

 n = 101 c = 2 isPrime = [False,False] for i in range(2,n): isPrime.append(i) def ifInt(isPrime): for item in isPrime: if type(item) == int: return item for d in range(2,n): if c != None: for i in range(c,n,c): isPrime[i] = False isPrime[c] = True c = ifInt(isPrime) 

Then I found it here . It has faster execution time, but it only displays a list of primes, not a list of n lengths, while list [n] returns True for primes and false otherwise.

I'm not sure that this second bit of code really contains the key to creating a main list of 1000.001 in less than 7 seconds, but this was the most important thing I could find in my research.

 def primes_sieve1(limit): limitn = limit+1 primes = dict() for i in range(2, limitn): primes[i] = True for i in primes: factors = range(i,limitn, i) for f in factors[1:]: primes[f] = False return [i for i in primes if primes[i]==True] print primes_sieve1(101) 

CS scripts are pretty widely used, but I could not find a working version of this for Python. Hope this will be easy to solve for someone.

This question is different from this because I am not trying to quickly create a list of primes, but to create a list of all positive integers from 0 to n, which are denoted as "true" and not "just" False.

Answer 1

Edit: I realized that there are a lot of optimizations on SO, but they are rarely ever explained by others for a simple sieve algorithm, so they make it difficult for beginners or creators of the first algorithm to approach. All solutions will be in python to be on the same page for speed and optimization. These decisions will gradually become faster and more difficult. :)

Vanilla solution

 def primesVanilla(n): r = [True] * n r[0] = r[1] = False for i in xrange(n): if r[i]: for j in xrange(i+i,n,i): r[j]=False return r 

This is a very simple sieve implementation. Please make sure that you understand what is happening above before continuing. The only thing to note is that you start to mark non-prime numbers in i + i instead of i, but this is pretty obvious. (Because you assume that I myself am simple). To make the tests honest, all numbers will be up to 25 million for the list.

 real 0m7.663s user 0m7.624s sys 0m0.036s 

Small improvement 1 (square roots):

I will try to sort them in terms of a direct transition to less direct changes. Please note that we do not need to iterate over n, but only go to the square root of n. The reason is that any composite number under n must have a prime coefficient equal to or equal to the square root of n. When you sit by hand, you will notice that all “invalid” numbers above the square root of n are prime by default.

One more note: you have to be a little careful when the square root is integer, so you have to add it in this case so that it covers it. IE, with n = 49 you want to loop up to 7 inclusive, or you can conclude that 49 is simple.

 def primes1(n): r = [True] * n r[0] = r[1] = False for i in xrange(int(n**0.5+1)): if r[i]: for j in xrange(i+i,n,i): r[j]=False return r real 0m4.615s user 0m4.572s sys 0m0.040s 

Please note that this is pretty fast. When you think about it, you only loop to the square root, so now it will take 25 million top-level iterations to only 5,000 upper levels.

Minor improvement 2 (Skip in the inner loop):

Note that in the inner loop, instead of the beginning of i + i, we can start with i * i. This follows from a very similar argument with respect to the square root, but the big idea is that any composites between me and I * I were already marked with lesser touches.

 def primes2(n): r = [True] * n r[0] = r[1] = False for i in xrange(int(n**0.5+1)): if r[i]: for j in xrange(i*i,n,i): r[j]=False return r real 0m4.559s user 0m4.500s sys 0m0.056s 

Good thing a little disappointing. But hey, it's even faster.

Some significant improvement 3 (even skipping):

The idea here is that we can present all even indices and then skip iterations on 2 in the main loop. After that, we can start the outer loop at 3, and the inner loop can skip instead of 2 * i. (Since switching to i instead means that it will be even, (i + i) (i + i + i + i), etc.)

 def primes3(n): r = [True] * n r[0] = r[1] = False for i in xrange(4,n,2):r[i]=False for i in xrange(3,int(n**0.5+1),2): if r[i]: for j in xrange(i*i,n,2*i): r[j]=False return r real 0m2.916s user 0m2.872s sys 0m0.040s 

Cool improvements 4 (wim idea):

This solution is a pretty advanced trick. The purpose of the slice is faster than the cycle, so in this case the notation of the python fragment is used. r [start: end: skip]

 def primes4(n): r = [True] * n r[0] = r[1] = False r[4::2] = [False] * len(r[4::2]) for i in xrange(3,int(1 + n**0.5),2): if r[i]: r[i*i::2*i] = [False] * len(r[i*i::2*i]) return r 10 loops, best of 3: 1.1 sec per loop 

Minor improvement 5

Note that python resells r [4 :: 2] when calculating the length, so it takes quite a bit of extra time, since all we need is a length calculation. To do this, we use some nasty math.

 def primes5(n): r = [True] * n r[0] = r[1] = False r[4::2] = [False] * ((n+1)/2-2) for i in xrange(3,int(1 + n**0.5),2): if r[i]: r[i*i::2*i] = [False] * ((n+2*i-1-i*i)/(2*i)) return r 10 loops, best of 3: 767 msec per loop 

Assign acceleration (Padraic Cunningham) Notice that we assign an array with all True, and then set half (evens) to False. In fact, we can start with a boolean array that alternates.

 def primes6(n): r = [False,True] * (n//2)+[True] r[1],r[2]=False, True for i in xrange(3,int(1 + n**0.5),2): if r[i]: r[i*i::2*i] = [False] * ((n+2*i-1-i*i)/(2*i)) return r 10 loops, best of 3: 717 msec per loop 

Don't quote me on this, but I think that without any nasty mathematical methods there are no obvious improvements to this latest version. One of the nice features that I tried, but did not turn out to be faster, notes that primes other than 2,3 should look like 6k + 1 or 6k-1. (Note that if it is 6k, then divisible by 6, 6k + 2 | 2, 6k + 3 | 3, 6k + 4 | 2, 6k + 5, compare with -1 mod 6. This suggests that we can skip 6 each time and check both sides, either due to a poor implementation on my side or from inside python elements I couldn’t find any significant speed increase .:(

Answer 2

The first thing I saw is how you create the source list (loop and add) is inefficient and unnecessary. You can simply add lists instead of loops and add each item.

The second thing I saw is that the type check you do is not needed, this function call is expensive, and you can reorganize to avoid it.

Finally, I think that the “big thing” that you can get in any sieve implementation uses the slice assignment. You must cross out all the factors in one stroke, and not loop. Example:

 from math import sqrt def primes(n): r = [True] * n r[0] = r[1] = False r[4::2] = [False] * len(r[4::2]) for i in xrange(int(1 + sqrt(n))): if r[i]: r[3*i::2*i] = [False] * len(r[3*i::2*i]) return r 

Note. I also have a few other tricks:

• Avoid half the work by immediately striking out even numbers.
• only iteration is required up to sqrt length

On my shitty low-power macbook, this code can generate a list of 1000,001 in about 75 milliseconds:

 >>> timeit primes(1000001) 10 loops, best of 3: 75.4 ms per loop 

Answer 3

Some timings are shown in python2, and the 3 wim-approach is much faster, it can be slightly optimized by the way the list is created:

 def primes_wim_opt(n): r = [False, True] * (n // 2) r[0] = r[1] = False r[2] = True for i in xrange(int(1 + n ** .5)): if r[i]: r[3*i::2*i] = [False] * len(r[3*i::2*i]) return r 

Python2 timeouts:

 In [9]: timeit primesVanilla(100000) 10 loops, best of 3: 25.7 ms per loop In [10]: timeit primes_wim(100000) 100 loops, best of 3: 3.59 ms per loop In [11]: timeit primes1(100000) 100 loops, best of 3: 14.8 ms per loop In [12]: timeit primes_wim_opt(100000) 100 loops, best of 3: 2.18 ms per loop In [13]: timeit primes2(100000) 100 loops, best of 3: 14.7 ms per loop In [14]: primes_wim(100000) == primes_wim_opt(100000) == primes(100000) == primesVanilla(100000) == primes2(100000) Out[14]: True 

Timing for python3, where using the same functions just changes to a range:

 In [76]: timeit primesVanilla(100000) 10 loops, best of 3: 22.3 ms per loop In [77]: timeit primes_wim(100000) 100 loops, best of 3: 2.92 ms per loop In [78]: timeit primes1(100000) 100 loops, best of 3: 10.9 ms per loop In [79]: timeit primes_wim_opt(100000) 1000 loops, best of 3: 1.88 ms per loop In [80]: timeit primes2(100000) 100 loops, best of 3: 10.3 ms per loop In [81]: primes_wim(100000) == primes_wim_opt(100000) == primes(100000) == primesVanilla(100000) == primes2(100000) Out[80]: True 

it can be optimized further, instead of using the len range / xrange instead of slicing:

 def primes_wim_opt(n): is_odd = n % 2 & 1 r = [False, True] * (n // 2 + is_odd) r[0] = r[1] = False r[2] = True for i in range(int(1 + n ** .5)): if r[i]: r[3*i::2*i] = [False] * len(range(3*i,len(r), 2 * i)) return r 

Python3 removes a good piece:

 In [16]: timeit primes_wim_opt_2(100000) 1000 loops, best of 3: 1.38 ms per loop 

And the same for python2 using xrange:

 In [10]: timeit primes_wim_opt_2(100000) 1000 loops, best of 3: 1.60 ms per loop 

Using (((n - 3 * i) // (2 * i)) + 1) should also work:

 def primes_wim_opt_2(n): is_odd = n % 2 & 1 r = [False, True] * ((n // 2) + is_odd) r[0] = r[1] = False r[2] = True for i in range(int(1 + n ** .5)): if r[i]: r[3*i::2*i] = [False] * (((n - 3 * i) // (2 * i)) + 1) return r 

This is very slightly faster:

 In [12]: timeit primes_wim_opt_2(100000) 1000 loops, best of 3: 1.32 ms per loop In [6]: timeit primes5(100000) 100 loops, best of 3: 2.47 ms per loop 

You can also start with 3 and 2:

 def primes_wim_opt_2(n): r = [False, True] * (n // 2) r[0] = r[1] = False r[2] = True for i in range(3, int(1 + n ** .5),2): if r[i]: r[3*i::2*i] = [False] * (((n - 3 * i) // (2 * i)) + 1) return r 

What happens faster:

 In [2]: timeit primes_wim_opt_2(100000) 1000 loops, best of 3: 1.10 ms per loop 

python2:

 In [2]: timeit primes_wim_opt_2(100000) 1000 loops, best of 3: 1.29 ms per loop