Sum of primes under a certain number

Question

I have tried two different methods to find the sum of primes under a certain number. However both of these methods end up taking a very long time for large numbers. How can I optimise either of these two methods?

Method 1 (Sieve of Eratosthenes):

import time
import math

primes_in_range = [2]

n = int(input("Enter the range of your sieve"))

start = time.time()

for x in range(3, n+1, 2):
    primes_in_range.append(x)

d = 3

while d <= math.sqrt(n):
    if d in primes_in_range:
        for j in range(2 * d, n + 1, d):
            if j in primes_in_range:
                primes_in_range.remove(j)
    d += 2

print(sum(primes_in_range))
end = time.time()
print(end - start)

Method 2:

import time
import math

def is_prime(n):
  if n>1:
    for i in range (2, int(math.sqrt(n)) + 1):
      if (n % i) == 0:
        return False
    else:
      return True
  else:
    return False

n = int(input("Enter the range of your sieve"))

start = time.time()

total = 2

i = 3

while i < n:
    if is_prime(i):
        total += i
    i += 2

print(total)


end = time.time()
print(end -start)

Answer 1

Method 1 (Sieve of Eratosthenes) is definitely the faster method to use.

Starting point on my computer, summing primes less than 100,000 takes 39.916 seconds.

List to Set

The problem is basically this line:

primes_in_range = [2]

This creates a list, to which you add all odd numbers above 1, and then repeatedly search for existence of a number within (d in primes_in_range), then for all multiples of that number, search for the multiple (j in primes_in_range), and then search for it a second time in the process of removing it (primes_in_range.remove(j)). Searching a list is an \$O(n)\$ operation. When you do it approximately once for every number in the list, your algorithm devolves to \$O(n^2)\$.

On top of this, you’re computingmath.sqrt(n) many, many times - a relatively expensive operation - yet the value does not change. You only need to compute it once.

The second item you can solve by computing the \$\sqrt n\$ once, and store it in a local variable.

For the first issue, there are a few possibilities. One, you could use a set, instead of a list. A set has approximately \$O(1)\$ lookup time, which is much faster than \$O(n)\$.

primes_in_range = { 2 }

Since sets are not ordered, you don’t "append()" items to it, you "add" them instead.

Time has now dropped to 0.020511 seconds

Since the time has dropped to under a second, let's increase to primes less than 10,000,000 for better time measurements.

Time for summing primes less than 10,000,000: 3.519 seconds

Search to Remove

When removing, you are checking for the existence of the item in primes_in_range before removing it, since removing it when it is not present will raise an exception. The set class has a discard method, which will only remove an item if present. No need to check for it to be present ahead of time.

Time for summing primes less than 10,000,000: 3.453 seconds

Removing Multiples

When you find a prime, d, you cross off multiples of that prime, starting at 2*d, and going 3*d, 4*d, 5*d, and so on. But when you determined 3 was a prime number, you crossed off all multiples of 3, so crossing off 3*d is a waste of time, for d > 3. Similarly, crossing off multiples of the prime 5 means crossing off 5*d for d > 5 is also a waste of time. The first multiple of d that hasn't been crossed off will be d*d, so you can remove starting at that point:

        for j in range(d * d, n + 1, d):
            primes_in_range.discard(j)

Time for summing primes less than 10,000,000: 3.305 seconds

Similarly, since you are only considering odd prime candidates in your sieve, you don't need to cross off any even multiples. Starting at d*d, the next multiple you need to remove would be (d+2)*d, not (d+1)*d. Your step size can be 2*d.

        for j in range(d * d, n + 1, 2 * d):
            primes_in_range.discard(j)

Time for summing primes less than 10,000,000: 2.105 seconds

Loop Like a Native

Python has its own looping, which includes incrementing loop indexes. This native Python looping is usually faster than writing loops with their own index manipulation:

primes_in_range = {2}
primes_in_range.update(x for x in range(3, n+1, 2))

for d in range(3, math.isqrt(n) + 1, 2):
    if d in primes_in_range:
        primes_in_range.difference_update(j for j in range(d * d, n + 1, 2*d))

Here, constructing the set is done with a generator expression in primes_in_range.update(...), and removing prime multiples is done with a generator expression in primes_in_range.difference_update(...) which adds and removes entire swathes of items in a single operation:

Time for summing primes less than 10,000,000: 2.049 seconds

Direct Indexing

However, a set is still an expensive memory structure to maintain. It requires hashing the values which are being added, creating bins to store the values in, rebinning as the set size changes, and so on. An array of flags is much more efficient, time-wise:

