4
\$\begingroup\$

I wrote a simple prime number program that checks if a number is prime or not using sieve of eratosthenes, I also wrote an optimized version that only generates the sieve_arr when the val exceeds the length of sieve_arr.

sieve.py

# sieve.py

from array import array
from math import sqrt

def sieve(n=80):
    """Sieve of Erathosthenes"""

    sieve_arr = array('i', [True for _ in range(n+1)])
    sieve_arr_length = len(sieve_arr)
    
    for i in range(2, int(sqrt(sieve_arr_length)) + 1):
        if sieve_arr[i]:
            for j in range(i + i, n + 1, i):
                sieve_arr[j] = False;

    return sieve_arr


prime.py

# prime.py

import timeit
from array import array
from sieve import sieve

sieve_arr = array('i')

def isprime(val):
    global sieve_arr 
    if val >= len(sieve_arr):
        sieve_arr = sieve(val * 2)
    
    return sieve_arr[val]

def isprime2(val):
    return sieve(val)[val]

if __name__ == '__main__':
    TEST1 = '''
from prime import isprime
for i in range(1, 10001):
    isprime(i)
    '''

    TEST2 = '''
from prime import isprime2
for i in range(1, 10001):
    isprime2(i)
    '''

    tt1 = timeit.repeat(stmt=TEST2, repeat = 5, number=1)
    print('prime Time taken: ', min(tt1))
    tt2 = timeit.repeat(stmt=TEST1, repeat = 5, number=1)
    print('Optimized prime Time taken: ', min(tt2))

I ran the test and had the following results

prime Time taken:  9.446519022000302
Optimized prime Time taken:  0.002043106000201078

Note: I didn't check for 0 and 1 in isprime() and isprime2() which aren't prime number, I wrote the code with performance in mind.

\$\endgroup\$
7
  • 1
    \$\begingroup\$ Interesting. I didn't know that 0 and 1 are prime. \$\endgroup\$ – superb rain Nov 30 '20 at 17:38
  • \$\begingroup\$ They aren't, I was mostly writing it for performance. \$\endgroup\$ – theProgrammer Nov 30 '20 at 17:41
  • \$\begingroup\$ Why use array, instead of basic lists? \$\endgroup\$ – AMC Dec 1 '20 at 1:05
  • \$\begingroup\$ I don't understand the rationale behind this sieve_arr = sieve(val * 2), If you want to check val is prime or not, generate values until val why val*2? If you want to check if a number is prime or not, why use sieve of Eratosthenes in the first place? Check this out Primality test, 6k+1 optimization python code. \$\endgroup\$ – Ch3steR Dec 1 '20 at 5:27
  • 1
    \$\begingroup\$ Since Python 3.8 you can use math.isqrt instead of int(math.sqrt(...)) \$\endgroup\$ – Marc Dec 1 '20 at 7:41
3
\$\begingroup\$

Your sieve can be optimized further. Pull the case for i = 2 outside of the for loop and then change the for loop to only check odd indexes. Also use slice assignment instead of a nested loop. This code is 10x faster (<100ms to sieve primes to 1M compared to >1s for your sieve).

def sieve2a(n=80):
    """Sieve of Erathosthenes"""

    sieve_arr = [True]*n
    sieve_arr[0] = sieve_arr[1] = False
    
    sieve_arr[4::2] = [False]*len(range(4, len(sieve_arr), 2))
    for i in range(3, int(sqrt(len(sieve_arr))) + 1, 2):
        if sieve_arr[i]:
            sieve_arr[i*i::i] = [False]*len(range(i*i, len(sieve_arr), i))

    return sieve_arr
\$\endgroup\$
2
  • \$\begingroup\$ What's the purpose of sieve_arr[0] = sieve_arr[1] = False? Seems to be an attempt to correctly support n=0 and n=1, but if I try those, it of course crashes with an IndexError. \$\endgroup\$ – superb rain Dec 3 '20 at 14:22
  • \$\begingroup\$ @superbrain, Its just to keep the array in synch, so sieve_arr[n] is True if n is prime. \$\endgroup\$ – RootTwo Dec 4 '20 at 3:45
5
\$\begingroup\$

Built-in functions

As mentioned by Marc, Since Python 3.8 you can use math.isqrt. This is faster. It is also guaranteed to be accurate when the size of the argument exceeds the number of bits of mantissa in a float.

Type coercion

sieve_arr = array('i', [True for _ in range(n+1)])

Your array is declared to hold 'i' values: integers. But you are trying to storing boolean values in the array. This necessitates an unnecessary conversion from True and False to 1 and 0.

Underused variable

The value sieve_arr_length is used once. It effectively holds n + 1. In fact, sieve_arr_length = n + 1 might be a more efficient way to initialize it, rather than asking the array for its length. Then instead of recomputing n+1 each time the j loop starts, you could simply use sieve_arr_length.

Special case even numbers

Every even number greater than 2 is not prime. There is no need to check for these. You could simply start your i loop at 3, instead of 2, and go up by twos to check only the odd numbers greater than 1.

