Example code - Prime numbers optimization

This is going to be a long, but quite instructive rant about various optimization ideas and techniques.
Let us start with my "classic" prime number generator. The idea is to see how long time is takes to find/generate the primes between 0 and 100000000 (100 millions) and use that as a measure for performance. I need to point out that the basic speed of the computer plays a big role in getting these numbers.

We know that primes are integers which can only be divided by 1 and the primes itself. Therefore it is an obvious idea to evaluate if a number is a prime by starting to divide it by 2, 3, 4 ...n all the way up to the number itself and if it can be divided at any point, it is not a prime.
This basic idea can be improved in various ways.

• It is not necessary to divide with all the numbers up to n (the prospective prime). It is sufficient to divide with the numbers up to sqrt(n). This is because for any number where bignum * smallnum = n it also goes that smallnum * bignum = n, which means we only need to find the small number which is a divisor, to identify the number as not-prime. The small number will always be less or equal to sqrt(n), and the big number will always be greater or equal to sqrt(n). sqrt(n) is the equilibrium.
• All even numbers do not need to be checked. They are all not-primes, except 2. This means the sequence of numbers you test for division of n are odd, i.e. 3, 5, 7, 9, 11 ... sqrt(n).
• Actually, this can be improved as you only have to divide n by the primes in the testing, not all odd numbers. If you can divide n with an odd number, which is not a prime, then the odd number is a composite of at least two other smaller numbers, which you have tested for earlier.

Math definition: All non-primes are called composite numbers because they can be calculated as a multiplication of a number of primes.

In this PrimeGenerator class I am using the above tricks, remembering every prime as it is being calculated in order to efficiently calculate the following primes.
Time: 953 sec.

#!/usr/bin/env python3
# Prime number generator

class PrimeGenerator:
    # Class varible, known primes in consecutive order, can be extended, but must contain these
    knownprimes = [2, 3]
    # Highest tested number for prime
    highesttested = 3

    # Instatiation
    def __init__(self, number=None):
        if number is not None:
            if not isinstance(number, int):
                raise ValueError("Integer expected")
        self.target = number                
    
    # Initializing iteration
    def __iter__(self):
        if self.target is None:
            raise ValueError("No number specified")
        self.pos = 0
        return self
    
    # Find next prime
    def __next__(self):
        # Can we use the list of known primes to find the next?
        if self.pos < len(self.knownprimes):
            nextprime = self.knownprimes[self.pos]
            if nextprime >= self.target:
                raise StopIteration
            self.pos += 1
            return nextprime
        # No, start computing the next prime
        while self.target > PrimeGenerator.highesttested+1:
            PrimeGenerator.highesttested += 1
            if self._isprime(PrimeGenerator.highesttested):
                self.knownprimes.append(PrimeGenerator.highesttested)
                self.pos += 1
                return self.highesttested
        raise StopIteration

    # Private method for identifying a prime
    def _isprime(self, number):
        factor = 0
        pos = 0
        while factor*factor <= number:
            # find next potential factor either in known primes or odd numbers above last known prime
            if pos < len(self.knownprimes):
                factor = self.knownprimes[pos]
                pos += 1
            else:
                factor += 2
            # test if it truly is a factor
            if number % factor == 0:
                return False
        return True

    # It is nice be able to ask if an number is a prime
    def isprime(self, number=None):
        if number is None:
            number = self.target
        if not isinstance(number, int):
            raise ValueError("Integer expected")
        if number in PrimeGenerator.knownprimes:
            return True
        if number <= PrimeGenerator.highesttested:
            return False
        return self._isprime(number)

if __name__ == "__main__":
    for i in PrimeGenerator(100000000):
        print(i)
    # Big prime, don't try
    # 999296950101072104250052714631

The sieve of Eratosthanes

This old Greek figured out another method for computing the primes until "some number". You initialize a list n long with just True on every position. The True means 'this is a prime' for the number given by the position. Initially position 0 and 1 is set to False, as they are obviously not primes. Now you find the next True in the list, which would be a position 2. When a True is found, then go through the entire list in jumps of the number/position and the the value to False. The first time this happens is with 2, so position 4, 6, 8... will be set to False. These are the even numbers or more importantly the numbers that can be divided by 2 and are therefore not primes - they are composite numbers. The next position/number in the list is 3, so run through the list in jumps of 3 and set the value to False at the position, i.e. 6 (already False), 9, 12 (already False), 15 ...
Continue this pattern until the next True position is more than sqrt(n), at which point you can stop for the reasons already stated for the previous method.
The sieve is really about eliminating all composite numbers in the list more than finding primes. It is of course the same thing, but from different viewpoints. Part of what makes the sieve efficient is that we get rid of all the divisions, which take a lot of time.

Time: 18.6 sec.

def eratosthenes(size):
    size += 1  # Adjust for zero based list
    array = [True] * size
    array[0] = array[1] = False
    for i in range(int(size**0.5)):
        if array[i]:
            # A prime, now get rid of all composites where this prime is a divisor
            for j in range(i+i, size, i):
                array[j] = False
    return [i for i in range(size) if array[i]]

for i in eratosthenes(100000000):
    print(i)

This was a big performance boost. This is due to a simpler idea, which leads to simpler code.
What happens if we use a simpler bytearray instead of a list? The default init of bytearray is 0, so the logic must be reverted.

Time: 12.5 sec.

def eratosthenes2(size):
    size += 1  # Adjust for zero based list
    array = bytearray(size)
    array[0] = array[1] = 1
    for i in range(int(size**0.5)):
        if array[i] == 0:
            # A prime, now get rid of all composites where this prime is a divisor
            for j in range(i+i, size, i):
                array[j] = 1
    return [i for i in range(size) if array[i] == 0]

for i in eratosthenes2(100000000):
    print(i)

The return statement contains a comparison, which can be eliminated with an enumeration, or removed by a negation.

• return [i for i, is_prime in enumerate(array) if not is_prime]
• return [i for i in range(size) if not array[i]]

It comes out to approximately the same improvement.

Time: 12.2 sec.
Code not shown.

A common optimization is replacing

for j in range(i+i, size, i):

with

for j in range(i*i, size, i):

This is not extremely obvious. If, say 7, is the prime we found, then we start with eliminating 14, 21, 28, 35, 42, i.e. 7*n where 1 < n < 7. However, these numbers have already been eliminated earlier by other primes; 2, 3 & 5.

Time: 11.9 sec.

def eratosthenes3(size):
    size += 1  # Adjust for zero based list
    array = bytearray(size)
    array[0] = array[1] = 1
    for i in range(int(size**0.5)):
        if array[i] == 0:
            # A prime, now get rid of all composites where this prime is a divisor
            for j in range(i*i, size, i):
                array[j] = 1
    return [i for i in range(size) if not array[i]]

for i in eratosthenes3(100000000):
    print(i)

When we have lists or bytearrays, we can use slicing, or more specifically slice assignment. In other words, the inner loop can be replaced with a slice assignment. This trick works only because it is a very regular piece of the bytearray we replace. Good trick, nonetheless.

Time: 7.4 sec.

from math import ceil
def eratosthenes4(size):
    size += 1  # Adjust for zero based list
    array = bytearray(size)
    array[0] = array[1] = 1
    for i in range(int(size**0.5)):
        if array[i] == 0:
            # A prime, now get rid of all composites where this prime is a divisor
            array[i*i:size:i] = [1] * (((size-1) // i) - i +1)
    return [i for i in range(size) if not array[i]]