Iterate until condition is met

I want to find the first 5 prime numbers of the form $n^6 + 1091$. I have used this code:

Timing[Select[Table[n^6 + 1091, {n, 10000}], PrimeQ, 5]]

Which gives the desired answer, and takes 0.0156 seconds. What if I didn't have a rough idea of what values of $n$ to check?

Timing[Select[Table[n^6 + 1091, {n, 100000}], PrimeQ, 5]]

Takes 10 times as long.

So, my question is: Is there a simple(ish) way to make Mathematica stop when the desired number of solutions are found?

• Please don't be so hasty to Accept answers. I am working on one now and it is a discouragement to see an Accept given before I even have a chance to finish. Mar 18 '13 at 11:32
• Sorry, I am fairly new to this. I have unaccepted. Mar 18 '13 at 11:34
• Says the man who has the fastest gun in the west :) (I agree though not to accept answers too quickly) Mar 18 '13 at 11:34
• @IstvánZachar Fast? I just spent 45 minutes writing what is really a very simple answer. Mar 18 '13 at 12:17
• Thanks to everybody for the help, I have learnt something from each answer, and I would have chosen more than one as the official answer, if I could have. Mar 18 '13 at 17:12

A more efficient use of Select

If the likely bound of the problem is easily stored in memory it is practical to generate a Range, which is fast, and then Select from that. Since the Range will be unpacked by Select you must consider the unpacked size. For example:

ByteCount @ DeveloperFromPackedArray @ Range@1*^7

240 000 032

This is a reasonable starting size in many cases.

Finding five values take a quarter of a second:

