Welcome, guest | Sign In | My Account | Store | Cart

Fastest way to list all primes below N in python (Python recipe) by robert.william.hanks
ActiveState Code (http://code.activestate.com/recipes/577331/)

I was looking for a fast way to list all primes below n, so far i came up to this with the numpy solution the fastest. It does primes up to 10e6 in 15ms in my old machine, and it is capable of reaching 10e9.

Python, 33 lines
Download
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
    
#numpy solution
    
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=numpy.bool)
    for i in xrange(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)]

#pure-python solution

def primes1(n):
    """ 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 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 xrange(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 xrange(1,n/3-correction) if sieve[i]]

I am looking for some feedback,

  1. Does anyone know a faster way?
  2. Any recommendations in how to improve this code?

This code was first posted here: Fastest way to list all primes below N in python.
There you can find some benchmarks.

Tags: fast, fastest, math, numbers, numpy, primes, python

3 comments

Thomas Lehmann 13 years, 8 months ago  # | flag

primes1 does not return primes only!

  • The first output is by myself.
  • the second output is from your primes1 function -> primes1(100).
  • instead: primes2 is working fine.

[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]

[2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, 37, 41, 43, 47, 49, 53, 55, 59, 61, 65, 67, 71, 73, 77, 79, 83, 85, 89, 91, 95, 97]

robert.william.hanks (author) 13 years, 8 months ago  # | flag

@Thomas Lehmann: Strange, well i am using python 2.6.

To test your claim i just copy/paste the code above in ideone.com.

here is the link -> http://ideone.com/MyojM it contains code & output.

As you can see it is working fine there.

Lorenz Fischer 7 years, 6 months ago  # | flag 
pypy -m timeit -r10 -s"from sympy import sieve" "primes = list(sieve.primerange(1, 10**6))"
10 loops, best of 10: 2.03 msec per loop

So fast, so good. But if you are looking for the fastest isprime-function, you better use python with gmpy2.

python -m timeit -n10 -s"from gmpy2 import is_prime, mpz" "is_prime(mpz(2**521-1));is_prime(mpz(2**1279-1))"
10 loops, best of 3: 29.2 msec per loop
pypy -m timeit -n10 -s"from sympy import isprime" "isprime(2**521-1);isprime(2**1279-1)"
10 loops, best of 3: 214 msec per loop

I couldn't install gmpy2 for pypy, hence i had to compare python with gmpy2 against pypy with sympy.

Created by robert.william.hanks on Sat, 24 Jul 2010 (MIT)
◄ Python recipes (4591) ►
◄ robert.william.hanks's recipes (2) ►

Required Modules

  • (none specified)

Other Information and Tasks