Message Boards Message Boards

1
|
18109 Views
|
16 Replies
|
14 Total Likes
View groups...
Share
Share this post:

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

Posted 11 years ago
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
POSTED BY: H FR
16 Replies
Joe, there are some issues with the editor. Please see reply from Vitaliy here for a few hints on easier editor handling:
POSTED BY: Moderation Team
Posted 11 years ago
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
POSTED BY: Joe Gilray
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. 
POSTED BY: Sander Huisman
Posted 11 years ago
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
POSTED BY: Joe Gilray
@chip
Indeed that's a nice trick! Thanks! Didn't think about it ;)

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...
POSTED BY: Sander Huisman
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 emoticon
POSTED BY: Sander Huisman
Posted 11 years ago
@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.
POSTED BY: Greg Hurst
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... 
POSTED BY: Sander Huisman
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.
POSTED BY: Sander Huisman
@Sander, 

Good one! I think the nature of this problem is pro-loop. It avoids many unnecessary steps by using Range or Table. 
POSTED BY: Shenghui Yang
"A picture is worth a thousand words" -> "A Wolfram Alpha function is worth a thousand command lines" 
POSTED BY: Shenghui Yang
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]



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

POSTED BY: Vitaliy Kaurov
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 :-)

How about something like
Rest[NestList[(NestWhile[# + 1 &, #, ! PrimeQ[#1] || ! PrimeQ[#3] &, 3] - 1) &, 1, 100] - 1]
POSTED BY: Ilian Gachevski
Posted 11 years ago
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
POSTED BY: H FR
The third argument of Select could be useful, for example
Select[Prime[Range[1000]], PrimeQ[# + 2] &, 100]
POSTED BY: Ilian Gachevski
Posted 11 years ago
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.
POSTED BY: H FR
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract