Project Euler #10 with Python, Numpy sum summing incorrectly

Question

Project Euler #10

I know this question has been asked, but I am having trouble finding out why I am getting the wrong answer, and the other posts about this have not been helpful to me. The code is to find the sum of the prime numbers below 2,000,000.

import math
import numpy as np

def sum_primes(limit):      
    potential = 2
    primes = [potential]
    potential += 1
    primes.append(potential)

    while potential < limit:
        potential+=2
        test = True
        sqrt_potential = math.sqrt(potential)
        for a in primes:
            if a > sqrt_potential:
                break
            if potential%a == 0:
                test = False
                break
        if test and potential <= limit:
            primes.append(potential)

    print np.sum(primes)
    return

It appears that I am hitting all the primes but for some reason the sum isn't coming out correctly. Anyone see a problem here?

Here is some output I am getting:

>>> sum_primes(20000)
21171191
>>> sum_primes(200000)
1709600813
>>> sum_primes(2000000)
1179908154

Answer 1

I gave a fairly naive implementation on this, the last iteration took a while to return results.

def return_primes(upto=100):
    primes = []
    for i in range(2, upto+1):
        if not any(not i % p for p in primes):
            primes.append(i)
    return primes

And usage:

>>> sum(return_primes(upto=20000))
21171191
>>> sum(return_primes(upto=200000))
1709600813
>>> sum(return_primes(upto=2000000))
142913828922

Here's another implementation using the Sieve of Eratosthenes:

def return_primes(upto=100):
    primes = []
    sieve = set()
    for i in range(2, upto+1):
        if i not in sieve:
            primes.append(i)
            sieve.update(range(i, upto+1, i))
    return primes

This is much faster, running in about a second or two, compared to minutes for the above:

>>> sum(return_primes(200000))
1709600813
>>> sum(return_primes(2000000))
142913828922

---

To attempt to better diagnose your problem with this information, we see that

>>> 142913828922 % 2**32
1179908154

Which (thank you Will Ness) leads one to conclude that you're summing with 32bit integers. If we change your code to the following:

import math
import numpy as np

def sum_primes(limit):      
    potential = 2
    primes = [potential]
    potential += 1
    primes.append(potential)
    while potential < limit:
        potential+=2
        test = True
        sqrt_potential = math.sqrt(potential)
        for a in primes:
            if a > sqrt_potential:
                break
            if potential%a == 0:
                test = False
                break
        if test and potential <= limit:
            primes.append(potential)
    print np.sum(primes, dtype=np.int64)
    return

then you should get the correct output.

And I do replicate your behavior if I instead use:

    print np.sum(primes, dtype=np.int32)

Answer 2

In the initialization section of your script, you set potential = 3. Then the first statement in the while loop increments potential to 5. You never consider 3 as a potential prime. Fix it by setting potential -= 1 in the initialization section of your script.

A better algorithm uses the Sieve of Eratosthenes:

def sumPrimes(n):
    b, p, sum = [True] * (n+1), 2, 0
    for p in xrange(2, n+1):
        if b[p]:
            sum += p
            for i in xrange(p, n+1, p):
                b[i] = False
    return sum

If you're going to be doing the prime-number problems at Project Euler, you might want to read Programming with Prime Numbers at my blog.