Fastest way to list all primes below N

Question

This is the best algorithm I could come up.

def get_primes(n):
    numbers = set(range(n, 1, -1))
    primes = []
    while numbers:
        p = numbers.pop()
        primes.append(p)
        numbers.difference_update(set(range(p*2, n+1, p)))
    return primes

>>> timeit.Timer(stmt='get_primes.get_primes(1000000)', setup='import   get_primes').timeit(1)
1.1499958793645562

Can it be made even faster?

This code has a flaw: Since numbers is an unordered set, there is no guarantee that numbers.pop() will remove the lowest number from the set. Nevertheless, it works (at least for me) for some input numbers:

>>> sum(get_primes(2000000))
142913828922L
#That's the correct sum of all numbers below 2 million
>>> 529 in get_primes(1000)
False
>>> 529 in get_primes(530)
True

Solution 1

Warning: timeit results may vary due to differences in hardware or version of Python.

Below is a script which compares a number of implementations:

• ambi_sieve_plain,
• rwh_primes,
• rwh_primes1,
• rwh_primes2,
• sieveOfAtkin,
• sieveOfEratosthenes,
• sundaram3,
• sieve_wheel_30,
• ambi_sieve (requires numpy)
• primesfrom3to (requires numpy)
• primesfrom2to (requires numpy)

Many thanks to stephan for bringing sieve_wheel_30 to my attention. Credit goes to Robert William Hanks for primesfrom2to, primesfrom3to, rwh_primes, rwh_primes1, and rwh_primes2.

Of the plain Python methods tested, with psyco, for n=1000000, rwh_primes1 was the fastest tested.

+---------------------+-------+
| Method              | ms    |
+---------------------+-------+
| rwh_primes1         | 43.0  |
| sieveOfAtkin        | 46.4  |
| rwh_primes          | 57.4  |
| sieve_wheel_30      | 63.0  |
| rwh_primes2         | 67.8  |    
| sieveOfEratosthenes | 147.0 |
| ambi_sieve_plain    | 152.0 |
| sundaram3           | 194.0 |
+---------------------+-------+

Of the plain Python methods tested, without psyco, for n=1000000, rwh_primes2 was the fastest.

+---------------------+-------+
| Method              | ms    |
+---------------------+-------+
| rwh_primes2         | 68.1  |
| rwh_primes1         | 93.7  |
| rwh_primes          | 94.6  |
| sieve_wheel_30      | 97.4  |
| sieveOfEratosthenes | 178.0 |
| ambi_sieve_plain    | 286.0 |
| sieveOfAtkin        | 314.0 |
| sundaram3           | 416.0 |
+---------------------+-------+

Of all the methods tested, allowing numpy, for n=1000000, primesfrom2to was the fastest tested.

+---------------------+-------+
| Method              | ms    |
+---------------------+-------+
| primesfrom2to       | 15.9  |
| primesfrom3to       | 18.4  |
| ambi_sieve          | 29.3  |
+---------------------+-------+

Timings were measured using the command:

python -mtimeit -s"import primes" "primes.{method}(1000000)"

with {method} replaced by each of the method names.

primes.py:

#!/usr/bin/env python
import psyco; psyco.full()
from math import sqrt, ceil
import numpy as np

def rwh_primes(n):
    # https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    sieve = [True] * n
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i]:
            sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
    return [2] + [i for i in xrange(3,n,2) if sieve[i]]

def rwh_primes1(n):
    # https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns  a list of primes < n """
    sieve = [True] * (n/2)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i/2]:
            sieve[i*i/2::i] = [False] * ((n-i*i-1)/(2*i)+1)
    return [2] + [2*i+1 for i in xrange(1,n/2) if sieve[i]]

def rwh_primes2(n):
    # https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns a list of primes, 2 <= p < n """
    correction = (n%6>1)
    n = {0:n,1:n-1,2:n+4,3:n+3,4:n+2,5:n+1}[n%6]
    sieve = [True] * (n/3)
    sieve[0] = False
    for i in xrange(int(n**0.5)/3+1):
      if sieve[i]:
        k=3*i+1|1
        sieve[      ((k*k)/3)      ::2*k]=[False]*((n/6-(k*k)/6-1)/k+1)
        sieve[(k*k+4*k-2*k*(i&1))/3::2*k]=[False]*((n/6-(k*k+4*k-2*k*(i&1))/6-1)/k+1)
    return [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) if sieve[i]]

