# Memory limit hit: optimize code for finding twin primes

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.

## migrated from math.stackexchange.comFeb 15 '16 at 16:47

This question came from our site for people studying math at any level and professionals in related fields.

• Where the memory maxes out is when I enter 'testPrime[113, 1] for cell 2, or i = 29, i <= 29 in cell 3 for the 29th prime, 113. I can even have cell 3 run 'i = 4, i <= 28 with no issues, but for some reason it can't do the 29th prime or above without maxing out memory. I'm thinking there may be a way to have it run cell 2 in a fragmented manner where each prime product p is run independently, spits out a true or false for whether it creates a twin prime, then checks the next. Perhaps then a memory limit can be set on each individual fragment? – user3363795 Feb 15 '16 at 15:53
• Here's a link to the Mathematica Notebook: dl.dropboxusercontent.com/u/76769933/TwinPrimes.nb – user3363795 Feb 15 '16 at 16:05
• maybe instead of storing the primes in full form, store n where Prime[n] is the prime you are interested in. – Manuel --Moe-- G Feb 15 '16 at 17:54
• But I think the multiplying out is the memory killer. Just have Prime[n] and Prime[n+1], where Prime[n] + 2 == Prime[n+1], and do division tests using Mod for the list of factors you require – Manuel --Moe-- G Feb 15 '16 at 17:56
• forgive me if I have a poor understanding of what your code does – Manuel --Moe-- G Feb 15 '16 at 17:57

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} *)

• YOU...ARE...AWESOME! This is sick! If this becomes useful, I'll credit you! – Elem-Teach-w-Bach-n-Math-Ed Feb 16 '16 at 1:26
• Couple thoughts... Elsewhere it's been mentioned, "The product of all twin primes is an even perfect square minus 1. Any factor (prime or not) less than the square root of X^2 - 1 eliminates the composite as a candidate. Furthermore, these factors are usually trivial - that is, they're found long before X. Demonstration: naturalnumbers.org/TwinPrimeCalc.xlsm" Don't know what the code uses to check, but this may be handy. – Elem-Teach-w-Bach-n-Math-Ed Feb 16 '16 at 1:46
• Also, I added this to a cell at the bottom, and look out! 8000 different twin prime pairs all mapped to a single prime! EPIC! dl.dropboxusercontent.com/u/76769933/8000%20twin%20primes.txt – Elem-Teach-w-Bach-n-Math-Ed Feb 16 '16 at 1:48
• Also, how about an option to run a set number of threads each testing a separate prime? – Elem-Teach-w-Bach-n-Math-Ed Feb 16 '16 at 1:56
• One more request... How about another cell at the bottom that would list all matching twin primes for a prime input rather than stopping at the first. – Elem-Teach-w-Bach-n-Math-Ed Feb 16 '16 at 6:13