3
\$\begingroup\$

I wanted to test my performance optimizing skills, and so wanted to find how fast I could get the first \$n\$ prime numbers. I limited myself to only the standard library, as I'm sure numpy or another library written in C has a prime generator way faster than Python ever will be - and offloading to a library doesn't really improve my ability to improve performance.

I've implemented both the Sieve of Eratosthenes and the Sieve of Sundaram. The SoE was based off my answer here, and the SoS was based off Wikipedias definition. These are available at the end of the answer.

I improved performance by:

  • SoE: Vectorizing the creation of primes.

    primes[base*2:limit:base] = [False]*(ceil(limit / base) - 2)
    
  • SoE: Change the start of the slice from \$2b\$ to \$b^2\$. [1]

    primes[base*base:limit:base] = [False]*((((limit - base*base) - 1) // base) + 1)
    
  • SoE: Simplify the calculations - addition seems to be faster than multiplication.

    primes[base * base::base] = [False] * ((limit - 1) // base - base + 1)
    
  • SoE: Use itertools.compress, rather than a comprehension.

  • SoS: Vectorize the inner loop.

    start = 1 + 3*j
    step = 1 + 3*j
    primes[start::step] = [False] * ceil((n - start) / step)
    
  • SoS: Vectorize the creation of values that only have one value in the sequence.

    When \$\frac{n - \text{start}}{\text{stop}} = \frac{n - (1 + 3j)}{1 + 2j} \le 1\$ is equivalent to \$n \le 2 + 5j\$ we know we can stop at \$j = \frac{n - 2}{5}\$.

    multi_stop = (n - 2) // 5
    for j in range(1, multi_stop):
        start = 1 + 3*j
        step = 1 + 2*j
        primes[start::step] = [False] * ceil((n - start) / step)
    
    if multi_stop >= 1:
        single_start = multi_stop * 3 + 1
        primes[single_start::3] = [False] * ceil((n - single_start) / 3)
    
  • SoS: It doesn't seem like you need the if created above, and so you can just save wasted cycles.

I tried defining false = [False]*limit and slice it, but found it to be slower than creating new lists in the loop.

This got the following prime sieves:

from math import ceil
from itertools import compress


def sieve_eratosthenes(limit):
    if limit <= 1:
        return []

    primes = [True] * limit
    for base in range(2, int(limit**0.5 + 1)):
        if primes[base]:
            primes[base * base::base] = [False] * ((limit - 1) // base - base + 1)

    primes[0] = primes[1] = False
    return list(compress(range(limit), primes))


def sieve_sundaram(limit):
    if limit <= 1:
        return []

    n = (limit - 1) // 2
    primes = [True] * n
    for j in range(1, (n - 2) // 5):
        start = 1 + 3*j
        step = 1 + 2*j
        primes[start::step] = [False] * ceil((n - start) / step)

    return [2] + [2*i + 1 for i, p in enumerate(primes) if p][1:]

Both are faster than the both the original functions.

enter image description here enter image description here

Code to generate graphs:

from math import ceil
from itertools import compress

import numpy as np
import matplotlib.pyplot as plt
from graphtimer import Plotter, MultiTimer


def sieve_eratosthenes_orig(limit):
    if limit <= 1:
        return []

    primes = [True] * limit
    for base in range(2, int(limit**0.5 + 1)):
        if primes[base]:
            for composite in range(base * 2, limit, base):
                primes[composite] = False
    return [num for num, is_prime in enumerate(primes) if is_prime][2:]


def sieve_eratosthenes(limit):
    if limit <= 1:
        return []

    primes = [True] * limit
    for base in range(2, int(limit**0.5 + 1)):
        if primes[base]:
            primes[base * base::base] = [False] * ((limit - 1) // base - base + 1)

    primes[0] = primes[1] = False
    return list(compress(range(limit), primes))


def sieve_sundaram_orig(limit):
    if limit <= 1:
        return []

    n = (limit - 1) // 2
    primes = [True] * n
    for j in range(1, n):
        for i in range(1, j + 1):
            value = i + j + 2*i*j
            if value < n:
                primes[value] = False

    return [2] + [2*i + 1 for i, p in enumerate(primes) if p][1:]


def sieve_sundaram(limit):
    if limit <= 1:
        return []

    n = (limit - 1) // 2
    primes = [True] * n
    for j in range(1, (n - 2) // 5):
        start = 1 + 3*j
        step = 1 + 2*j
        primes[start::step] = [False] * ceil((n - start) / step)

    return [2] + [2*i + 1 for i, p in enumerate(primes) if p][1:]


def sieve_test(limit):
    if limit <= 1:
        return []

    n = (limit - 1) // 2
    primes = [True] * n
    multi_stop = (n - 2) // 5
    for j in range(1, multi_stop):
        start = 1 + 3*j
        step = 1 + 2*j
        primes[start::step] = [False] * ceil((n - start) / step)

    return [2] + [2*i + 1 for i, p in enumerate(primes) if p][1:]


def test():
    for exp in range(6):
        limit = 10 ** exp
        assert sieve_test(limit) == sieve_eratosthenes(limit)


def main():
    fig, axs = plt.subplots()
    axs.set_yscale('log')
    axs.set_xscale('log')
    (
        Plotter(MultiTimer([
            sieve_eratosthenes_orig,
            sieve_eratosthenes,
            sieve_sundaram,
            sieve_sundaram_orig,
            # sieve_test,
        ]))
            .repeat(5, 5, np.logspace(0.35, 2), args_conv=int)
            .min()
            .plot(axs, x_label='limit')
    )
    fig.show()


if __name__ == '__main__':
    test()
    main()

To use the above code snippet you need to install numpy, matplotlib and graphtimer. All should be available via pypi.

Can they be made faster, or is a different sieve faster?

\$\endgroup\$
  • \$\begingroup\$ Did you try the Sieve of Atkin? Something tells me it's faster than the Sieve of Erathosthenes. \$\endgroup\$ – Justin Jun 4 at 15:24
  • \$\begingroup\$ @Justin No, I didn't. I'm not too sure about it being faster than SoE in Python. But I'd be happy to be proven wrong. \$\endgroup\$ – Peilonrayz Jun 4 at 15:26
  • \$\begingroup\$ Here, I might have just found something - programmingpraxis.com/2010/02/19/sieve-of-atkin-improved. Does this help? \$\endgroup\$ – Justin Jun 4 at 15:30
  • \$\begingroup\$ @Justin Since it was written for Python 2 it may be unfair for it to be performance tested in Python 3. However with xrange = range and ensuring no IndexErrors I get the following graph. Which makes it faster than SoS at 1000 and always slower than SoE. \$\endgroup\$ – Peilonrayz Jun 4 at 15:39
  • \$\begingroup\$ Does this help - stackoverflow.com/questions/2068372/…? \$\endgroup\$ – Justin Jun 16 at 14:01

Your Answer

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

Browse other questions tagged or ask your own question.