# Goldbach Partition

I want to check the Goldbach conjecture for a big number of $n$, but I don't know how to define this in Mathematica.

These are my questions:

1. Find a pair of primes $(p,q)$ for every even integer $n$, such $p+q=n$, using Mathematica.

and

1. How to calculate number of ways to write an even number $n$ as the sum of two primes?

Take a look at IntegerPartitions, although it relies on brute-force enumeration that is unlikely to scale well.

f1[n_] := IntegerPartitions[n, {2}, Prime @ Range @ PrimePi @ n, 1]

f2[n_] := Length @ IntegerPartitions[n, {2}, Prime @ Range @ PrimePi @ n]


Test:

f1[3412]

{{3407, 5}}

f2[3412]

43


### Further experiments

Anton Antonov's recent answer surprised me, in that for larger n values his use of Select[FrobeniusSolve[{1, 1}, n], And @@ Map[PrimeQ, #] &] is faster than f2 above. (In version 10.1 under Windows.) It seems that long lists for the third parameter of IntegerPartitions is slow.

At the cost of increased memory consumption over f2 but less memory consumption than FrobeniusSolve we may enumerate with IntegerPartitions and then filter:

f3[n_] := Cases[IntegerPartitions[n, {2}], {__?PrimeQ}]


Timings on a fairly large integer:

f2[1787834]          // AbsoluteTiming

Length @ f3[1787834] // AbsoluteTiming

{7.73629, 6643}

{0.846121, 6643}


And Anton's method with the benefit of somewhat faster filtering (as used above):

Length @ Cases[FrobeniusSolve[{1, 1}, 1787834], {__?PrimeQ}] // AbsoluteTiming

{2.2008, 13285}


Note that in this output each solution is enumerated twice which I conjecture is the reason for the lower performance relative to f3.

I propose to use FrobeniusSolve. It seems it gives results fairly quickly. How large is the number $n$ ?

AbsoluteTiming[
res = Select[FrobeniusSolve[{1, 1}, 120022], And @@ Map[PrimeQ, #] &];
]
(* {0.340291, Null} *)

Row[{"Number of prime partitions: ", Length[res]}]
(* "Number of prime partitions: ", 1668 *)

Row[{"Sample: ", Take[SortBy[res, Abs[Subtract @@ #] &], 4]}]
(* "Sample: ", {{59981, 60041}, {60041, 59981}, {59921,
60101}, {60101, 59921}} *)


Here is another computation with a larger number:

In[9]:= AbsoluteTiming[
res = Select[FrobeniusSolve[{1, 1}, 7878344],
And @@ Map[PrimeQ, #] &];
]
Out[9]= {25.3882, Null}

• Interesting variation. I didn't expect this to be competitive with IntegerPartitions, but it seems that on larger numbers it is. Also I believe that your filter can be improved; try: AbsoluteTiming[ Cases[FrobeniusSolve[{1, 1}, 120022], {__?PrimeQ}]; ] – Mr.Wizard Dec 4 '16 at 11:52
• @Mr.Wizard Thanks for pointing out the speed-up with Cases / __?PrimeQ] ! I was not aware that it can be that significant. – Anton Antonov Dec 4 '16 at 22:27

As @Mr.Wizard showed, IntegerPartitions answers both your questions directly, and he warned that it will be slow for large $n$ because it calculates all possible partitions. There is a faster answer to your first question of finding just one partition of even $n=p+q$. Set $p\le q$, and note that usually $p$ is a small prime. The function GoldbachTest uses a While loop to check a sequence of small primes in the input list of candidates $p$. If the list is exhausted, the failure condition $\{0,0\}$ is returned.

GoldbachTest[n_?EvenQ, p_List] :=
Block[{m = Length[p], i = 1},
While[i <= m && CompositeQ[n - p[[i]]], i += 1];
If[i > m, {0, 0}, {#, n - #} &[p[[i]]]]
]


Given a list $p$ of the first $k$ primes, there is a smallest even $n$ which cannot be represented as $n=p+q$, with prime $q$. The sequence begins $\{6,12,30,98,98,98,98,220,308,...\}$, which is Sloane's A152522. This page links to a paper by Granville, Van de Lune, and te Riele, where they conjecture that the smallest prime $p$ in a Goldbach partition of even $n=p+q$ is $p<1.603*\log[n]^2 \log[\log[n]]$. They confirmed their conjecture for even $n \le 2*10^{10}$.

Thus, for an efficient test of Goldbach's conjecture from $n=n_{min}$ to $n=n_{max}$ try the following. Using ParallelTable will be even faster.

GoldbachTestList[nmin_?EvenQ, nmax_Integer] :=
With[{p = Prime[Range[PrimePi[1.603*Log[nmax]^2 Log[Log[nmax]]]]]},
Table[GoldbachTest[n, p], {n, nmin, nmax, 2}]
]


Timings show that GoldbachTestList is over 100 times faster than the equivalent IntegerPartitions formulation for $n_{max} \ge 10^5$.

• That's nice, +1. – ciao Dec 4 '16 at 4:20

I wanted to summarize solutions to this question based on IntegerPartitions, FrobeniusSolve, and a new hybrid method.

The first method f1[n] based on IntegerPartitions finds pairs of primes summing to n. @Mr.Wizard found this method to be slow for long lists in the third argument. Hence, he proposed function f3[n] to enumerate all pairs then filter.

f1[n_] := IntegerPartitions[n, {2}, Prime@Range@PrimePi@n]

f3[n_] := Cases[IntegerPartitions[n, {2}], {__?PrimeQ}]


For reasonably large n=1787834, function f3[n] is much faster.

With[{n = 1787834}, {AbsoluteTiming[Length@f1[n]], AbsoluteTiming[Length@f3[n]]}]
(* {{10.5924, 6643}, {0.809354, 6643}} *)


The first two methods, f4[n] and f5[n], based on FrobeniusSolve are from @Mr.Wizard and @AntonAntonov, respectively. The third method f6[n] adapts f5[n] by using Pick instead of Select.

f4[n_] := Cases[FrobeniusSolve[{1, 1}, n], {__?PrimeQ}]

f5[n_] := Select[FrobeniusSolve[{1, 1}, n], And @@ Map[PrimeQ, #] &]

f6[n_] := Pick[#, And @@@ PrimeQ[#]] &[FrobeniusSolve[{1, 1}, n]]


Pick is faster than Select, but the fastest formulation is f4[n] with Cases. The number of solutions is twice as large because both solutions {p,q} and {q,p} are returned. Since half of 1787834 is prime, there is also a solution {p,p} which is not duplicated.

With[{n = 1787834}, {
AbsoluteTiming[Length@f4[n]],
AbsoluteTiming[Length@f5[n]],
AbsoluteTiming[Length@f6[n]]}]
(* {{2.55104, 13285}, {5.35851, 13285}, {3.78575, 13285}} *)


The hybrid solution recognises that FrobeniusSolve[{1,1},n] simply returns pairs of integers {k,n-k} for k=0,1,...,n, and tests each pair for primality. Function f7[n] begins with a list of primes p about half as long as for f1[n], then tests for prime n-p.

f7[n_] := Block[{p = Prime[Range[PrimePi[Quotient[n, 2]]]]},
Transpose[{#, n - #} &[Pick[p, PrimeQ[n - p]]]]
]


Reversed solutions are not returned by f7[n], so the count is again 6643. The timing is about 10 times faster than the best solution so far, f3[n].

AbsoluteTiming[Length@f7[1787834]]
(* {0.062481, 6643} *)