def sieve_wheel_30(N):
    # http://zerovolt.com/?p=88
    ''' Returns a list of primes <= N using wheel criterion 2*3*5 = 30

Copyright 2009 by zerovolt.com
This code is free for non-commercial purposes, in which case you can just leave this comment as a credit for my work.
If you need this code for commercial purposes, please contact me by sending an email to: info [at] zerovolt [dot] com.'''
    __smallp = ( 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
    61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
    149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
    229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311,
    313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401,
    409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491,
    499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
    601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683,
    691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
    809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887,
    907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997)

    wheel = (2, 3, 5)
    const = 30
    if N < 2:
        return []
    if N <= const:
        pos = 0
        while __smallp[pos] <= N:
            pos += 1
        return list(__smallp[:pos])
    # make the offsets list
    offsets = (7, 11, 13, 17, 19, 23, 29, 1)
    # prepare the list
    p = [2, 3, 5]
    dim = 2 + N // const
    tk1  = [True] * dim
    tk7  = [True] * dim
    tk11 = [True] * dim
    tk13 = [True] * dim
    tk17 = [True] * dim
    tk19 = [True] * dim
    tk23 = [True] * dim
    tk29 = [True] * dim
    tk1[0] = False
    # help dictionary d
    # d[a , b] = c  ==> if I want to find the smallest useful multiple of (30*pos)+a
    # on tkc, then I need the index given by the product of [(30*pos)+a][(30*pos)+b]
    # in general. If b < a, I need [(30*pos)+a][(30*(pos+1))+b]
    d = {}
    for x in offsets:
        for y in offsets:
            res = (x*y) % const
            if res in offsets:
                d[(x, res)] = y
    # another help dictionary: gives tkx calling tmptk[x]
    tmptk = {1:tk1, 7:tk7, 11:tk11, 13:tk13, 17:tk17, 19:tk19, 23:tk23, 29:tk29}
    pos, prime, lastadded, stop = 0, 0, 0, int(ceil(sqrt(N)))
    # inner functions definition
    def del_mult(tk, start, step):
        for k in xrange(start, len(tk), step):
            tk[k] = False
    # end of inner functions definition
    cpos = const * pos
    while prime < stop:
        # 30k + 7
        if tk7[pos]:
            prime = cpos + 7
            p.append(prime)
            lastadded = 7
            for off in offsets:
                tmp = d[(7, off)]
                start = (pos + prime) if off == 7 else (prime * (const * (pos + 1 if tmp < 7 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 11
        if tk11[pos]:
            prime = cpos + 11
            p.append(prime)
            lastadded = 11
            for off in offsets:
                tmp = d[(11, off)]
                start = (pos + prime) if off == 11 else (prime * (const * (pos + 1 if tmp < 11 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 13
        if tk13[pos]:
            prime = cpos + 13
            p.append(prime)
            lastadded = 13
            for off in offsets:
                tmp = d[(13, off)]
                start = (pos + prime) if off == 13 else (prime * (const * (pos + 1 if tmp < 13 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 17
        if tk17[pos]:
            prime = cpos + 17
            p.append(prime)
            lastadded = 17
            for off in offsets:
                tmp = d[(17, off)]
                start = (pos + prime) if off == 17 else (prime * (const * (pos + 1 if tmp < 17 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 19
        if tk19[pos]:
            prime = cpos + 19
            p.append(prime)
            lastadded = 19
            for off in offsets:
                tmp = d[(19, off)]
                start = (pos + prime) if off == 19 else (prime * (const * (pos + 1 if tmp < 19 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 23
        if tk23[pos]:
            prime = cpos + 23
            p.append(prime)
            lastadded = 23
            for off in offsets:
                tmp = d[(23, off)]
                start = (pos + prime) if off == 23 else (prime * (const * (pos + 1 if tmp < 23 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # 30k + 29
        if tk29[pos]:
            prime = cpos + 29
            p.append(prime)
            lastadded = 29
            for off in offsets:
                tmp = d[(29, off)]
                start = (pos + prime) if off == 29 else (prime * (const * (pos + 1 if tmp < 29 else 0) + tmp) )//const
                del_mult(tmptk[off], start, prime)
        # now we go back to top tk1, so we need to increase pos by 1
        pos += 1
        cpos = const * pos
        # 30k + 1
        if tk1[pos]:
            prime = cpos + 1
            p.append(prime)
            lastadded = 1
            for off in offsets:
                tmp = d[(1, off)]
                start = (pos + prime) if off == 1 else (prime * (const * pos + tmp) )//const
                del_mult(tmptk[off], start, prime)
    # time to add remaining primes
    # if lastadded == 1, remove last element and start adding them from tk1
    # this way we don't need an "if" within the last while
    if lastadded == 1:
        p.pop()
    # now complete for every other possible prime
    while pos < len(tk1):
        cpos = const * pos
        if tk1[pos]: p.append(cpos + 1)
        if tk7[pos]: p.append(cpos + 7)
        if tk11[pos]: p.append(cpos + 11)
        if tk13[pos]: p.append(cpos + 13)
        if tk17[pos]: p.append(cpos + 17)
        if tk19[pos]: p.append(cpos + 19)
        if tk23[pos]: p.append(cpos + 23)
        if tk29[pos]: p.append(cpos + 29)
        pos += 1
    # remove exceeding if present
    pos = len(p) - 1
    while p[pos] > N:
        pos -= 1
    if pos < len(p) - 1:
        del p[pos+1:]
    # return p list
    return p

def sieveOfEratosthenes(n):
    """sieveOfEratosthenes(n): return the list of the primes < n."""
    # Code from: <[email protected]>, Nov 30 2006
    # http://groups.google.com/group/comp.lang.python/msg/f1f10ced88c68c2d
    if n <= 2:
        return []
    sieve = range(3, n, 2)
    top = len(sieve)
    for si in sieve:
        if si:
            bottom = (si*si - 3) // 2
            if bottom >= top:
                break
            sieve[bottom::si] = [0] * -((bottom - top) // si)
    return [2] + [el for el in sieve if el]

def sieveOfAtkin(end):
    """sieveOfAtkin(end): return a list of all the prime numbers <end
    using the Sieve of Atkin."""
    # Code by Steve Krenzel, <[email protected]>, improved
    # Code: https://web.archive.org/web/20080324064651/http://krenzel.info/?p=83
    # Info: http://en.wikipedia.org/wiki/Sieve_of_Atkin
    assert end > 0
    lng = ((end-1) // 2)
    sieve = [False] * (lng + 1)

    x_max, x2, xd = int(sqrt((end-1)/4.0)), 0, 4
    for xd in xrange(4, 8*x_max + 2, 8):
        x2 += xd
        y_max = int(sqrt(end-x2))
        n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
        if not (n & 1):
            n -= n_diff
            n_diff -= 2
        for d in xrange((n_diff - 1) << 1, -1, -8):
            m = n % 12
            if m == 1 or m == 5:
                m = n >> 1
                sieve[m] = not sieve[m]
            n -= d

    x_max, x2, xd = int(sqrt((end-1) / 3.0)), 0, 3
    for xd in xrange(3, 6 * x_max + 2, 6):
        x2 += xd
        y_max = int(sqrt(end-x2))
        n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
        if not(n & 1):
            n -= n_diff
            n_diff -= 2
        for d in xrange((n_diff - 1) << 1, -1, -8):
            if n % 12 == 7:
                m = n >> 1
                sieve[m] = not sieve[m]
            n -= d

    x_max, y_min, x2, xd = int((2 + sqrt(4-8*(1-end)))/4), -1, 0, 3
    for x in xrange(1, x_max + 1):
        x2 += xd
        xd += 6
        if x2 >= end: y_min = (((int(ceil(sqrt(x2 - end))) - 1) << 1) - 2) << 1
        n, n_diff = ((x*x + x) << 1) - 1, (((x-1) << 1) - 2) << 1
        for d in xrange(n_diff, y_min, -8):
            if n % 12 == 11:
                m = n >> 1
                sieve[m] = not sieve[m]
            n += d

    primes = [2, 3]
    if end <= 3:
        return primes[:max(0,end-2)]

    for n in xrange(5 >> 1, (int(sqrt(end))+1) >> 1):
        if sieve[n]:
            primes.append((n << 1) + 1)
            aux = (n << 1) + 1
            aux *= aux
            for k in xrange(aux, end, 2 * aux):
                sieve[k >> 1] = False

    s  = int(sqrt(end)) + 1
    if s  % 2 == 0:
        s += 1
    primes.extend([i for i in xrange(s, end, 2) if sieve[i >> 1]])

    return primes

def ambi_sieve_plain(n):
    s = range(3, n, 2)
    for m in xrange(3, int(n**0.5)+1, 2): 
        if s[(m-3)/2]: 
            for t in xrange((m*m-3)/2,(n>>1)-1,m):
                s[t]=0
    return [2]+[t for t in s if t>0]

def sundaram3(max_n):
    # https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/2073279#2073279
    numbers = range(3, max_n+1, 2)
    half = (max_n)//2
    initial = 4

    for step in xrange(3, max_n+1, 2):
        for i in xrange(initial, half, step):
            numbers[i-1] = 0
        initial += 2*(step+1)

        if initial > half:
            return [2] + filter(None, numbers)

################################################################################
# Using Numpy:
def ambi_sieve(n):
    # http://tommih.blogspot.com/2009/04/fast-prime-number-generator.html
    s = np.arange(3, n, 2)
    for m in xrange(3, int(n ** 0.5)+1, 2): 
        if s[(m-3)/2]: 
            s[(m*m-3)/2::m]=0
    return np.r_[2, s[s>0]]

def primesfrom3to(n):
    # https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Returns a array of primes, p < n """
    assert n>=2
    sieve = np.ones(n/2, dtype=np.bool)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i/2]:
            sieve[i*i/2::i] = False
    return np.r_[2, 2*np.nonzero(sieve)[0][1::]+1]    

def primesfrom2to(n):
    # https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
    sieve[0] = False
    for i in xrange(int(n**0.5)/3+1):
        if sieve[i]:
            k=3*i+1|1
            sieve[      ((k*k)/3)      ::2*k] = False
            sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
    return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]

if __name__=='__main__':
    import itertools
    import sys

    def test(f1,f2,num):
        print('Testing {f1} and {f2} return same results'.format(
            f1=f1.func_name,
            f2=f2.func_name))
        if not all([a==b for a,b in itertools.izip_longest(f1(num),f2(num))]):
            sys.exit("Error: %s(%s) != %s(%s)"%(f1.func_name,num,f2.func_name,num))

    n=1000000
    test(sieveOfAtkin,sieveOfEratosthenes,n)
    test(sieveOfAtkin,ambi_sieve,n)
    test(sieveOfAtkin,ambi_sieve_plain,n) 
    test(sieveOfAtkin,sundaram3,n)
    test(sieveOfAtkin,sieve_wheel_30,n)
    test(sieveOfAtkin,primesfrom3to,n)
    test(sieveOfAtkin,primesfrom2to,n)
    test(sieveOfAtkin,rwh_primes,n)
    test(sieveOfAtkin,rwh_primes1,n)         
    test(sieveOfAtkin,rwh_primes2,n)

Running the script tests that all implementations give the same result.

Solution 2

Faster & more memory-wise pure Python code:

def primes(n):
    """ Returns  a list of primes < n """
    sieve = [True] * n
    for i in range(3,int(n**0.5)+1,2):
        if sieve[i]:
            sieve[i*i::2*i]=[False]*((n-i*i-1)//(2*i)+1)
    return [2] + [i for i in range(3,n,2) if sieve[i]]

or starting with half sieve

def primes1(n):
    """ Returns  a list of primes < n """
    sieve = [True] * (n//2)
    for i in range(3,int(n**0.5)+1,2):
        if sieve[i//2]:
            sieve[i*i//2::i] = [False] * ((n-i*i-1)//(2*i)+1)
    return [2] + [2*i+1 for i in range(1,n//2) if sieve[i]]

Faster & more memory-wise numpy code:

import numpy
def primesfrom3to(n):
    """ Returns a array of primes, 3 <= p < n """
    sieve = numpy.ones(n//2, dtype=bool)
    for i in range(3,int(n**0.5)+1,2):
        if sieve[i//2]:
            sieve[i*i//2::i] = False
    return 2*numpy.nonzero(sieve)[0][1::]+1

a faster variation starting with a third of a sieve:

import numpy
def primesfrom2to(n):
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    sieve = numpy.ones(n//3 + (n%6==2), dtype=bool)
    for i in range(1,int(n**0.5)//3+1):
        if sieve[i]:
            k=3*i+1|1
            sieve[       k*k//3     ::2*k] = False
            sieve[k*(k-2*(i&1)+4)//3::2*k] = False
    return numpy.r_[2,3,((3*numpy.nonzero(sieve)[0][1:]+1)|1)]

A (hard-to-code) pure-python version of the above code would be:

def primes2(n):
    """ Input n>=6, Returns a list of primes, 2 <= p < n """
    n, correction = n-n%6+6, 2-(n%6>1)
    sieve = [True] * (n//3)
    for i in range(1,int(n**0.5)//3+1):
      if sieve[i]:
        k=3*i+1|1
        sieve[      k*k//3      ::2*k] = [False] * ((n//6-k*k//6-1)//k+1)
        sieve[k*(k-2*(i&1)+4)//3::2*k] = [False] * ((n//6-k*(k-2*(i&1)+4)//6-1)//k+1)
    return [2,3] + [3*i+1|1 for i in range(1,n//3-correction) if sieve[i]]

Unfortunately pure-python don't adopt the simpler and faster numpy way of doing assignment, and calling len() inside the loop as in [False]*len(sieve[((k*k)//3)::2*k]) is too slow. So I had to improvise to correct input (& avoid more math) and do some extreme (& painful) math-magic.

Personally I think it is a shame that numpy (which is so widely used) is not part of Python standard library, and that the improvements in syntax and speed seem to be completely overlooked by Python developers.

Solution 3

There's a pretty neat sample from the Python Cookbook here -- the fastest version proposed on that URL is:

import itertools
def erat2( ):
    D = {  }
    yield 2
    for q in itertools.islice(itertools.count(3), 0, None, 2):
        p = D.pop(q, None)
        if p is None:
            D[q*q] = q
            yield q
        else:
            x = p + q
            while x in D or not (x&1):
                x += p
            D[x] = p

so that would give

def get_primes_erat(n):
  return list(itertools.takewhile(lambda p: p<n, erat2()))

Measuring at the shell prompt (as I prefer to do) with this code in pri.py, I observe:

$ python2.5 -mtimeit -s'import pri' 'pri.get_primes(1000000)'
10 loops, best of 3: 1.69 sec per loop
$ python2.5 -mtimeit -s'import pri' 'pri.get_primes_erat(1000000)'
10 loops, best of 3: 673 msec per loop

so it looks like the Cookbook solution is over twice as fast.

Solution 4

Using Sundaram's Sieve, I think I broke pure-Python's record:

def sundaram3(max_n):
    numbers = range(3, max_n+1, 2)
    half = (max_n)//2
    initial = 4

    for step in xrange(3, max_n+1, 2):
        for i in xrange(initial, half, step):
            numbers[i-1] = 0
        initial += 2*(step+1)

        if initial > half:
            return [2] + filter(None, numbers)

Comparasion:

C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.get_primes_erat(1000000)"
10 loops, best of 3: 710 msec per loop

C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.daniel_sieve_2(1000000)"
10 loops, best of 3: 435 msec per loop

C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.sundaram3(1000000)"
10 loops, best of 3: 327 msec per loop

Solution 5

The algorithm is fast, but it has a serious flaw:

>>> sorted(get_primes(530))
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73,
79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163,
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 527, 529]
>>> 17*31
527
>>> 23*23
529

You assume that numbers.pop() would return the smallest number in the set, but this is not guaranteed at all. Sets are unordered and pop() removes and returns an arbitrary element, so it cannot be used to select the next prime from the remaining numbers.

View more solutions

Comments

[Alex Martelli]

@jbochi, you're welcome -- but do look at that URL, including the credits: it took ten of us to collectively refine the code to this point, including Python-performance luminaries such as Tim Peters and Raymond Hettinger (I wrote the final text of the recipe since I edited the printed Cookbook, but in terms of coding my contribution was on a par with the others') -- in the end, it's really subtle and finely tuned code, and that's not surprising!-)

[jbochi]

@Alex: Knowing that your code is "only" twice as fast as mine, makes me pretty proud then. :) The URL was also very interesting to read. Thanks again.

[jbochi]

That's not pure Python, but it's the fastest version so far. thanks!

[Alex Martelli]

If you're interested in non-pure-Python code, then you should check out gmpy -- it has pretty good support for primes, via the next_prime method of its mpz type.

[Daniel G]

To be fair, though, you'd have to count the time downloading, unzipping, and formatting the primes and compare that with the time to generate primes using an algorithm - any one of these algorithms could easily write the results to a file for later use. I think in that case, given enough memory to actually calculate all primes less than 982,451,653, the numpy solution would still be faster.

[Kimvais]

@Daniel correct. However the store what you have and continue whenever needed still stands...

[Kimvais]

just for correctness, the code example should have import numpy as np

[Alok Singhal]

int(n ** 0.5) should be int(math.ceil(n ** 0.5)) or int(n ** 0.5) + 1. ambi_sieve(10) gives wrong answer otherwise.

[unutbu]

@Alok: Very true. Thanks for catching this.

[mthurlin]

I managed to speed up your function about 20% by adding "zero = 0" at the top of the function and then replacing the lambda in your filter with "zero.__sub__". Not the prettiest code in the world, but a bit faster :)

[jbochi]

@truppo: Thanks for your comment! I just realized that passing None instead of the original function works and it's even faster than zero.__sub__

[jbochi]

Thank's for the Miller-Rabin primality test, but this code is actually slower and is not providing the correct results. 37 is prime and does not pass the test.

[jbochi]

Or download them from here primes.utm.edu/lists/small/millions, but the idea is to test python's limit and see if beautiful code emerge from optimization.

[Ruggiero Spearman]

I guess 37 is one of the special cases, my bad. I was hopeful about the deterministic version though :)

[stephan]

You might want to add wheel factorization to your comparison. The version on zerovolt.com/?p=88 (which happened to be the first that I found googling) is faster on my machine than the other pure Python versions you list.

[unutbu]

@stephan: Thank you for the link. Are you sure you would not like to throw your hat into the ring?

[stephan]

@~unutbu: go ahead if you want. To me, your comparison of various very fast algorithms seems to be the perfect place to add this specific implementation. I have used your script to compare the speed and verify correctness.

[jfs]

@~unutbu: sieveOfEratosthenes() is faster if we allow only pure Python than any other variant (without psyco, numpy). Python3: sieveOfEratosthenes: 0.055, sieve_wheel_30: 0.061 ambi_sieve_plain: 0.130, sieveOfAtkin: 0.270; Python2: sieveOfEratosthenes: 0.061, sieve_wheel_30: 0.086, ambi_sieve_plain: 0.137, sieveOfAtkin: 0.216. If we allow C extensions then gmpy2-based variant (0.041) is the fastest in Python3 stackoverflow.com/questions/2897297/… ambi_sieve: 0.015 is the fastest on Python2 if C extensions are allowed.

[unutbu]

@J.F Sebastian: Thanks for the info. I've rechecked my timeit runs and at least on my machine (Intel E2160 running Ubuntu 9.10) sieveOfEratosthenes is definitely not the fastest. Do you think our differing results could be due to hardware? Also, have you checked Robert William Hanks' primesfrom3to1? stackoverflow.com/questions/2068372/… Using numpy I'm getting primesfrom3to1 (19.6) versus ambi_sieve (29.4).

[jfs]

@~unutbu: Indeed primesfrom3to1() (corrected to return r_[...) takes 17 seconds for n=1e9 vs. 22.2 seconds for ambi_sieve() (for n=1e6 both are almost the same (15ms) ambi_sieve() being slightly faster). Among pure Python versions sieveOfEratosthenes is the fastest on my machine (Intel i7, Ubuntu 10.04)

[jfs]

New @Robert William Hanks's primes() function stackoverflow.com/questions/2068372/… is the fastest (for n=1e6) now among pure Python solutions ideone.com/weu23

[unutbu]

@J.F. Sebastian: Neat program. I love the decorators. Using your program on my machine, primes1 is still sometimes faster than primes (well, on 2 out of 3 runs...). Without psycho the difference between those two seems to be on the order of the standard error between runs.

[Robert William Hanks]

in my machine times for primes & primes1 (without psycho) are exactly the same,i just posted as an example because it uses slightly less memory. i just posted a variation of primesfrom3to with the same upside that at least for 1e7 is faster in may machine.

[jfs]

@Robert William Hanks: primesfrom2to is slightly faster than ambi_sieve for n=1e6 ideone.com/9MOfq

[Robert William Hanks]

@J.F. Sebastian: sorry to ask you that, would you mind to test it again? i made a change in the translation part of the code primesfrom2to and i am getting better results in my machine.

[jfs]

@Robert William Hanks: new version of primesfrom2to is definitely faster than ambi_sieve for sufficiently large n. ideone.com/ddiU2

[Robert William Hanks]

@~unutbu: just added a new version of primesfrom2to & primes2, if you have time you may want to timeit those.

[jfs]

upto=10**6: primesfrom2to() - 7 ms; prime6() - 12 ms ideone.com/oDg2Y

[Ishbir]

And it can be made even faster with a minor change: see stackoverflow.com/questions/2211990/…

[Will Ness]

... And it can be made yet faster with additional ~1.2x-1.3x speedup, drastic reduction in memory footprint from O(n) to O(sqrt(n)) and improvement in empirical time complexity, by postponing the addition of primes to the dict until their square is seen in the input. Test it here.

[Will Ness]

for a postponed addition of primes into the dict (until the square of a prime is seen in the input) see stackoverflow.com/a/10733621/849891 .

[Adam]

Numpy is now compatible with Python 3. The fact that it's not in the standard library is good, that way they can have their own release cycle.

[Will Ness]

posting one timing for an unknown machine says nothing. The accepted answer here says "without psyco, for n=1000000, rwh_primes2 was the fastest". So if you'd provide your timings for that code as well as yours, on the same machine, and at 2, 4, 10 mln as well, then it'd be much more informative.

[GordonBGood]

-1, This code depends on special features of the bitarray implemented in C, which is why the code is fast as most of the work is being done in native code in the slice assignment. The bitarray package BREAKS the standard definition for proper slices (indexed over a range) for mutable sequences in that it allows assigning a single boolean 0/1 or True/False to all elements of the slice, whereas the standard behavior for pure Python seems to be to not allow this and only allow the assignment value of 0 in which case it is treated as a del of all of the slice elements from the sequence/array.

[GordonBGood]

cont'd: If calling non-standard native code were to be compared, we may as well write a "fastprimes" sequence generator package based on C code such as that of Kim Walisch's primesieve and generate all the primes in the four billion plus 32-bit number range in just a few seconds with a single call to the sequence generator. This would also use almost no memory as the linked code is based on a segmented Sieve of Eratosthenes and thus only uses a few ten's of Kilobytes of RAM, and if a sequence were generated there would be no list storage required.

[wrhall]

Did you know that if you pass sundaram3(9) it will return [2, 3, 5, 7, 9]? It seems to do this with numerous -- perhaps all -- odd numbers (even when they aren't prime)

[Misguided]

There isn't any special case for rabin miller.

[Scott Tesler]

Dead links are dead :(

[kmonsoor]

but, for noobies, which one is simplest, yet reasonably fast ?

[unutbu]

Conceptually, the simplest algorithm (to understand) is the Sieve of Eratosthenes. I would start there. If it is fast enough, great. If not, use the benchmarks posted above to guide you to a faster algorithm.

[Logan]

You misread the article. It is 31, not 37. This is why your implementation fails.

[Admin]

Using itertools.compress makes rwh_primes1 the fastest in the category "plain Python without Psycho" in my testing. I.e. replace the last line with return [2] + list(compress(xrange(3,n,2), sieve[1:]))

[Ehsan Kia]

If you're using pypy, these benchmarks (the psyco ones) seem fairly off. Surprisingly enough, I found sieveOfEratosthenes and ambi_sieve_plain to be the fastest with pypy. This is what I found for the non-numpy ones gist.github.com/5bf466bb1ee9e5726a52

[Ross]

@Daniel G I think download time is irrelevant. Isn't it really about generating the numbers, so you would want to take into account the algorithm used to create that list you're downloading. And any time complexity would ignore the once of file transfer given it O(n).

[jfs]

it is incorrect e.g., ajs_primes3a(10) -> array([2, 3, 5, 7, 9]). 9 is not a prime

[Alan James Salmoni]

You spotted an edge case I hadn't – well done! The problem was in 'for idx in range(3, int(upto ** 0.5), 2):' which should be 'for idx in range(3, int(upto ** 0.5) + 1, 2):'. Thanks but it works now.

[Alan James Salmoni]

The reason was that the idx loop went up to 'upto ** 05' which for cases up to and including 15. From 16 onwards, it works fine. This was a set of edge cases I hadn't tested for. Adding 1 means it should work for all numbers.

[jfs]

It seems to work now. It is the slowest among numpy-based solutions that return an array. Note: no true Sieve of Eratosthenes implementation uses modulo -- no need to mention it. You could use mat[idx*idx::idx] instead of mat[idx*2::idx]. And np.nonzero(mat)[0] instead of np.where(mat == True)[0].

[Alan James Salmoni]

Thanks JF. I tested against prime6() and got a result faster up to (IIRC) about 250k when prime6() took over. primesfrom2to() was faster. At up to 20m, ajs_primes3a() took 0.034744977951ms, prime6() took 0.0222899913788ms and primesfrom2to() took 0.0104751586914ms (same machine, same load, best of 10 timings). It's honestly better than I thought it would be!

[gaborous]

If someone wonders how the functions here fare against PG7.8 of Wikibooks for pure python without psyco nor pypy: for n = 1000000: PG7.8: 4.93 s per loop ; rwh_primes1: 69 ms per loop ; rwh_primes2: 57.1 ms per loop

[ljetibo]

It's not exactly what OP ordered but I fail to see why the downvote. It's a 2.8sec solution unlike some other outside modules. I've noticed in the source that it's threaded, got any tests on how well it scales?

[Colonel Panic]

@ljetibo cheers. The bottleneck seems to be copying C++ vector to Python list, thus the count_primes function is much faster than generate_primes

[Colonel Panic]

On my computer it can comfortably generate primes up to 1e8 (it gives MemoryError for 1e9) , and count primes up to 1e10. @HappyLeapSecond above compares algorithms for 1e6

[hiro protagonist]

to just store binary values in an array i suggest bitarray - as used here (for the simplest prime sieve; not a contender in the race here!) stackoverflow.com/questions/31120986/…

[355durch113]

When casting in the primesfrom2to() method, should the division be inside of the brackets?

[mndrake]

It would be faster if mat[idx*2::idx] = False was checked if needed. if mat[idx]: mat[idx*2::idx] = False

[Assem]

it has an issue: sundaram3(7071) includes 7071 while it is not prime

[nbro]

@unutbu Please, fix the links to the different algorithms.

[unutbu]

@nbro: Which links are broken?

[unutbu]

@nbro: The sieveOfAtkin like has been updated with a web.archive.org cached page.

[Will Ness]

to visualize your test results, plot them on log-log scale, to see and compare the empirical orders of growth.

[ABri]

@ Will thanks for the input, i'll have this in mind the next time i need such comparison

[Moebius]

For a pure python version compatible with python 3, follow this link : stackoverflow.com/a/33356284/2482582

[noɥʇʎԀʎzɐɹƆ]

Can you update this with PyPy, now that psyco is dead and PyPy has superseded it?

[PM 2Ring]

FWIW, the version in your 1st code block is the topic of this question. I should also mention that I used a Python 3 version of your code here.

[Batman]

The FAQ for the UTM prime page suggests calculating small primes is faster than reading them off a disk (the question is what small means).

[samerivertwice]

In Python 3 I used this function stackoverflow.com/a/3035188/7799269 but replaced / with // and xrange with range and they seemed much faster than this.

[hiro protagonist]

the zerovolt.com link is broken... but there is a wheel implementation in python right here: stackoverflow.com/a/30563958/4954037

[cs95]

Would be great if these functions and timings could be updated for python3.

[Him]

Holy buttsnacks this sucker is fast.

[Will Ness]

to compare algorithmic performances, it's better to plot at a log-log scale.

[lukejanicke]

I notice this prints the whole list, whereas from the community wiki answer primesfrom2to(10000) returns [ 2 3 5 ... 9949 9967 9973]. Is that shortening a NumPy nd.array thing?

[ack]

Link to example is dead.

[Alexander]

Thanks, I've removed the link.

[Will Ness]

mmm, log-log plots... :)

[Pierre D]

2.47ms for 1e6 using numba and a simple (older) Sieve of Atkin.

[Pierre D]

How about numba versions of these algos? For motivation, here I use numba to implement an older version of Atkin's sieve. It takes 2.5ms to get the primes up to 10**6 and 10.6s up to 10**9.

[ishandutta2007]

@PierreD that's as fast as my c++ implementation of Atkin's segmented sieve. But in c++ I currently use min_25 sieve which can sieve up to 10**9 in less than 1 sec. Can it be made avaiable in python?

[ishandutta2007]

In c++ I use min_25 sieve.I takes less than 1 sec to generate primes till 1 Billion. Which of the above catagories would it fall into for python ?