15
\$\begingroup\$

I've implemented the sieve of Eratosthenes in Python 3 as follows:

def sieve_of_eratosthenes(n):
    is_prime = [True] * (n+1)
    is_prime[0] = False
    is_prime[1] = False
    p = 0

    while True:        
        for i, prime in enumerate(is_prime):
            if prime and i > p:
                p = i
                break
        else:
            break

        multiple = p + p

        while multiple <= n:
            is_prime[multiple]= False
            multiple = multiple + p

    r = []
    for i, prime in enumerate(is_prime):
        if prime: r.append(i)
    return r

Running this to 100,000 takes ~25 seconds.

I felt that was on the slow side, so I decided to take a different approach. The Sieve works prospectivly: it works ahead to clear out all multiples. I made a retrospective function, trial-dividing everything lower than the current prime-candidate:

def retrospect(n):
    primes = [2, 3]
    i = 5

    while i <= n:
        isprime = True
        lim = math.sqrt(n)+1
        for p in primes:
            if p > lim:
                break
            if i % p == 0:
                isprime = False
                break
        if isprime:
            primes.append(i)
        i += 2

    return primes

This is way faster!

Questions:

  • Are there any obvious shortcomings to the Sieve implementation?
  • Is there a point where the retrospect gets slower than the sieve?
\$\endgroup\$
27
\$\begingroup\$
    while True:        
        for i, prime in enumerate(is_prime):
            if prime and i > p:
                p = i
                break
        else:
            break

This is the cause of the slowness: track the state between loops and the 25 seconds drop to 0.05 seconds.


        multiple = p + p

        while multiple <= n:
            is_prime[multiple]= False
            multiple = multiple + p

It's probably more Pythonic to use range here:

        for multiple in range(p + p, n + 1, p):
            is_prime[multiple] = False

When handling larger inputs, it becomes worthwhile starting at p * p instead of p + p (on the basis that every composite multiple of p smaller than p * p also has a smaller prime factor).


    r = []
    for i, prime in enumerate(is_prime):
        if prime: r.append(i)
    return r

IMO it makes more sense to inline this into the main loop.


Edited code at this point:

def sieve_of_eratosthenes(n):
    is_prime = [True] * (n+1)
    is_prime[0] = False
    is_prime[1] = False

    primes = []

    for p in range(2, n+1):
        if is_prime[p]:
            primes.append(p)

            for multiple in range(p * p, n + 1, p):
                is_prime[multiple] = False

    return primes

Now if you want more speed at the cost of complicating things ever so slightly, you can apply a simple wheel by special-casing the only even prime:

def sieve_of_eratosthenes(n):
    is_prime = [True] * (n+1)
    primes = [2]

    for p in range(3, n + 1, 2):
        if is_prime[p]:
            primes.append(p)

            for multiple in range(p * p, n + 1, p + p):
                is_prime[multiple] = False

    return primes
\$\endgroup\$
  • 1
    \$\begingroup\$ Basically, your for loop only selects every number in range once and discards it if it isn't prime. The total number of iterations to find the next prime when n=1000 is 831, which when added to the number of primes < n is 999. In my original function, every number is tested again and again, leading to 77296 iterations. \$\endgroup\$ – steenbergh Jul 19 '18 at 11:04
8
\$\begingroup\$

Your prime sieve is overly complicated. You only need a single loop (and can also make it a generator to potentially save some memory):

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

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

This also uses the fact that p + p is already known to be composite for most numbers (since it is 2 * p and all multiples of 2 have already been marked composite once you have added 2 to the primes) and indeed all multiples of n * p for n < p have already been excluded. So it is sufficient to start at p * p, because that is the first multiple that has not been marked of as composite.

This is strictly faster than your retrospect function:

enter image description here

It is also linear with n, whereas retrospect is roughly \$\mathcal{O}(n^2)\$, due to having to iterate over all already found primes.


In your other function I would also use a for loop. It has the nice feature that it has an else clause that gets executed if the loop was not interrupted with a break statement:

def retrospect(n):
    primes = [2, 3]
    lim = math.sqrt(n) + 1

    for i in range(5, n + 1, 2):
        for p in primes:
            if p > lim:
                break
            if i % p == 0:
                break
        else:
            primes.append(i)
    return primes

It is also sufficient to calculate lim only once.

\$\endgroup\$
  • \$\begingroup\$ Couldn't you take your last 2 lines and combine them into a slice assignment? prime[p*p:limit:p] = False should speed it up even more. \$\endgroup\$ – Sunny Patel Jul 19 '18 at 21:17
  • \$\begingroup\$ @SunnyPatel That works only with advanced broadcasting rules like in numpy. Vanilla Python needs an iterable of correct length on the right side for that to work. \$\endgroup\$ – Graipher Jul 19 '18 at 21:20
  • \$\begingroup\$ Oooh, I thought I had implemented one like this before in the same fashion. \$\endgroup\$ – Sunny Patel Jul 19 '18 at 21:22
8
\$\begingroup\$

Your sieve function is rather convoluted and slow. If we think carefully about what's actually required, for each identified prime number, we cross off multiples. We can stop when we get to \$\sqrt{n}\$ and each prime cross-out operation for a prime \$i\$ can start with \$i^2\$.