Select[Range@1*^7, PrimeQ[#^6 + 1091] &, 5] // Length // Timing
{0.249, 5}

Fining all values with a seed <= 10^7:

Select[Range@1*^7, PrimeQ[#^6 + 1091] &] // Length // Timing
{19.797, 3338}

You can see that this scales pretty well, without excessive overhead in the case that values are found quickly, yet for your search 3338 results can be found before the limit is reached.

(The seeds found with this method need to be converted to primes with #^6 + 1091 &.)

A method for larger problems

If this will not work, either because you have no idea what the likely upper bound is, or it is too high to hold a Range in memory, then it will be most efficient to operate in blocks, due to Mathematica's vector optimizations.

First, there is a more efficient way to build the list of candidates:

Range^6 + 1091

This takes full advantage of the vectorized operations available.

Pick a large enough block size that the average element processing time is relatively low, but not so large as to process more elements than are likely needed. I will pick a block size of 100,000 and I will try to find 5,000 solutions:

block = 100000;

result = {};

find = 5000;

hits = 0;

Do[
If[hits >= find, Return["Done!"]];
Select[(n block + Range[block])^6 + 1091, PrimeQ] //
(result = {result, #}; hits += Length@#) &,
{n, 0, 1*^9}
] ~Monitor~ n

The results are stored in a linked list rather than using Sow. This has the advantage of letting you open a sub-session and examine the results up to that point. For example, in a separate cell enter:

Flatten @ result // Short

And use F7 to show the results, then resume calculation.

The Monitor lets you see how many blocks have been processed.

Timings

I made the claim that these methods are efficient. Let me give some comparative timings to support my position. WReach's lazy lists code which Rojo used in his answer is a wonderful approach. It is not however, as written, fast. My methods are by comparison clunky but they are also more practical.

I will use a variation of Timo's timeAvg function:

SetAttributes[timeAvg, HoldFirst]

timeAvg[func_] := Do[If[# > 0.3, Return[#/5^i]] & @@ Timing@Do[func, {5^i}], {i, 0, 15}]

Compared to my first method:

Table[
{n,
Select[Range@1*^7, PrimeQ[#^6 + 1091] &, n] // timeAvg,
sIntegers[] ~sMap~ (#^6 + 1091 &) ~sFilter~ PrimeQ ~sTake~ n // sList // timeAvg
},
{n, {5, 15, 50, 150, 500}}
] // TableForm[#, TableHeadings -> {None, {"n", "Select", "Lazy"}}] & My Select method shows the overhead of generating and unpacking the Range but it soon catches up and exceeds the lazy lists method in performance. Remember also the human overhead of writing this if it is not going to be used many times. The Select method is very simple and direct.

Now, for my second method there is a tuning parameter: the block size. It could be argued that changing this parameter mid-test is not fair play so I will use a fixed block size of 1000.

finder[find_, block_: 1000] := Module[{result = {}, hits = 0},
Do[
If[hits >= find, Return[Flatten@result ~Take~ find]];
Select[(n block + Range[block])^6 + 1091, PrimeQ] //
(result = {result, #}; hits += Length@#) &,
{n, 0, 1*^9}
]
]

Table[
{n,
finder[n] // timeAvg,
sIntegers[] ~sMap~ (#^6 + 1091 &) ~sFilter~ PrimeQ ~sTake~ n // sList // timeAvg
},
{n, {5, 15, 50, 150, 500, 1500, 5000}}
] // TableForm[#, TableHeadings -> {None, {"n", "Finder", "Lazy"}}] & Here the superiority of the block-based method is fully apparent.

• +1 for the useful tutorial! May I ask why you are using nested lists in result= {result, #} instead of AppendTo? Is it to keep track of how far you've gone by looking at the depth of result?
– gpap
Mar 18 '13 at 12:37
• @gpap Linked lists are faster than AppendTo, yet provide the option of viewing intermediate results (unlike Sow). See the bottom of this page for a reference. Faster is Sow/Reap. The ultimate solution is InternalBag and friends, which is slightly faster than Sow/Reap and yet makes it possible to view results, but I didn't think that level of optimization was needed here. Mar 18 '13 at 12:46
• thanks for that - favorited a couple for further reading.
– gpap
Mar 18 '13 at 12:59
• Splitting it into blocks like that is clever, and, I think, very readable, so I am selecting this as the answer. Mar 18 '13 at 17:11
• @nick Thank you, I hope this method serves you well. If you are handling a very large problem you'll want to save results as you go, like this. Mar 19 '13 at 0:54

Using @WReach's lazy lists, you can neatly solve this

sIntegers[] ~sMap~ (#^6 + 1091 &) ~sFilter~ PrimeQ ~sTake~ 5 // sList
• Shiny. I need to make better use of those constructs. +1 Mar 18 '13 at 12:50
• By the way, what does the "s" stand for? Slothful? :^) Mar 18 '13 at 13:00
• +1. I did somewhat similar thing, but conceptually yor way is better - if one goes with lazy lists, one should do it consistently. Mar 18 '13 at 13:13
• @Mr.Wizard haha, he called them streams
– Rojo
Mar 18 '13 at 13:20
• Ah, I missed that. Also, just a reminder, please add extra space around ~infix~ operators, especially when they are not built-in functions. This greatly improves the legibility of the code. Mar 18 '13 at 13:21

Use While:

primes = {}; i = 1;
AbsoluteTiming[While[Length@primes < 10,
i++; n = i^6 + 1091; If[PrimeQ@n, AppendTo[primes, n]]]; primes]
{0.046800, {3551349655007944406147, 9724154565432384001091,
19007561792329201176707, 24141016911038304486467,
26494507823224896001091, 301789003173921081001091,
4871535877925849664001091, 7484471576562177870314627,
9255722232902778801001091, 33270780762388607695743107}}

For comparison:

AbsoluteTiming[Select[Table[n^6 + 1091, {n, 100000}], PrimeQ, 10]]
{0.234001, {3551349655007944406147, 9724154565432384001091,
19007561792329201176707, 24141016911038304486467,
26494507823224896001091, 301789003173921081001091,
4871535877925849664001091, 7484471576562177870314627,
9255722232902778801001091, 33270780762388607695743107}}

You can also break out of inbuilt functions, though in that case I suggest using Do instead of Table (no use of building a table if it will be aborted before completion), and sowing found solutions on the fly with Sow & Reap:

c = 0;
AbsoluteTiming[
First@Last@
Reap@Do[If[PrimeQ[n^6 + 1091], Sow[n^6 + 1091]; c++;
If[c == 10, Break[]]], {n, 100000}]]
{0.062400, {3551349655007944406147, 9724154565432384001091,
19007561792329201176707, 24141016911038304486467,
26494507823224896001091, 301789003173921081001091,
4871535877925849664001091, 7484471576562177870314627,
9255722232902778801001091, 33270780762388607695743107}}
• Thanks István! I was looking for a more "Mathematica" way though. Using procedural approaches tends to slow things down a lot ie. This takes 3 times as long as the initial solution. Mathematica clearly has faster inbuilt code, it's just a case of getting it to stop. Mar 18 '13 at 11:19
• Well, your problem specifically ask for a procedural break when condition is met. If you check timings, you will see that it depends on your range of scope (10000), and thus While is indeed faster. Mar 18 '13 at 11:21
• Yep. You're right. Is there a way to get the inbuit algorithms to stop, though? Mar 18 '13 at 11:22
• @nickjamesuk since you have a PrimeQ in the loop, it's very unlikely that an extra, slow addition will be more expensive than that. Using Do or While or whatever should not make a difference, speed-wise.
– acl
Mar 18 '13 at 11:23
• I like the Sow and Reap thing, I have always wanted to work out how they are used. Thanks so much for all your work. Mar 18 '13 at 11:30

You can use NestWhile but I doubt it will be much (any?) faster than the one @IstvánZachar gave you.

prime = NestWhile[# + 1 &, # + 1, ! PrimeQ@((# + 1)^6 + 1091) &] &
((# + 1)^6 + 1091) & /@ NestList[prime, 1, 6]~Drop~1

The pure function prime is easier to read if you change the middle argument of NestWhile i.e.

NestWhile[# + 1 &, 2, ! PrimeQ@((# + 1)^6 + 1091) &]

looks for the number (after 2) for which $\text{number}^6 + 1091$ is prime. So, turning this into a pure function, you can then apply it repeatedly for every number it spits out using NestList.

For some reason NestList produces its zeroth evaluation (the number 1) and I want to get rid of that from the list using Drop[ListofNumbersINeed,1] which in Infix notation is ListofNumbersINeed~Drop~1. (I learnt that from @MrWizard and it sometimes makes code easier to understand - not the case here:) ).

I then apply the pure function ((# + 1)^6 + 1091) & to the list to get the actual primes. There is bound to be a better way to do the whole thing in one go or with much more obscure notation or in a hundredth of the speed - give it a couple of hours/days and you'll see what I mean.

• I think my Mathematica knowledge is not yet at the stage where I can parse that. What does the ~Drop~ do at the end? Mar 18 '13 at 11:33
• I edited to make it a little clearer - NestWhile, NestList, Function and Infix are the relevant pages.
– gpap
Mar 18 '13 at 11:54
• Ah, so the ~Drop~ is Infix notation. Thanks very much for the clarification. Mar 18 '13 at 14:06
• @nickjamesuk I'm largely the one to blame for the proliferation of that syntax around here. It apparently trips up a lot of people upon first seeing it, but after extensive use I find that it is often easier to read than the nested-bracket nightmares that can be the alternative. If you use a string of infix, as shown in Rojo's answer, it is very important IMHO to add space around each operator to distinguish these from arguments. Once this is done however code reads left to right in a logical fashion, e.g. a ~p~ b ~q~ c ~r~ d, which is equivalent to r[q[p[a, b], c], d]. Mar 18 '13 at 14:35
• That doesn make sense. I'll be looking out for ways to utilise it. Mar 18 '13 at 17:06

Inspired by gpap's answer here is a use of NestWhileList that I cannot recommend it practice but which I had fun writing. (It's much slower than the vectorized methods in my other answer.)

n = 5;

NestWhileList[
If[PrimeQ[(++i)^6 + 1091], i, ## &[]] &,
i = 1,
Length@{##} < n &,
All
]

The "vanishing function" ##&[] is used to silently remove any steps that don't yield a solution, and the All parameter passes all found solutions to the length check function, which stops the loop when n solutions are found.