Goldbach’s Conjecture

March 2, 2010

We present two solutions. The first solution generates a list of primes using the Sieve of Eratosthenes from a previous exercise, then for each prime checks if the other number is prime:

(define (goldbach1 n)
  (let loop ((ps (primes n)))
    (if (member (- n (car ps)) ps)
        (list (car ps) (- n (car ps)))
        (loop (cdr ps)))))

In most cases, one of the two primes is quite small. For instance, considering all the even primes to a million, only 484 have both primes greater than 200, and the maximum small-prime occurs at 503222 = 523 + 502699. We can exploit that fact to make a function that works much faster for large numbers:

(define (goldbach2 n)
  (let ((ps '(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61
         67 71 73 79 83 89 97 101 103 107 109 113 127 131 137
         139 149 151 157 163 167 173 179 181 191 193 197 199)))
    (let loop ((ps ps))
      (cond ((null? ps)
              (let loop ((p 211))
                (cond ((< n p) (error 'goldbach "conjecture fails"))
                      ((and (prime? p) (prime? (- n p))) (list p (- n p)))
                      (else (loop (+ p 2))))))
            ((prime? (- n (car ps)))
              (list (car ps) (- n (car ps))))
            (else (loop (cdr ps)))))))

Here we test each prime less than 200. If we still haven’t found the two primes, we start from 211 (the first prime larger than 200), test for primality using the Baillie-Wagstaff test from a previous exercise, test the other number for primality, and either report success or add two and loop. Goldbach2 is faster than goldbach1 because we don’t need to generate and store the list of primes, most of which are never used.

Here are some examples:

> (goldbach1 28))
(5 23)
> (goldbach2 28)
(5 23)
> (goldbach2 986332)
(353 985979)

You can run the program at http://programmingpraxis.codepad.org/B0jzdRps.

About these ads

Pages: 1 2

