8
\$\begingroup\$

I wrote a version of sieve of eratosthenes:

number =100000000
def eratos_sieve(n:"find prime below this"=None):
        if not(n):
            n = number
        numbers = [x for x in range(0, n+1)]
        current_prime = 2
        sqrt_limit = int(math.sqrt(n))+1
        while numbers[current_prime] <= sqrt_limit:
            #(len(numbers)-1) is the actual size since
            #since, len([0,...,n])=n+1
            for i in range(current_prime, ((len(numbers)-1)//current_prime)+1):
                numbers[current_prime*i]=-1
            current_prime += 1
            for i in range(current_prime, sqrt_limit):
                if numbers[i]>0:
                    current_prime=i
                    break


        #this part is optional
        j = 0
        for i in range(len(numbers)):
            if numbers[i]>0:
                numbers[j]=numbers[i]
                j += 1
        #uptil here if later primes are to be found using indices
        #else all primes upto n are returned as a single list
        return numbers[1:j]

It seems to perform fast enough for my purposes. I'd like to get some review about it's performance since I'm a novice. Also, I'd like someone with enough memory (if possible) to run the below two tests:

import time
import math
num1 = 100000000 #100 million
num2 = 1000000000 #1 billion
start = time.perf_counter()
eratos_sieve(num1)
total_100_m = time.perf_counter()-start()

start = time.perf_counter()
eratos_sieve(num2)
total_1_b = time.perf_counter()-start()

print("For 100m it takes: ", total_100_m)
print("For 1b   it takes: ", total_1_b)
\$\endgroup\$
8
\$\begingroup\$

Four* letters to rule them all

The Sieve of Eratosthenes is an algorithm which heavily relies on loops. Unfortunately, Python's convenient scripty nature comes at a cost: it's not terribly fast when it comes to loops, so it's best to avoid them.

However, this is not always possible, or one, like me in this case, is not so much into algorithms to transform them into another form. Enter numba. numba is a just-in-time compiler for Python code. The just-in-time compiler can transform plain Python code into "native" code for your system. So what are those magical four letters I was talking about? It's @jit.

This simple decorator and no further code changes brought the time for n = 100 000 000 from 36s down to about 3s on my machine. The larger test case did not work out to well with my 16GB RAM and some other things to do apart from Code Review.

from numba import jit

@jit
def eratos_sieve(n:"find prime below this"=100000000):
    # ... the rest of your code here

You can also have a look at this generator based implementation of the algorithm presented on Stack Overflow. The implementation can benefit from just-in-time compilation as well. A quick test using the SO approach and @jit for the previous n brought the execution time down to just below 2s on the same machine as before.

What numba cannot fix for you

There are a few aspects with respect to style and idiomatic language features (often called pythonic code) that you have to/can improve yourself.

While

number =10000000
def eratos_sieve(n:"find prime below this"=None):
    if not(n):
        n = number
    ...

might be working and syntactically valid, I find it very, very uncommon and would not recommend it. A more "pythonic" way to express this would be

def eratos_sieve(n: int = 100000000):
    """Find all primes below a given whole number n"""
    ...

This makes use of Python 3's type hinting and the recommended way to write function documentation using docstrings.

The way you construct your candidate list using numbers = [x for x in range(0, n+1)] is also a little bit complicated, and may be simplified to numbers = list(range(n+1)). 0 is the implicit starting value for range(...) and since range(...) already returns a generator you can simply pass it to a list constructor and work with it.

You should also have a look at the official Style Guide for Python Code (PEP8) and follow it's recommendations regarding whitespace within expressions and statements.

If you intend to supply your testing code together with your function in a single script, wrap it in the well known if __name__ == "__main__": construct to ensure it is only run if your code is used as a script.


* Well not strictly speaking, but it's close.

\$\endgroup\$
  • \$\begingroup\$ Could you also give feedback about the performance of this implementation? And, also, did you run the tests? \$\endgroup\$ – mathmaniage May 10 at 2:22
  • \$\begingroup\$ I mention your tests directly above the first code block. The first case took about 36s on average and was the base for all further testing. The second one didn't finish because my 16GB of RAM (or maybe rather 13GB of it) don't seem to be sufficient. \$\endgroup\$ – AlexV May 10 at 6:32
  • \$\begingroup\$ oh, sorry I didn't see the line, thank you for testing it. But it seems on average I bet the record of using set implementation given here: idlecoding.com/making-eratosthenes-go-faster-1 , and also with your idea of using jit, it surpasses it waaaay past. :-) \$\endgroup\$ – mathmaniage May 10 at 7:41
  • \$\begingroup\$ @jit doesn't work in class methods, right? There seems to be @jitclass but, that seems to be in development. \$\endgroup\$ – mathmaniage May 14 at 16:12
  • \$\begingroup\$ This not really suited to be answered in a comment, but this SO post might help. \$\endgroup\$ – AlexV May 14 at 17:37
5
\$\begingroup\$

At the end of their answer, @AJNeufeld has the right idea to improve the algorithm even further (in pure Python). You want to minimize the amount of numbers you have to mark off as being composite for every new prime you find. In order to do this, you can actually start at current_prime * current_prime, since you should already have marked up all smaller multiples of current_prime when you found all primes smaller than current_prime.

In addition I would use a slightly different data structure. Instead of having a list of all numbers up to limit, just have a list of True/False, indicating whether or not the number is (potentially) prime. You can immediately mark off 0 and 1 as not prime (by the usual definition) and then just need to proceed to the next True in the array, the index of which is the next prime number. Then mark off all multiples off that prime number as composite (False).

def prime_sieve(limit):
    prime = [True] * limit
    prime[0] = prime[1] = False

    for i, is_prime in enumerate(prime):
        if is_prime:
            yield i
            for n in range(i * i, limit, i):
                prime[n] = False

This function is a generator. You can just iterate over it to get one prime number at a time, or consume the whole generator using list(prime_sieve(1000)).

Using a for loop to mark off composites can be improved further using slicing:

prime[i*i:limit:i] = [False] * len(range(i*i, limit, i))

This uses the fact that the range object in Python 3 is a generator-like object which has a quick calculation of len implemented.

As for timings, here are a few. In the end, compiling it with numba is usually going to win, but you can get a factor of two just by making the algorithm a bit easier. And this factor of two carries through when jitting this answer as well.

enter image description here

\$\endgroup\$
  • \$\begingroup\$ thats a great comparison chart. Thanks! :-) \$\endgroup\$ – mathmaniage May 10 at 16:24
  • \$\begingroup\$ I've been thinking about the generator function above, even we do cut down the memory cost, It still builds a large prime[limit] list right? If we take a number say n in the order of $10^30$, this would mean infeasible memory capacity requirement, but without having to build a list, how can we generate primes? \$\endgroup\$ – mathmaniage May 16 at 16:03
  • 2
    \$\begingroup\$ @mathmaniage Yes, it would. Apparantly you can get the space requirements down to O(sqrt(n)) with algorithms more clever than this, though: stackoverflow.com/a/10733621/4042267 (and the other answers to that question). If you want O(1) space requirements, you are probably stuck with explicitly checking all relevant divisors for each number, which is quite a lot slower. \$\endgroup\$ – Graipher May 16 at 16:16
  • 1
    \$\begingroup\$ @mathmaniage Apparently yes, not sure how practical it is going to be, though: stackoverflow.com/questions/5602155/… The accepted answer won't directly work, though, since here we do need to modify the contents of the array. \$\endgroup\$ – Graipher May 16 at 16:36
  • 1
    \$\begingroup\$ I think your 3500000GB of memory is overstating the requirements by quite a bit. I've added memory usage analysis to my review. \$\endgroup\$ – AJNeufeld May 17 at 22:51
5
\$\begingroup\$

AlexV has made some great points. I’ll try not to duplicate them.

while numbers[current_prime] <= sqrt_limit:

This seems like odd code. numbers[x] == x is true, or numbers[x] == -1 is true. You should save the indirect lookup in the numbers array, and just loop:

while current_prime <= sqrt_limit:

You take great pains to explain in a comment that len(numbers) == n+1, but then you evaluate len(numbers)-1 anyway. It is just n. So you could simply write:

for i in range(current_prime, n // current_prime + 1):
    numbers[current_prime * i] = -1

Now, i starts at current_prime, and goes up to the last multiplier that results in current_prime * i still being within the bounds of the numbers list. You could simply use the range() method to count by multiples of current_prime, and let it figure out when it has gone beyond the limit for you:

for index in range(current_prime*current_prime, n+1, current_prime):
    numbers[index] = -1

Your summary loop at the end is not Pythonic:

for i in range(len(numbers)):
    # code which only uses numbers[i], never i on its own.

Instead of looping over the indices of the numbers list, you simply want to loop over the values of the list.

j = 0
for value in numbers:
    if value > 0:
        numbers[j] = value
        j += 1

But this loop is simply filtering values from one list into another list (although you are reusing the original list to store the new list). Python excels at filtering lists, using list comprehension. The above code can be written simply as:

numbers = [ value for value in numbers if value > 0 ]

Which you could then follow by:

return numbers[1:]

to trim off the first value of the list.


There are a number of things which you can do to further improve the performance. The easiest would be to special case “2” as prime, and then start at 3 and count by 2’s to skip over the even numbers. When clearing out multiples of prime_number, you would again skip the even multiples by incrementing the index by 2*prime_number, approximately doubling the speed of the algorithm by cutting the work in half.


Memory

Python is a great scripting language, but that ease of use comes with a cost. Everything is an object, and the Python environment manages the object allocation for you. This results in a lot of memory usage per object. With a few thousand objects, that isn't a problem. With hundreds of millions, that overhead starts to have a serious impact.

>>> sys.getsizeof(0)
24
>>> sys.getsizeof(1)
28
>>> sys.getsizeof(0x3fffffff)
28
>>> sys.getsizeof(0x40000000)
32

Positive integers between 1 and 1.07 billion take up 28 bytes each. If you have a list of the first billion integers, you're going to need 28 GB of memory for those integers. Plus ...

>>> sys.getsizeof([])
64
>>> sys.getsizeof([0])
72
>>> sys.getsizeof([0,1])
80
>>> sys.getsizeof([0,1,2])
88

... a list takes 64 bytes, plus 8 per object in the list. So the list of a billion entries is going to take an addition 8 GB of memory, for a total of 36 GB!

Instead of storing integer objects in a list, we can be way more efficient if we ask Python to store integers directly in an array.

>>> import array
>>> a = array.array('L', [0]*1000000000)
>>> sys.getsizeof(a)
4000000064

a is an array of 4-byte unsigned integers. And it only requires a mere 4 GB of memory.

As pointed out by @Graipher, you don't need to store the integers in the array for the sieve; you can simply store boolean flags. The array.array object is designed for numbers, not booleans, but we can store True and False as a 1 and a 0, in which case we only need 1 byte integers, so we can reduce the memory requirements down to 1 GB.

>>> a = array.array('B', [0, 1]*500000000)
>>> a[0:10]
array('B', [0, 1, 0, 1, 0, 1, 0, 1, 0, 1])
>>> sys.getsizeof(a)
1000000064

Instead of filling in the memory with a billion 0's, I've initialized all the odd elements to 1, for possibly prime, and 0 for the even elements, for not possibly prime. We just need a minor correction, so a[1] is flagged as not prime, and a[2] is prime.

array.array is overly complex for this issue; it has to do calculation based on the size of the elements. Since we are storing individual flags in bytes, a bytearray should be faster:

>>> a = bytearray(1000000000)
>>> a[1::2] = [1]*500000000         # Mark odd numbers as possible prime
>>> a[1:3] = [0, 1]                 # Mark 1 as not prime, 2 as prime
>>> list(a[0:20])
[0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1]
>>> sys.getsizeof(a)
1000000057

Now you can proceed with your sieve, testing only odd numbers, and marking as not prime only the odd multiples of a prime number, start with its square, up to the end of the sieve:

for i in range(3, int(sqrt(1000000000)), 2):
    if a[i]:
        a[i*i::2*i] = bytes(len(a[i*i::2*i]))

How many primes do you have? Just add up all of the 1 flags!

n = sum(a)

Where are you going to store these? How about in an array object, for efficiency? Since the prime number will go up to 1 billion, we'll again need 4-byte integer storage.

primes = array.array('L', (i for i, flag in enumerate(a) if flag))

import array, math, time

#n = 100
#n = 1000000
#n = 100000000
n = 1000000000

start = time.perf_counter()
a = bytearray(n)
a[1::2] = [1]*len(a[1::2])
a[1:3] = [0, 1]

for i in range(3, int(math.sqrt(n)), 2):
    if a[i]:
        a[i*i::2*i] = bytes(len(a[i*i::2*i]))

primes = array.array('L', (i for i, flag in enumerate(a) if flag))

end = time.perf_counter()
mins, secs = divmod(end - start, 60)

count = len(primes)
print(f"Found {count} primes below {n} in {int(mins)}m:{secs:06.3f}s")
if count < 30:
    print(primes)

Results:

Found 78498 primes below 1000000 in 0m:00.071s
Found 5761455 primes below 100000000 in 0m:06.623s
Found 50847534 primes below 1000000000 in 1m:09.585s

If our machine speeds are comparable, this sample point would fall between the eratos_sieve_jit (orange) and prime_sieve2_ (red) in @Graipher's plot.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.