3
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – user3363795 Feb 15 '16 at 15:53
  • $\begingroup$ Here's a link to the Mathematica Notebook: dl.dropboxusercontent.com/u/76769933/TwinPrimes.nb $\endgroup$ – user3363795 Feb 15 '16 at 16:05
  • $\begingroup$ maybe instead of storing the primes in full form, store n where Prime[n] is the prime you are interested in. $\endgroup$ – Manuel --Moe-- G Feb 15 '16 at 17:54
  • $\begingroup$ 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 $\endgroup$ – Manuel --Moe-- G Feb 15 '16 at 17:56
  • $\begingroup$ forgive me if I have a poor understanding of what your code does $\endgroup$ – Manuel --Moe-- G Feb 15 '16 at 17:57
6
$\begingroup$

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} *)
$\endgroup$
  • $\begingroup$ YOU...ARE...AWESOME! This is sick! If this becomes useful, I'll credit you! $\endgroup$ – Elem-Teach-w-Bach-n-Math-Ed Feb 16 '16 at 1:26
  • $\begingroup$ 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. $\endgroup$ – Elem-Teach-w-Bach-n-Math-Ed Feb 16 '16 at 1:46
  • $\begingroup$ 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 $\endgroup$ – Elem-Teach-w-Bach-n-Math-Ed Feb 16 '16 at 1:48
  • $\begingroup$ Also, how about an option to run a set number of threads each testing a separate prime? $\endgroup$ – Elem-Teach-w-Bach-n-Math-Ed Feb 16 '16 at 1:56
  • $\begingroup$ 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. $\endgroup$ – Elem-Teach-w-Bach-n-Math-Ed Feb 16 '16 at 6:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.