Starting the j loop at i+i, or i*2 is an optimization, but not enough of one. By the time you're checking 13, you've already eliminated 3*13, 5*13, 7*13, and 11*13 because you've eliminated multiples of smaller primes. You can simply start at i*i.

Plus, if you're only checking the odd prime candidates (because you've special-cased the even), you don't need to eliminate 14*13, 16*13, 18*13 and so on, because those are guaranteed to be even. Your step size can be 2*i, and eliminate half of the redundant assignments.

    # Odd prime candidate sieve ...
    for i in range(3, isqrt(sieve_arr_length) + 1, 2):
        if sieve_arr[i]:
            for j in range(i * i, sieve_arr_length, 2 * i):
                sieve_arr[j] = 0;

Memory

The array('i', ...) uses 4 or 8 (depending on architecture) bytes per element in the array. RootTwo's solution uses a list structure, instead of the array. This will increase the amount of storage required for the sieve, since a list must store a pointer in each list element.

Using a bytearray instead of an array('i', ...) would reduce your memory usage, down to 1 byte per sieve element, allowing you to maintain a larger sieve.

def ajneufeld_bytearray(n=80):
    """Sieve of Erathosthenes"""

    prime = bytearray((0, 1)) * ((n + 1) // 2)
    prime[1] = 0
    prime[2] = 1
    for i in range(3, isqrt(n) + 1, 2):
        if prime[i]:
            prime[i*i::2*i] = bytes(len(range(i*i, n, 2*i)))

    return prime

You can further pack 8 flags into each byte, dramatically decreasing the size of your sieve even further, but this bit packing would be complex to implement yourself. Fortunately it has already been done as a nice package: bitarray.

bitarray implementation

The bitarray package has nice bitslice assignment operators. RootTwo's slice assignment:

sieve_arr[i*i::i] = [False]*len(range(i*i, len(sieve_arr), i))

is much easier with bitarray. The equivalent statement would be:

sieve_arr[i*i::i] = False

assigning False to each bit in the slice.

from bitarray import bitarray
from math import isqrt

def ajneufeld(n=80):
    """Sieve of Erathosthenes"""

    prime = bitarray(n)
    prime.setall(False)
    prime[3::2] = True

    prime[2] = True
    for candidate in range(3, isqrt(n) + 1, 2):
        if prime[candidate]:
            prime[candidate*candidate::2*candidate] = False

    return prime

Test code:

import matplotlib.pyplot as plt
from timeit import timeit
from sys import getsizeof

if __name__ == '__main__':

    candidates = (theprogrammer, roottwo, ajneufeld_bytearray, ajneufeld_bitarray)
    limits = [10, 30, 100, 300, 1000, 3000, 10_000, 30_000, 100_000, 300_000]

    print("Memory usage:")
    for candidate in candidates:
        result = candidate(100_000)
        print(f"  {candidate.__name__:19} {getsizeof(result):6}")

    fig, ax = plt.subplots()

    for candidate in candidates:
        times = [ timeit(lambda: candidate(limit), number=1) for limit in limits ]
        ax.plot(limits, times, '-+', label=candidate.__name__)

    ax.legend()
    plt.xscale('log')
    plt.yscale('log')
    plt.show()

Output:


Memory usage:
  theprogrammer       400068
  roottwo             800056
  ajneufeld_bytearray 100057
  ajneufeld_bitarray   12564

enter image description here

\$\endgroup\$
9
  • \$\begingroup\$ Woow. This is a great improvement. \$\endgroup\$ – theProgrammer Dec 2 '20 at 22:20
  • \$\begingroup\$ IIRC, using bitarray would nearly double the access time compared to a regular list or bytearray. You may want to explicitly mention that time-memory trade-off when you introduce bitarray. \$\endgroup\$ – GZ0 Dec 3 '20 at 23:08
  • 2
    \$\begingroup\$ @GZ0 bytearray is only better that bitarray for sieves from 30_000 to 3_000_000. bitarray leads at the smaller sizes and larger sizes (on my laptop, Python 3.8 64-bit, bitarray 1.6.1) . Agreed, access to individual bits should be slower than to bytes, but L1 & L2 cache sizes probably start to play a role, and using 1/8th the memory gives bitarray an advantage there. \$\endgroup\$ – AJNeufeld Dec 3 '20 at 23:58
  • \$\begingroup\$ Correct me if am wrong, I don't think in a list, elements holds a pointer to the next element. A list uses array as its internal structure. \$\endgroup\$ – theProgrammer Dec 4 '20 at 13:42
  • \$\begingroup\$ @theProgrammer A Python list is an array of object references. The same object can be stored in many different variables or structures (like lists, dictionaries, tuples). Internally, this is accomplished via pointers. A list is implemented (CPython) as an array of object pointers. Not to be confused with a linked list, where each object has a pointer to the next. Try calling sys.getsizeof(x) for lists of different sizes. The result will depend only on len(x), regardless of the size or kind of each element in the list. \$\endgroup\$ – AJNeufeld Dec 4 '20 at 18:26

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.