def sieve_of_eratosthenes(n):
    is_prime = [True] * (n+1)
    is_prime[0] = False
    is_prime[1] = False

    limit = int(math.sqrt(n))+1
    for i in range(2, limit):
        if is_prime[i]:
            for j in range(i*i, n+1, i):
                is_prime[j] = False

    r = []
    for i, prime in enumerate(is_prime):
        if prime: r.append(i)
    return r

A still faster way to do it would be to note that all even numbers except for 2 are composite (that is, non-prime). We can slightly tweak the code to account for that:

def sieve_of_eratosthenes(n):
    is_prime = [False, True] * int(n/2 + 1)
    del is_prime[n+1:]
    is_prime[1] = False
    is_prime[2] = True

    limit = int(math.sqrt(n))+1
    for i in range(3, limit):
        if is_prime[i]:
            for j in range(i*i, n+1, 2*i):
                is_prime[j] = False

    r = []
    for i, prime in enumerate(is_prime):
        if prime:
            r.append(i)

    return r

Performance

Just for fun, I decided to repeat and expand the performance measurement that @Graipher made. Here's the code I used to do it:

import time
import matplotlib.pyplot as plt

def time_algo(algo, times):
    result = []
    for t in times:
        start = time.time()
        algo(t)
        result.append((time.time() - start)*1000)

    return result


names = [
         #'retrospect',
         taylor1,
         graipher,
         edward1,
         taylor2,
         edward2,
         ]
numarray = [1000, 3000,
            10000, 30000,
            100000, 300000,
            1000000, 3000000,
            10000000, 30000000,]
trials = []
for algo in names:
    trials.append(time_algo(algo, numarray))
for t in trials:
    plt.plot(numarray, t, '.-')
plt.ylabel('time (ms)') 
plt.xlabel('n')
plt.title("Prime algorithm perfomance")
plt.legend([alg.__name__ for alg in names])
plt.show()

Here are the results from the test on my machine: enter image description here

\$\endgroup\$
  • \$\begingroup\$ I've corrected the limit, as noted by @ack. \$\endgroup\$ – Edward Jul 19 '18 at 14:54
4
\$\begingroup\$

@Edward: corrected limit

def sieve_of_eratosthenes(n):
    is_prime = [False, True] * int(n/2 + 1)
    n_plus_one = n + 1  # save n+1 for later use
    del is_prime[n_plus_one:]  #
    is_prime[1] = False
    is_prime[2] = True

    #limit = int(math.sqrt(n))
    limit = int(n**.5) + 1  # import math  not required, limit corrected!
    for i in range(3, limit):
        if is_prime[i]:
            for j in range(i*i, n_plus_one, 2*i):  # limit corrected!
                is_prime[j] = False

    r = []
    for i, prime in enumerate(is_prime):
        if prime:
            r.append(i)

    return r

without corrected limit s_o_e(9) yields [2,3,5,7,9], so 9 would be prime (as would 11,13,17,...)

\$\endgroup\$
  • \$\begingroup\$ I don't think this constitutes an answer, it should probably be a comment on @Edward answer. Once you have enough reputation (50), you will be able to comment everywhere. \$\endgroup\$ – Graipher Jul 19 '18 at 7:36
  • \$\begingroup\$ @Graipher I agree, though this is a valid point; in fact I believe ack deserves rep for this find. It's a shame to see the user not rewarded for their effort because of a lack of reputation-points... \$\endgroup\$ – steenbergh Jul 19 '18 at 7:56
  • \$\begingroup\$ @steenbergh I agree that it is a good find. But if they had had more rep and commented this instead, they would not have gotten any rep from upvotes on that comment either. \$\endgroup\$ – Graipher Jul 19 '18 at 7:58
  • 1
    \$\begingroup\$ @Graipher I think this is a significant improvement enough to Edward's answer that is possible to stand as its own answer. (For this time) \$\endgroup\$ – Simon Forsberg Jul 19 '18 at 8:07
  • 1
    \$\begingroup\$ I suspect there is little difference among those options, but the way to find out is to measure. \$\endgroup\$ – Edward Jul 20 '18 at 23:08
2
\$\begingroup\$

On the performance side, using numpy arrays makes the computations 10 times faster of course.

Using Edward's code for performance, together with the following code ("Edward1" algorithm adapted to numpy) :

def sieve_of_eratosthenes(n):
    import numpy as np

    is_prime = np.ones(n,dtype=bool)
    is_prime[0:2] = False

    limit = int(math.sqrt(n))+1
    for i in range(2, limit):
        if is_prime[i]:
            is_prime[np.arange(i*i,n,i)] = False

    r = np.where(is_prime == True)[0]
    return r

I get on my machine : Prime numbers

\$\endgroup\$
  • \$\begingroup\$ Will definitely look into Numpy then! \$\endgroup\$ – steenbergh Jul 19 '18 at 13:39
  • 1
    \$\begingroup\$ This reproduces the bug in Edward's code. It's also not actually reviewing the code in the question, so by the standards of this site it shouldn't be an answer. \$\endgroup\$ – Peter Taylor Jul 19 '18 at 13:40
  • \$\begingroup\$ Sorry about that. This has just been corrected for. \$\endgroup\$ – D. Le Borgne Jul 19 '18 at 13:41

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.