Iterate until condition is met

Question

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?

Answer 1

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 @ Developer`FromPackedArray @ 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[100]^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.

Answer 2

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

sIntegers[] ~sMap~ (#^6 + 1091 &) ~sFilter~ PrimeQ ~sTake~ 5 // sList

Answer 3

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

Answer 4

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.

Answer 5

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.