Memory limit hit: optimize code for finding twin primes

Question

I need some help optimizing a Mathematica code so that it'll not max out RAM. Here's the code:

Cell 1

twinPrimesQ[tp_] := 
  tp[[1]] + 2 == tp[[2]] && PrimeQ[tp[[1]]] && PrimeQ[tp[[2]]];
primesList[p_] := Module[{out = {Prime[3]}, i},
   For[i = 4, Prime[i] <= p, i = i + 1,
    out = Append[out, Prime[i]];
    ];
   out
   ];
testPrime[p_, pl_] := 
  Module[{out, found = False, twinPrimes, primeFactors, 
    primeFactorsPowers, i},
   twinPrimes = {
     {3 5 prod p - 4, 3 5 prod p - 2},
     {3 5 prod p + 2, 3 5 prod p + 4}
     };
   primeFactors = primesList[p];
   primeFactorsPowers = Tuples[Range[0, pl], primeFactors // Length];
   For[i = 1, i <= Length[primeFactorsPowers], i = i + 1,
    out = 
     twinPrimes /. 
      prod -> Product[
        primeFactors[[k]]^primeFactorsPowers[[i]][[k]], {k, 1, 
         primeFactors // Length}];
    found = twinPrimesQ[out[[1]]] || twinPrimesQ[out[[2]]];
    If[found, Break[]];
    ];
   If[found, out~Select~(twinPrimesQ[#] &) // First, False]
   ];

Cell 2

testPrime[109, 1]

Cell 3

out = List[]; For[j = 30, j <= 30, j = j + 1, 
out = out~
   Append~{j, Prime[j], testPrime[Prime[j], 1]}]; out // TableForm

The code returns the first twin prime found that fits {3 5 prod p - 4, 3 5 prod p - 2}, {3 5 prod p + 2, 3 5 prod p + 4} where $p$ is a product of primes between 3 and the first input in cells 2 and 3.

Thanks for your help in advance!

Answer 1

Here is a rewritten version. I hesitate to call it optimised as I suspect there are significant further improvements that could be made.

The major change is to remove the Tuples bit (which was using all the memory) and replace it with a call to IntegerDigits inside the loop. I also changed the way prod is computed, using the Listable attribute of Power and also replacing Product with Times.

There are a couple of other tweaks, but nothing important for performance. I replaced For with Do. I prefer Throw & Catch over Break, and I added a maximum number of iterations for the loop, to avoid silliness.

I also altered primesList to make use of the Listable attribute of Prime.

twinPrimesQ[{a_, b_}] := a + 2 == b && PrimeQ[a] && PrimeQ[b];

primesList[p_] := Prime@Range[3, PrimePi@p]

testPrime[p_, pl_] := 
  Module[{twinPrimes, primeFactors, n, prod, maxiters = 10^7}, 
   twinPrimes = {{3 5 prod p - 4, 3 5 prod p - 2}, {3 5 prod p + 2, 3 5 prod p + 4}};
   primeFactors = primesList[p];
   n = Length[primeFactors];
   Catch[Do[
     prod = Times @@ (primeFactors^IntegerDigits[i, pl + 1, n]);
     If[twinPrimesQ[#], Throw[#]] & /@ twinPrimes,
     {i, Min[(pl + 1)^n - 1, maxiters]}]]];


(* example *)

testPrime[Prime[100], 5]
(* {26993431542303411984765525701, 26993431542303411984765525703} *)