The Goldbach distribution function

Goldbach's conjecture states that every even number is the sum of two primes. Computation suggests that this is true, and indeed that every even number seems to be the sum of two primes in many ways. The next obvious questions are "how many ways", and "how small can the first prime be" – it can't be arbitrarily small, since (n!+n-p) is not prime for any p between 1 and n, but it's certainly conceivable that the first prime can be less than 97 for more than 97% of even numbers.

So: define R(n) as the number of ways in which n can be written as the sum of two distinct primes (with order significant; for example R(10) = 3 since 10 = 3+7 = 5+5 = 7+3), and define p(n) as the smallest prime with n-p(n) also prime. It's not difficult to compute small tables of R and p; there are various increasingly-interesting ways of computing R quickly, and you may be surprised to learn that the most efficient one I can think of involves finite-field fast Fourier transforms.  If those scare you off, see here for direct approaches; if not, turn your attention here. 

This page is mostly about computing the R(n); see to see what R(n) looks like and what its properties include.

Computing R(n) by lists-of-primes methods

Finding R(n) for the first few thousand n is trivial: produce an array whose entries are one for primes and zero otherwise, then run i through the odd numbers between 1 and n, counting how often P[i] and P[n-i] are simultaneously one. This clearly takes time O(n), even if you store the array efficiently (one bit per odd number) and take advantage of symmetry by something like

if ((i%4 == 2) && !get_bit(prime_array, (i-2)/4)) r=1; // 2p 

for (unsigned int j=1;j<i/2;j+=2)
  if (!get_bit(prime_array,(j-1)/2) && !get_bit(prime_array, (i-j-1)/2)) 
    r+=2;
		

and so is infeasible for large N; it takes almost exactly 45 seconds for N=3x10⁵ on my newly-accelerated desktop PC, or four and a quarter hours for N=6x10⁶. This is the function goldbach_slowest in the source code.

Getting to O(N²(log N)^-1)

You'll notice that the inner loop above is iterating over primes by iterating over odd integers and checking which are prime. If instead you construct a table of (odd) primes, and replace the loop by

for (primi = prima.begin(); 2*(*primi)<i; primi++)
  if (!get_bit(prime_array, (i-(*primi)-1)/2)) 
    r+=2;

the inner loop will run p(i) times, and by the Prime Number Theorem the total run-time will be reduced by about a factor log N. On implementation, this takes nine seconds for N=3x10⁵. This is the functiongoldback_slower. But I didn't bother timing an overnight run, because ...

Getting to O(N²(log N)^-2)

If we consider the improved inner loop, we are basically iterating backwards over a list of integers by checking which are prime. We could instead try to use the table of primes; the cost is O(N) memory usage (for an array of N counters, at 32 bits each, rather than N is-it-prime flags at one bit each), but generally, when N² is a feasible amount of computation, N is a feasible amount of storage.

In pseudo-code, the loop looks like

• Set R[i] to zero for i = 0 .. N/2
• Let i run through odd primes less than or equal to N
  • If i <= N/2 then increment R[i]
  • Let j run through odd primes less than or equal to min(i,N-i)
    • Add 2 to R[(i+j)/2]

and in C++ it comes out as

vector<unsigned int> R(N/2+1);
for (i=0; i<np; i++)
{
  unsigned int pi = prime[i];
  if (pi+pi <= N) R[pi]++;
  int pj;
  for (int j=0; j<i && (pj=prime[j])<=(N-pi); j++)
  {
    R[(pi+pj)/2]+=2;
    if (first[(pi+pj)/2]==-1) first[(pi+pj)/2]=pj;
  }
}

This takes 1.6 seconds for N=3x10⁵ and just over 22 minutes for N=6x10⁶; we have an interestingly-shaped graph if we plot i against time, with the calculation running fairly slowly in the middle and then speeding up at the end as N-p_i gets small.

The overnight run mentioned above produced the quite stunningly irrelevant fact that the smallest number that can be written in more than half a million ways (counting order; a quarter-million, if you don't) as a sum of two primes is 23130030.

Calculating R values by convolution

Let P be the polynomial
x + x³ + x⁵ + x⁷+ x¹¹+ x¹³ + x¹⁷ + ...
whose exponents are the primes, and consider Q = P². The coefficient of x^q in Q will be sum(i= 0 .. q) (a_i a_q-i) where the a_i are the coefficients in P, and in particular are zero when i is not prime and one when it is. So the product is one if and only if i and q-i are both prime, in which case qis the sum of two primes. Since we sum over all i from zero to q, the coefficient of xⁿwill be R(n).