goldbach.py

#!/usr/bin/env python
#
# Author:  Yotam Medini  [email protected] -- Created: 2006/April/23
#
# Find Goldbach sums
#

import sys

M = 2006  # check limit


# Eratosthenes - helping function.
# Mark by False all multipliers of n in the sieve.
# Note that sieve is a list that is being _changed_ here for the caller.
def mark_mults(sieve, n):
    mult = n*n
    while mult < len(sieve):
        sieve[mult] = False
        mult += n

def eratosthenes(M):
    primes = []
    sieve = (M + 1)*[True]
    n_primes = 0
    n = 2
    while n <= M:
        if sieve[n]:
            primes += [n]
            mark_mults(sieve, n)
        n += 1
    return primes


# Program begins
primes = eratosthenes(M)
n_primes = len(primes)
n = 4
while n <= M:
    found = False
    pindex1 = 0
    while (pindex1 < n_primes) and not found:
        pindex2 = pindex1
        while (pindex2 < n_primes) and (primes[pindex1] + primes[pindex2] < n):
            pindex2 += 1
        if (pindex2 < n_primes):
            p1 = primes[pindex1]
            p2 = primes[pindex2] 
            found = (p1 + p2 == n)
            if found:
                sys.stdout.write("%4d = %3d + %4d\n" % (n, p1, p2))
        pindex1 += 1
    if not found:
        sys.stderr.write("Goldbach is wrong! n=%d\n" % n)
        sys.exit(1)
    n += 2
sys.exit(0)