10 Responses to “Goldbach’s Conjecture”

  1. […] Praxis – Goldbach’s Conjecture By Remco Niemeijer In today’s Programming Praxis problem we have to test Golbach’s conjecture that every even number […]

  2. Remco Niemeijer said

    My Haskell solution (see http://bonsaicode.wordpress.com/2010/03/02/programming-praxis-goldbach%E2%80%99s-conjecture/ for a version with comments):

    import Data.Numbers.Primes
    
    goldbach :: Integer -> (Integer, Integer)
    goldbach n = head [(p, n - p) | p <- takeWhile (< n) primes, isPrime (n - p)]
    
  3. Jason said

    Here’s my quick and naive C implementation.

    #include <stdio.h>
    #include <stdlib.h>
    
    void show_usage();
    unsigned long *create_sieve_to_number(unsigned long number);
    
    int main (int argc, const char * argv[]) {
    
        if (argc != 2) {
    		show_usage();
    		return -1;
    	}
    	
    	
    	unsigned long number = atol(argv[1]);
    	
    	if ((number % 2) != 0) {
    		show_usage();
    		return -1;
    	}
    	
    	unsigned long *sieve = create_sieve_to_number(number);
    	
    	for (unsigned long i = 2; i < number; i++) {
    		if (sieve[i] == 1) {
    			for (unsigned long j = i; j < number; j++) {
    				if (sieve[j] == 1) {
    					if (i + j == number) {
    						printf("Solution found: %d + %d\n", i, j);
    						return 0;
    					}
    				}
    			}
    		}
    	}
    	
    	printf("no solution found!  pick up your Fields Medal!\n");
    	
    		
    	return 0;
    	
    }
    
    void show_usage(void) {
    	printf("usage: goldbach [even number]\n");
    }
    
    unsigned long *create_sieve_to_number(unsigned long number) {
    	unsigned long *sieve;
    	
    	sieve = (unsigned long *)malloc(sizeof(unsigned long) * (number + 1));
    	
    	for (int i = 0; i < number; i++) {
    		sieve[i] = 1;
    	}
    	
    	for (unsigned long i = 2; i < number; i++) {
    		if (sieve[i] == 1) {
    			for (unsigned long j = i * i; j < number; j = j + i) {
    				sieve[j] = 0;
    			}
    		}
    	}
    	
    	return sieve;
    }
    
  4. […] Today’s Praxis is on the Goldbach Conjecture which states that any even number greater than 2 can be expressed as the sum of two primes. The challenge is to write a program that will take in an even number, and spit out the two primes that can be added together to make it. […]

  5. Jason said

    D’oh! The prime flags of the sieve should obviously be bool, and not ulong. Way to waste memory there, Jason. :)

  6. Dave said

    Here’s a managed code solution. No sieve is used, because it seemed like that would be a waste in most cases. A more optimized IsPrime() function could be swapped in.

    public static class Goldbach
    {
    private static List primes = new List() { 2, 3 };

    public static void GetGoldbachPrimes(int value, out int prime1, out int prime2)
    {
    Debug.Assert(value > 2 && value % 2 == 0, “value > 2 && value % 2 == 0″);
    if (value <= 2 || value % 2 != 0)
    {
    throw new ArgumentException("value must be even number greater than 2", "value");
    }

    foreach (var prime in GetPrimes())
    {
    int difference = value – prime;
    if (IsPrime(difference))
    {
    prime1 = prime;
    prime2 = difference;
    return;
    }
    }

    throw new InvalidOperationException("Primes could not be found");
    }

    private static IEnumerable GetPrimes()
    {
    for (int i = 0; i < primes.Count; i++)
    {
    yield return primes[i];
    }

    for (int i = primes[primes.Count – 1] + 2; i < int.MaxValue; i += 2)
    {
    if (IsPrime(i))
    {
    primes.Add(i);
    yield return i;
    }
    }
    }

    private static bool IsPrime(int value)
    {
    if (value = 0)
    {
    return true;
    }

    int sqrt = Convert.ToInt32(Math.Ceiling(Math.Sqrt(value)));

    for (int i = 3; i <= sqrt; i += 2)
    {
    if (value % i == 0)
    {
    return false;
    }
    }

    return true;
    }
    }

  7. I just came across this site and had done some related work with GC, in terms of symmetric prime pair solutions
    eg for 24, we have (11/13), (7/17),(5/19) as solutions, where each pair is symmetric about N/2=12
    it uses Pari/GP built-in function is ispseudoprime

    mirrors(N)=
    {
    \\ N is the number of interest to find symmetric prime pairs
    \\ written as pari/GP script
    i=1;cnt=0;Q1=999;

    if(N%2!=0|N5! “);return(2) );

    print(“selected N = “,N);

    while(Q1>5,

    Q1=N-i;Q2=N+i;

    if(ispseudoprime(Q1)&&ispseudoprime(Q2),
    \\## print(” mirror pair at “,Q1,” / “,Q2); \\disable for speed or if N> BB1
    cnt++ ); \\ end IF

    i=i+2; \\ skip multiples if 5! for speed
    if(i%5==0,i=i+2) \\ BB1

    ); \\ end WHILE

    print(“# pairs found = “,cnt);
    print(“other candidate N: 6,12,30,60,180,210,360,420,1260,2310,2520,4620,”);

    return(0);
    }

  8. Mike said

    Two python versions.

    For many values of n, it is faster to generate all the primes less than n than it is to generate each prime (p) up to n/2 and test whether n-p is prime. So the second routine runs faster than the first.

    Reuses is_prime and primes_to from earlier problems.

    from primes import primes_to, is_prime
    
    def prime_pairs1( x ):
        return ((p,x-p) for p in primes_to(x/2) if is_prime( x-p ))
    
    
    def prime_pairs2( x ):
        prime = list( primes_to( x ) )
        
        lo, hi = 0, len(prime)-1
        while lo <= hi:
    	n = prime[hi] + prime[lo]
    	if n < x:
    		lo += 1
    	elif n == x:
    		yield prime[lo], prime[hi]
    		lo += 1
    	else:
    		hi -= 1
    
  9. Not as elegant but I think this will do for now.

    def goldbach(n):
    	prime = getPrime(n)
    	try:
    		return [(prime[i], prime[j]) for i in range(len(prime)) for j in range(i, len(prime)) if prime[i] + prime[j] == n][0]
    	except IndexError:
    		pass
    	
    #get prime less than n
    def getPrime(n):
    	#list odd numbers
    	odd = []
    	a = 3
    	while a <= n: 
    		odd.append(a)
    		a = a + 2
    
    	#list prime numbers via sieve of eratosthenes
    	for p in odd:
    		q = 2*p 
    		while p <= odd[-1]:
    			p = p + q
    			if p in odd: 
    				odd.remove(p)
    
    	odd.insert(0,2)
    	return odd
    

    Sample output:

    4 (2, 2)
    6 (3, 3)
    8 (3, 5)
    10 (3, 7)
    12 (5, 7)
    14 (3, 11)
    16 (3, 13)
    18 (5, 13)
    20 (3, 17)
    22 (3, 19)
    24 (5, 19)
    26 (3, 23)
    28 (5, 23)
    30 (7, 23)
    32 (3, 29)
    34 (3, 31)
    36 (5, 31)
    38 (7, 31)
    40 (3, 37)
    42 (5, 37)
    44 (3, 41)
    46 (3, 43)
    48 (5, 43)
    50 (3, 47)
    52 (5, 47)
    54 (7, 47)
    56 (3, 53)
    58 (5, 53)
    60 (7, 53)
    62 (3, 59)
    64 (3, 61)
    66 (5, 61)
    68 (7, 61)
    70 (3, 67)
    72 (5, 67)
    74 (3, 71)
    76 (3, 73)
    78 (5, 73)
    80 (7, 73)
    82 (3, 79)
    84 (5, 79)
    86 (3, 83)
    88 (5, 83)
    90 (7, 83)
    92 (3, 89)
    94 (5, 89)
    96 (7, 89)
    98 (19, 79)
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 645 other followers

%d bloggers like this: