Question about how to speed up Mathematica code

Question

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.

Answer 1

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]

RuntimeTools`Profile[
  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]

RuntimeTools`Profile[
  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.

Answer 2

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}

Answer 3

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.

Answer 4

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

Answer 5

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!