# Question about how to speed up Mathematica code

When looking at the Minimal Goldbach prime partition point {p,q} for each n; where n=10^i and i = 2,3,4,...,10; I notice that these points reside in an interval with center n/2 and radius 250:

Clear[mgppp, lowerbound, upperbound]

mgppp[n_?EvenQ] :=
Last[ Transpose@{n - #, #} &@
Select[n - Prime@Range@PrimePi[n/2], PrimeQ]]

lowerbound[n_?EvenQ] := Module[{d}, {m = n/2, d = (m - 250)}; d]

upperbound[n_?EvenQ] := Module[{d}, {m = n/2, d = (m + 250)}; d]

hello = Table[{n, lowerbound[n], mgppp[n][[1]], n/2, mgppp[n][[2]],
upperbound[n]}, {n, {10^2, 10^3, 10^4, 10^5, 10^6, 10^7, 10^8,
10^9, 10^10}}];TableForm[hello,
TableHeadings -> {None, {"n", "(n/2)-250", "MGPPPp", "n/2", "MGPPPq",
"(n/2)+250"}}, TableAlignments -> Center]


Mathematica takes about 45 minutes to build the above table. Notice that each lower bound value (n/2)-250 is less than p and each upper bound value (n/2)+250 is greater than q.

Now I want to verify (or not) that the same happens for ALL even n's between 10^2 and 10^10. To simplify the number of evaluations; I compare only the lowerbound values with the p values as follows:

Clear[mgppp, lowerbound]

mgppp[n_?EvenQ] :=
Last[ Transpose@{n - #, #} &@
Select[n - Prime@Range@PrimePi[n/2], PrimeQ]]

lowerbound[n_?EvenQ] := Module[{d}, {m = n/2, d = (m - 250)}; d]

bag = {}; For[n = 10^2, n <= 10^10, n += 2,
If[lowerbound[n] < mgppp[n][[1]], Nothing,
AppendTo[bag, {n, lowerbound[n], mgppp[n][[1]]}]]]; TableForm[bag,
TableHeadings -> {None, {"n", "(n/2)-250", "MGPPPp"}},
TableAlignments -> Center]


The above evaluation is very slow and it takes days to be evaluated using my personal computer. I would appreciate any suggestions about how to speed up the calculations; perhaps using "Compile" or any other method to do this. Thank you.

• One of the problems, perhaps the problem, is PrimePi. Your algorithm as it is written will never be faster than PrimePi[Range[1, 10^10/2, 2]]; // AbsoluteTiming. Feb 23, 2018 at 20:27
• You are right C.E.; thank you; that version of mpppp is slow. The following version is much faster: mgppp2[n_?EvenQ] /; n > 3 := Block[{m = PrimePi[n/2], p}, While[! PrimeQ[q = n - (p = Prime[m])], m--]; {p, q}] But, I still have the same evaluation issues after using this second version. Feb 23, 2018 at 21:28
• What I am concerned about is actually the For loop. It will call mgpppp or mgppp2 once for every n between 10^2 and 10^8, so it will not be faster than PrimePi[Range[10^2, 10^10, 2]/2]; because those functions evaluate PrimePi for n divided by two. PrimePi appears to be slow for large n and that is a problem. Feb 24, 2018 at 7:52
• n is even. m=n/2 is an integer and PrimePi[m] is then calculated. These PrimePi[m] values fall in an interval with center n/2 and radius 250; so, computationally this shouldn’t be time exhausting. The problem lies someplace else C.E. Feb 25, 2018 at 15:22
• You're right, I made a mistake. I posted a more careful analysis. Feb 25, 2018 at 21:54

The combination of For and AppendTo in your code is an inefficient way of substituting Table, and although it is not a major culprit we might as well change that. What's more serious is that if mgppp[n][[1]] is larger than the lower bound then you evaluate it a second time. This is an expensive function, so it shouldn't be evaluated unless it has to be evaluated, so that should also be fixed. We end up with this code:

Clear[mgppp, lowerbound]

mgppp[n_?EvenQ] := Last[Transpose@{n - #, #} &@Select[n - Prime@Range@PrimePi[n/2], PrimeQ]]
lowerbound[n_?EvenQ] := Module[{d}, {m = n/2, d = (m - 250)}; d]

RuntimeToolsProfile[
res = Table[With[{y = First@mgppp[n]},
If[lowerbound[n] < y, Nothing, {n, lowerbound[n], y}]
], {n, 10^2, 10^5, 2}];
];


Where I have changed the upper limit from 10^10 to 10^5 and have wrapped the whole thing with a function that does profiling. You can read about profiling Mathematica code in the Q&A Profiling from Mathematica.

This code will not only run the computation but will also report (in debug mode – this is explained in the other question) on how much time each part of the code took. I changed the upper limit so that it would finish in a reasonable time. It outputs a report like this:

These reports are usually straightforward to interpret, in this case, however, I'm not sure. It says that Prime@Range@PrimePi[n/2] took 254.6 seconds, whereas Prime took 88 seconds. The difference should be in the time it took to evaluate PrimePi, but it says that PrimePi only took 0.46 seconds to evaluate. I'll ask about this in the chat room to see if anyone understands this.

In either case, I believe that the time is spent evaluating built-in functions. There isn't much overhead that can be optimized. Memoization can be used if we suspect that the code is doing the same work several times over, so I tried doing this:

Clear[prime, primeQ]

prime[x_] := prime[x] = Prime[x]
primeQ[x_] := primeQ[x] = PrimeQ[x]

RuntimeToolsProfile[
res = Table[With[{y = First@mgppp[n]},
If[lowerbound[n] < y, Nothing, {n, lowerbound[n], y}]
], {n, 10^2, 10^5, 2}];
];


But it didn't make it faster. So at this point, I'm wondering if this can really be optimized at the code level.

• I think Prime and PrimeQ cache results. My guess is that the 135,000,000 calls to Prime didn't add to the "time" of Prime when retrieving a cached result, but they still took around 135000000 * 10^-6 sec. to execute the function calls. Feb 25, 2018 at 22:06
• @MichaelE2 That makes sense, thank you. Feb 26, 2018 at 4:49
• I apologize for not not returning fast and proper feedback to all of your suggestions above but, my personal computer was attacked by a virus today at work and our IT technicians had to clean my PC and I had to re-install Mathematica and get use to a new Windows 10 Operations System. Feb 26, 2018 at 21:11
• Clear[mgppp, lowerbound, upperbound] mgppp[n_?EvenQ] /; n > 3 := Block[{m = PrimePi[n/2], p}, While[! PrimeQ[q = n - (p = Prime[m])], m--]; {p, q}] lowerbound[n_?EvenQ] := Module[{d}, {m = n/2, d = (m - 250)}; d] (** 10^2 to 10^3 range: **) hello = Table[{n, lowerbound[n], mgppp[n][[1]], n/2, mgppp[n][[2]], upperbound[n], If[lowerbound[n] < mgppp[n][[1]], "ok", "Not bounded"]}, {n, 10^2, 10^3, 2}]; TableForm[hello, TableHeadings -> {None, {"n", "(n/2)-250", "MGPPPp", "n/2", "MGPPPq", "(n/2)+250", "Bounding Status"}}, TableAlignments -> Center] Feb 26, 2018 at 21:12

The main problem is you are generating a large number of prime pairs then using Last to throw most of them away. This starts in the middle and stops as soon as we find a pair: This runs in a fraction of a second..

myp[n_ /; PrimeQ[n/2]] := {n/2, n/2}  (* special case if n/2 is prime *)
myp[n_?EvenQ] := Module[{y = Prime@(PrimePi[n/2]+1)},
While[! AllTrue[{y, n - y}, PrimeQ], y = NextPrime[y]];
{y,n-y}]
hello = Table[{n, lowerbound[n], (mm = myp[n])[[1]], n/2, mm[[2]],
upperbound[n]}, {n, {10^2, 10^3, 10^4, 10^5, 10^6, 10^7, 10^8,
10^9, 10^10}}]; // AbsoluteTiming
TableForm[hello,
TableHeadings -> {None, {"n", "(n/2)-250", "MGPPPp", "n/2", "MGPPPq",
"(n/2)+250"}}, TableAlignments -> Center]


also in your table you were calling mgppp[n] twice for every n. Take note how I fixed that.. (factor of two timing savings right there )

aside, we can easily find numbers that break your bounds..

myp[#] - #/2 &@ 12573226
mgppp[#] - #/2 &@ 12573226


{-1140, 1140}

{-1140, 1140}

If we evaluate:

Clear[mgppp, lowerbound, upperbound]

mgppp[n_?EvenQ] /; n > 3 :=
Block[{m = PrimePi[n/2], p},
While[! PrimeQ[q = n - (p = Prime[m])], m--]; {p, q}]

lowerbound[n_?EvenQ] := Module[{d}, {m = n/2, d = (m - 250)}; d]


(** 10^2 to 10^3 range: **)

hello = Table[{n, lowerbound[n], mgppp[n][[1]], n/2, mgppp[n][[2]],
upperbound[n],
If[lowerbound[n] < mgppp[n][[1]], "ok", "Not bounded"]}, {n, 10^2,
10^3, 2}];

TableForm[hello,
TableHeadings -> {None, {"n", "(n/2)-250", "MGPPPp", "n/2", "MGPPPq",
"(n/2)+250", "Bounding Status"}}, TableAlignments -> Center]


(** 10^3 to 10^4 range: **)

hello2 = Table[{n, lowerbound[n], mgppp[n][[1]], n/2, mgppp[n][[2]],
upperbound[n],
If[lowerbound[n] < mgppp[n][[1]], "ok", "Not bounded"]}, {n, 10^3,
10^4, 2}];

TableForm[hello2,
TableHeadings -> {None, {"n", "(n/2)-250", "MGPPPp", "n/2", "MGPPPq",
"(n/2)+250", "Bounding Status"}}, TableAlignments -> Center]


we find that the lower and upper bounds work for n in the above ranges.

If we evaluate the following:

(** 10^4 to 10^5 range: **)

hello3 = Table[
If[lowerbound[n] < mgppp[n][[1]],
Nothing, {n, lowerbound[n], mgppp[n][[1]], n/2, mgppp[n][[2]],
upperbound[n], "Unbounded"}], {n, 10^4, 10^5, 2}];

TableForm[hello3,
TableHeadings -> {None, {"n", "(n/2)-250", "MGPPPp", "n/2", "MGPPPq",
"(n/2)+250", "Bounding Status"}}, TableAlignments -> Center]


we find the the bounds fail for all the n values listed on the hello3 table.

The evaluations take longer when we reach the 10^5 to 10^6 range:

(** 10^5 to 10^6 range: **)

hello4 = Table[
If[lowerbound[n] < mgppp[n][[1]],
Nothing, {n, lowerbound[n], mgppp[n][[1]], n/2, mgppp[n][[2]],
upperbound[n], "Unbounded"}], {n, 10^5, 10^6, 2}];

TableForm[hello4,
TableHeadings -> {None, {"n", "(n/2)-250", "MGPPPp", "n/2", "MGPPPq",
"(n/2)+250", "Bounding Status"}}, TableAlignments -> Center]


For larger ranges the evaluations take too long.

The n values below correspond to sequence A065978 at: http://oeis.org/A065978

The lowerbound values below correspond to sequence A155764 at: http://oeis.org/A155764

 mgppp[n_?EvenQ] /; n > 3 :=
Block[{m = PrimePi[n/2], p},
While[! PrimeQ[n - (p = Prime[m])], m--]; p]

lowerbound[n_?EvenQ] := Module[{d}, {m = n/2, d = (m - mgppp[n])}; d]

bag = {}; n = 4; max = -1; While[n <= 10^5,
If[lowerbound[n] > max, AppendTo[bag, {n, lowerbound[n]}];
max = lowerbound[n]]; n++]; bag


Notice that I was attempting to write a program that verifies whether pre-selected lower and upper bounds for the MGPPP work for ALL even n between 10^2 and 10^10. The purpose of the examples with tables listing the lower and upper bounds for n's that are powers of 10, was to illustrate and give a simpler starting point to get at the ultimate task of doing this for ALL even n's in the above range. This is not the best approach to tackle this problem; a better approach is to collect the n's corresponding to the strictly-increasing distances n/2-MGPPP[n] as done in A065978 at http://oeis.org/A065978 and A155764 at http://oeis.org/A155764 where this approach is applied. Thank you all for your valuable help!