How to find the first 100 prime twins without a do loop or iteration

Questions
1

I’m fairly new to programming with Mathematica and I’m trying to use Do loops and iteration as little as possible. One typical problem I run into is this:

The following code quickly finds all prime twins ( the first member of them) that are smaller than 100

Select[Range[100], (PrimeQ[#] && PrimeQ[# + 2]) &]

I would like to know if it is possible to use a similar simple line of code to find the first hundred prime twins without using a do or for loop and iteration, as I would normally do with a procedural language. I studied the manuals and several books, but I haven’t found a way to do this.

Any help would be appreciated.

Harrie

2

The third argument of Select could be useful, for example

Select[Prime[Range[1000]], PrimeQ[# + 2] &, 100]
3

Thanks for the reply, Ilian. Hadn’t thought of that.

Still, the code still needs some “advance” knowledge, because you have to tell it to generatie 1,000 primes first. But that assumes 1,000 consecutive primes will contain 100 prime twins and that may not be known in advance. In addition, it requires extra computation time. Do you think it would be possible to have a solution without the Range[1000], or is that impossible with a functional language?

Thanks again,

Harrie.

4

Right, to figure out the correct range to use we would need a lower bound on the twin prime counting function, which isn’t known or else we could have proved a famous conjecture :slight_smile:

How about something like

Rest[NestList[(NestWhile[# + 1 &, #, ! PrimeQ[#1] || ! PrimeQ[#3] &, 3] - 1) &, 1, 100] - 1]
5

I understand that this is a very targeted discussion to avoid usage of Range. But still just for the sake of general information I will notice the following. There is a curious coincidence - we had a very popular post about computing twin primes - maybe you can find something useful there: ** Ulam spiral to appear on Sept. 2013 cover of Math Horizons **. You could also use a function that calls Wolfram|Alpha to get those numbers in increments of 20:

TPP[k_] := 
 WolframAlpha["twin prime pair", {{"Result", 1}, "ComputableData"}, 
  PodStates -> Table["Result__More", {k}]]

TPP[13]

2013-10-28_21-10-14.jpg
2013-10-28_21-10-14.jpg545×477 204 KB

Which was basically generated from the following Wolfram|Alpha query: twin prime pair

2013-10-28_21-04-10.jpg

6

“A picture is worth a thousand words” → “A Wolfram Alpha function is worth a thousand command lines”

7

Thanks again, Ilian.

This one is still a bit over my head, but well worth studying. It does achieve my goal of avoiding a for or do loop and it also doesn’t use the range command.

Well done!

Regards,

Harrie

8

Intelligent Do loop works really fast:

Reap[
 Do[
  If[PrimeQ[Prime[i] + 2], Sow[{Prime[i], Prime[i] + 2}]]
  , {i, 1, PrimePi[100]}
 ]
]

First it figures out how many primes there are up to a 100 (PrimePi[100]) then iterates just over the primes and checks if nth prime + 2 is also a prime. Only if so the results are sown and reaped afterwards.

9

@Sander,

Good one! I think the nature of this problem is pro-loop. It avoids many unnecessary steps by using Range or Table.

10

Here is the full solution that will find the first m twins:

m=100;
n=1;
f=0;
res=Reap[
 While[f<m,
  If[PrimeQ[Prime[n]+2],Sow[Prime[n]];f++];
  n++;
 ];
][[2,1]];
twins={res,res+2}\[Transpose];

for m=10000 it takes 0.27 seconds on my computer, the solution by Ilian Gachevski takes 3.26 on my computer (so my code is about 12X faster).
How many are you looking for? I’m pretty sure you can’t go much faster in Mathematica in just a few lines. Parallelisation could speed things up quite a bit more…

11

O! In addition you could look at functions like NextPrime and check if if it two more than the previous one. I believe, though, that the function is not very fast. But you could try it for your self :slight_smile:

12

@Sander

One can speed up the loop by realizing there can be no twin prime triples other than (3, 5, 7).

That is, in your loop when we find a twin prime pair, do n += 2 rather than n++.

This offers about an 8% speed up for m=10000.

13

@chip
Indeed that’s a nice trick! Thanks! Didn’t think about it :wink:

In addition, once ould pre-allocate the results, and not use sow/reap, because the length is known beforehand. Most likely it will be slightly faster…

14

Harrie,

Here is an idiom that is useful when you don’t have explicit limits:

{{}, 3, 0} //. {r_, n_, c_ /; c < 100} :> If[PrimeQ[n] && PrimeQ[n + 2], 
  {Append[r, {n, n + 2}], n + 4, c + 1}, {r, n + 2, c}] // First

I’ve also found it to be extremely fast.

Looking up //. and /; and :> in the docs will also help with your general understanding of how Mathematica works.
-Joe

15

Dear Joe,

There is an error in your submission, it should be PrimeQ [ n ] , and for large number of primes, replaceall will give an error becuase MaxIterations is by default 2^16. So I changed the code of the 3 solutions now to measure the time needed to generate the primes and to measure computational duration for generating 100000 primes:

AbsoluteTiming[ReplaceRepeated[{{},3,0},{r_,n_,c_/;c<num}:>If[PrimeQ[n]&&PrimeQ[n+2],{Append[r,{n,n+2}],n+4,c+1},{r,n+2,c}],MaxIterations->\[Infinity]]//First;]

m=num;
n=1;
f=0;
AbsoluteTiming[res=Reap[While[f<m,If[PrimeQ[Prime[n]+2],Sow[Prime[n]];f++];n++;];][[2,1]];]
twins={res,res+2}\[Transpose];


AbsoluteTiming[Rest[NestList[(NestWhile[#+1&,#,!PrimeQ[#1]||!PrimeQ[#3]&,3]-1)&,1,num]-1];]

times:
{74.419518, Null}{5.666652, Null}{50.630102, Null}

So although the code with replaceall is very neat, it is slow compared to the others. The main problem with Joe’s solution is that primes get more and more separated for large numbers, so the stepping with 2 or 4 makes it slow.

16

Hi Sander,

Thanks for the education!

There is something funky about “publishing” code on this forum. I cut and pasted working code, but when published the was gone. I edited my post and it appears when editing, but not when published!? I finally got it to “take” by replacing with @n.

I guess I need to learn how to use the tags to get a code box.

-joe

17

Joe, there are some issues with the editor. Please see ** reply from Vitaliy here ** for a few hints on easier editor handling: