1.11M Members

Prime Number Function Improvements (Python)

 
0
 

This Python snippet shows you how to take typical beginner's example of prime number function and work on improving its speed. The decorator function for timing is ideal tool to follow your progress you make. The improvements are commented and you are encouraged to make further improvements. Python is fun!

# time a number of progressively imrove prime number functions
# a decorator timing function is used
# declare the @decorator just above the function to invoke print_timing()
# tested with Python25    HAB      03feb2007

import time

def print_timing(func):
    """set up a decorator function for timing"""
    def wrapper(*arg):
        t1 = time.clock()
        res = func(*arg)
        t2 = time.clock()
        print '%s took %0.3f ms' % (func.func_name, (t2-t1)*1000.0)
        return res
    return wrapper


@print_timing
def get_primes1(r=10):
    """
    unimproved original version
    """
    primes = []
    for i in range(2, r+1):
        prime = 1
        for divisor in range(2, i):
            if i % divisor == 0:
                prime = 0
                break
        if prime:
            primes.append(i)
    return primes

@print_timing
def get_primes2(r=10):
    """
    improved version, loop only over odd numbers starting with 3
    also limit the divisor range
    """
    primes = [2]
    # now ranges return odd numbers only ...
    for i in range(3, r+1, 2):
        prime = 1
        for divisor in range(3, i//2, 2):
            if i % divisor == 0:
                prime = 0
                break
        if prime:
            primes.append(i)
    return primes

@print_timing
def get_primes3(r=10):
    """
    improved version, loop only over odd numbers starting with 3
    and limit devisor range even more
    """
    primes = [2]
    # now ranges return odd numbers only ...
    for i in range(3, r+1, 2):
        prime = 1
        for divisor in range(3, i//7+3, 2):
            if i % divisor == 0:
                prime = 0
                break
        if prime:
            primes.append(i)
    return primes

@print_timing
def get_primes4(r=10):
    """
    improved version, loop only over odd numbers starting with 3
    and limit devisor range even more and using xrange()
    """
    primes = [2]
    # now ranges return odd numbers only ...
    for i in xrange(3, r+1, 2):
        prime = 1
        for divisor in xrange(3, i//7+3, 2):
            if i % divisor == 0:
                prime = 0
                break
        if prime:
            primes.append(i)
    return primes

@print_timing
def get_primes5(r=10):
    """
    improved version, loop only over odd numbers starting with 3
    and limit devisor range further with **0.5 and using xrange()
    """
    primes = [2]
    # now ranges return odd numbers only ...
    for i in xrange(3, r+1, 2):
        prime = 1
        for divisor in xrange(3, int(i ** 0.5)+1, 2):
            if i % divisor == 0:
                prime = 0
                break
        if prime:
            primes.append(i)
    return primes

@print_timing
def get_primes7(n):
    """
    standard optimized sieve algorithm to get a list of prime numbers
    --- this is the function to compare your functions against! ---
    """
    if n < 2:  return []
    if n == 2: return [2]
    # do only odd numbers starting at 3
    s = range(3, n+1, 2)
    # n**0.5 simpler than math.sqr(n)
    mroot = n ** 0.5
    half = len(s)
    i = 0
    m = 3
    while m <= mroot:
        if s[i]:
            j = (m*m-3)//2  # int div
            s[j] = 0
            while j < half:
                s[j] = 0
                j += m
        i = i+1
        m = 2*i+3
    return [2]+[x for x in s if x]

r = 9000
p1 = get_primes1(r)           
p2 = get_primes2(r)
p3 = get_primes3(r)
p4 = get_primes4(r)
p5 = get_primes5(r)
p7 = get_primes7(r)

# print the first 15 primes in the returned prime list for testing
print p1[:15]
print p2[:15]
print p3[:15]
print p4[:15]
print p5[:15]
print p7[:15]
print
# print the last 5 primes in the returned prime list for testing
print p1[-5:]
print p2[-5:]
print p3[-5:]
print p4[-5:]
print p5[-5:]
print p7[-5:]

"""
timing results for r=9000
get_primes1 took 2365.512 ms
get_primes2 took 467.231 ms
get_primes3 took 148.986 ms
get_primes4 took 104.222 ms
get_primes5 took 25.376 ms
get_primes7 took 2.483 ms
"""
 
0
 

Added sieve algorithm function to compare times with.

 
0
 

If you're going to generate a list of primes, rather than testing whether a single integer is prime, then why not use the list itself for the divisors?

 
0
 

If you're going to generate a list of primes, rather than testing whether a single integer is prime, then why not use the list itself for the divisors?

Ok here the more advanced alternatives:

from __future__ import print_function
import time
from math import sqrt, ceil
try:
    import numpy as np
except:
    print('No numpy')
    def primes_np(n):
        return []
    # no numpy asume Python 3, no xrange
    xrange = range
else:
    def primes_np(n):
        # http://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
        # for intensive use with non-numpy code we better transfer the format to standard int list
        return map(int,np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)])

try:
    import psyco
except:
    print('No psyco available')
else:
    psyco.full()

def rwh_primes1(n):
    # http://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 primes(n):
    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    # recurrence formula for length by amount1 and amount2 Tony Veijalainen 2010
    """ Returns  a list of primes < n """
    sieve = [True] * (n//2)
    amount1 = n-10
    amount2 = 6

    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i//2]:
             ## can you make recurrence formula for whole reciprocal?
            sieve[i*i//2::i] = [False] * (amount1//amount2+1)
        amount1-=4*i+4
        amount2+=4

    return [2] + [2*i+1 for i in xrange(1,n//2) if sieve[i]]

def get_primes(n):
    """
    standard optimized sieve algorithm to get a list of prime numbers
    --- this is the function to compare your functions against! ---
    """
    if n < 2:  return []
    if n == 2: return [2]
    # do only odd numbers starting at 3
    s = list(xrange(3, n+1, 2))
    # n**0.5 simpler than math.sqr(n)
    mroot = n ** 0.5
    half = len(s)
    i = 0
    m = 3
    while m <= mroot:
        if s[i]:
            j = (m*m-3)//2  # int div
            s[j] = 0
            while j < half:
                s[j] = 0
                j += m
        i = i+1
        m = 2*i+3
    return [2]+[x for x in s if x]

def sieve_of_Atkin(end):
    end += 1
    lng = (end-1)>>1   # originally: lng = (end/2)-1+(end&1)
    sieve = [False]*(lng + 1)
    x_max, x2 = int(sqrt((end-1)/4.0)), 0
    for xd in xrange(4, 8*x_max + 2, 8):
        x2 += xd
        y_max = int(sqrt(end-x2))
        n = x2 + y_max*y_max
        n_diff = (y_max << 1) - 1
        if n&1 == 0: 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 = int(sqrt((end-1)/3.0)), 0
    for xd in xrange(3, 6*x_max + 2, 6):
        x2 += xd
        y_max = int(sqrt(end-x2))
        n = x2 + y_max*y_max
        n_diff = (y_max << 1) - 1
        if n&1 == 0: 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 = ((x*x + x) << 1) - 1
        n_diff = (((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]; append= primes.append
    if end <= 3 : return primes[:max(0,end-2)]
    sqrtN  = int(sqrt(end)) + 1
    for n in xrange(2, sqrtN >> 1):
        if sieve[n]:
            i = (n << 1) +1;  j=i*i;  append(i)
            for k in xrange(j, end, 2*j):
                sieve[k>>1] = False
    if sqrtN&1 == 0: sqrtN += 1
    primes.extend([i for i in xrange(sqrtN, end, 2) if sieve[i >> 1]])
    return primes

def SoZP5(val):
    # all prime candidates > 5 are of form  30*k+(1,7,11,13,17,19,23,29)
    # initialize sieve array with only these candidate values
    # where sieve contains the odd integers representatives
    # convert integers to array indices/vals by  i = (n-3)>>1
    n1, n2, n3, n4, n5, n6, n7, n8 = -1, 2, 4, 5, 7, 8, 10, 13
    lndx = (val-1)>>1; sieve = [False]*(lndx+15)
    while n8 < lndx:
        n1 += 15;   n2 += 15;   n3 += 15;   n4 += 15
        n5 += 15;   n6 += 15;   n7 += 15;   n8 += 15
        sieve[n1] = sieve[n2] = sieve[n3] = sieve[n4] = \
        sieve[n5] = sieve[n6] = sieve[n7] = sieve[n8] = True
    # now initialize sieve with (odd) primes < 30
    sieve[0] = sieve[1] = sieve[2]  = sieve[4]  = sieve[5] = \
    sieve[7] = sieve[8] = sieve[10] = sieve[13] = True
    n = 0;  rescnt = 8
    for i in xrange(2, ((int(sqrt(val))-3)>>1)+1, 1):
        if not sieve[i]:  continue
        # p1=7*i, p2=11*i, p3=13*i --- p6=23*i, p7=29*i, p8=31*i,  k=30*i
        # j  i  7i  11i  13i  17i  19i  23i  29i  31i  30i
        # 7->2  23   37   44   58   65   79  100  107  105
        j = (i<<1)+3; j2 = j<<1; p1 = j2+j+i; p2 = p1+j2; p3 = p2+j; p4 = p3+j2
        p5 = p4+j; p6 = p5+j2; p7 = p6+j2+j;  p8 = p7+j;  k = p8-i
        x = k*(n>>3)
        n += 1  # x = k*(n/rescnt)
        p1 += x; p2 += x; p3 += x; p4 += x; p5 += x; p6 += x; p7 += x; p8 += x
        while p8 < lndx:
           sieve[p1] = sieve[p2] = sieve[p3] = sieve[p4] = \
           sieve[p5] = sieve[p6] = sieve[p7] = sieve[p8] = False
           p1 += k; p2 += k; p3 += k; p4 += k; p5 += k; p6 += k; p7 += k; p8 += k
        if p1 < lndx:
            sieve[p1] = False
            if p2 < lndx:
                 sieve[p2] = False
            if p3 < lndx:
               sieve[p3] = False
            if p4 < lndx:
                sieve[p4] = False
            if p5 < lndx:
                sieve[p5] = False
            if p6 < lndx:
                sieve[p6] = False
            if p7 < lndx: sieve[p7] = False
    if val < 3 : return [2]
    if val < 5 : return [2,3]
    primes = [(i<<1)+3 for i in xrange(2, lndx, 1)  if sieve[i]]
    primes[0:0] = [2,3,5]
    return primes

if __name__ == '__main__':
    numprimes=10**7
    for test in (primes_np, rwh_primes1, primes, get_primes, sieve_of_Atkin, SoZP5):
        t0= time.clock()
        pr = test(numprimes)
        t0 -= time.clock()
        print(test.__name__, 'sum', sum(pr), 'time', -t0*1000, 'ms')

""" no psyco timings (Python 2.7.1):
No psyco available
primes_np sum 3203324994356 time 848.42263472 ms
rwh_primes1 sum 3203324994356 time 1851.15929539 ms
primes sum 3203324994356 time 1842.69648796 ms
get_primes sum 3203324994356 time 4751.71900339 ms
sieve_of_Atkin sum 3203324994356 time 5066.25283381 ms
SoZP5 sum 3203324994356 time 2690.95541472 ms
"""
""" psyco timings (2.6.6):
primes_np sum 3203324994356 time 852.706698756 ms
rwh_primes1 sum 3203324994356 time 1554.33693388 ms
primes sum 3203324994356 time 1500.75470486 ms
get_primes sum 3203324994356 time 3224.53579994 ms
sieve_of_Atkin sum 3203324994356 time 978.697571898 ms
SoZP5 sum 3203324994356 time 677.250320921 ms
"""
Isn't it about time forums rewarded their contributors?

Earn rewards points for helping others. Gain kudos. Cash out. Get better answers yourself.

It's as simple as contributing editorial or replying to discussions labeled or OP Kudos

You
This is an OP Kudos discussion and contributors may be rewarded
Post:
Start New Discussion
Tags Related to this Article