An Empirical Look at the Goldbach Conjecture

The Goldbach conjecture says that every even number bigger than 2 is the sum of two primes. I imagine he tried out his idea on numbers up to a certain point and guessed that he could keep going. He lived in the 18th century, so he would have done all his calculation by hand. What might he have done if he could have written a Python program?

Let's start with a list of primes, say the first 100 primes. The 100th prime is p = 541. If an even number less than p is the sum of two primes, it's the sum of two primes less than p. So by looking at the sums of pairs of primes less than p, we'll know whether the Goldbach conjecture is true for numbers less than p. And while we're at it, we could keep track not just of whether a number is the sum of two primes, but also how many ways it is a sum of two primes.

    from sympy import prime
    from numpy import zeros

    N = 100
    p = prime(N)

    primes = [prime(i) for i in range(1, N+1)]
    sums = zeros(p, int)

    for i in range(N):
        # j >= i so we don't double count
        for j in range(i, N):
            s = primes[i] + primes[j]
            if s >= p:
                break
            sums[s] += 1

    # Take the even slots starting with 4
    evens = sums[4::2]

    print( min(evens), max(evens) )

This prints 1 and 32. The former means that every even number greater than 4 and less than p was hit at least once, that every number under consideration was the sum of two primes. The latter means that at least one number less than p can be written as a sum of two primes 32 different ways.

According to the Wikipedia article on the Goldbach conjecture, Nils Pipping manually verified the Goldbach conjecture for even numbers up to 100,000 in 1938, an amazing feat.

There are 9,952 primes less than 100,000 and so we would need to take N = 9592 in our program to reproduce Pipping's result. This took about seven minutes.

Update: As suggested in the comments, nearly all of the time is being spent generating the list of primes. When I changed the line

primes = [prime(i) for i in range(1, N+1)]

to

primes = [x for x in primerange(1, p)]

the runtime dropped from 7 minutes to 18 seconds.