primes_in_range = [False] * n
primes_in_range[2] = True
for x in range(3, n, 2):
    primes_in_range[x] = True

for d in range(3, math.isqrt(n) + 1, 2):
    if primes_in_range[d]:
        for j in range(d * d, n + 1, 2 * d):
            primes_in_range[j] = False

total = sum(idx for idx, flag in enumerate(primes_in_range) if flag)

The first line, above, allocates an array of n items, each containing the flag value False. Then, indexing into the array is an \$O(1)\$ operation. Set the flag to True for all odd numbers above 3, and then proceed to set the appropriate flags to False based on the Sieve of Eratosthenes algorithm.

Time for summing primes less than 10,000,000: 1.398 seconds

Bit Array

That last approach is a little inefficient when it comes to memory. It requires a list of n 28-byte items to store one True/False flag. These flags could be packed 8 to a byte, resulting in a 99.5% reduction in required memory! The bitarray package provides such an array of flags.

from math import isqrt
from bitarray import bitarray

def sum_of_primes_below(n):
    primes_in_range = bitarray(n)
    primes_in_range.setall(False)

    primes_in_range[2] = True      # 2 is a prime
    primes_in_range[3::2] = True   # Odd numbers starting at 3 are prime candidates

    for d in range(3, isqrt(n) + 1, 2):
        if primes_in_range[d]:
            primes_in_range[d*d::2*d] = False    # Reset multiples of prime candidate

    return sum(idx for idx, flag in enumerate(primes_in_range) if flag)

Time for summing primes less than 10,000,000: 0.455 seconds

Answer 2

If you're making a text interface, remember to put a prompt character with spaces like > so I know where to type. Even better would be to get an argument from sys.argv but I digress. Also, you should specify that specifically an integer greater than 1 is required, and raise an Exception otherwise.

After reading the pseudocode on Wikipedia your implementation is perfect. However your code is hard to read: include some comments. primes_in_range - in range of what? d - what is "d"? Since we have the ability to write variable names we should take advantage of it.

As for speed optimizations: in your first program you are storing the primes into a list. If you read the time complexity for Python collections, set is more efficient for deleting items - O(1) instead of O(n). After using a set instead of a list, I managed to significantly reduce the time it took for your code to run by more than 1000X.

In fact, after using a set instead of a list I was able to calculate primes up to 1,000,000 in 0.30 seconds, whereas using the list did not finish even after 20 minutes on my Core i7.

Your second solution has nothing horribly wrong with it, it's simply just a brute force algorithm. That's why it takes forever. Consider trying the Sieve of Atkin algorithm.

Answer 3

The problem with your sieve implementation is that the operations are not primitive operations but are rather time consuming.

You create an array of integers. Let's say n = 1 billion. You add 500 million integers. Each operation of adding an integer checks if there is enough space, increases the array size if needed, and adds a number. That’s slow but reasonably harmless.

Removing the non-primes is the killer. The harmless looking if j in primes_in_range visits all elements of the array until it finds the matching one. That can be tens of millions array elements checked for each multiple you are removing. And the remove(j) will find j again, remove it, and then fill the gap in the array by moving all the following items one position forward.

What you do instead is create an array of Boolean values, and set these Boolean to true or false, which will be a primitive and therefore fast operation.

An implementation in C or C++ would be expected to handle a few hundred million integers per second.

Measure the time; as long as it grows quadratically with n, you are still doing many non-primitive operations.

Answer 4

For your second algorithm which I suppose is slow, but not fatally slow like the first one, you can actually perform some brute force optimisations:

You are trying out whether n is divisible by 2, 3, 4, 5, 6, 7, 8, 9, 10, etc. Now if you think about it, n cannot fail the test of divisibility by 4, 6, 8 or 10, because it would have been divisible by 2 already. n also cannot fail the division test for 9, 15, 21, 27 etc. because it would have been divisible by 3 already.

If you check first whether n is divisible by 2, 3 and 5, then you only need to check divisibility by other numbers that are not divisible by 2, 3 or 5. Only 8 out of every 30 consecutive numbers are not divisible by 2, 3, or 5, and that is the numbers 30k+1, 30k+7, 30k+11, 30k+13, 30k+17, 30k+19, 30k+23 and 30k+29. If you only check divisibility by these numbers, then only 8 out of 30 tests are needed. That will be about 4 times faster.

You can do even better by only checking divisibility by primes. Store all the small primes in an array, and only check divisibility by primes. If n > p^2 for the largest stored prime p, you have to find the next prime first.

And please make sure that your algorithm gives the correct result for n = 1, n = 0, n < 0.