Generate a list in Mathematica with a conditional

Question

Suppose we want to generate a list of primes p for which p + 2 is also prime.

A quick solution is to generate a complete list of the first n primes and use the Select function to return the elements which meet the condition.

Select[Table[Prime[k], {k, n}], PrimeQ[# + 2] &]

However, this is inefficient as it loads a large list into the memory before returning the filtered list. A For loop with Sow/Reap (or l = {}; AppendTo[l, k]) solves the memory issue, but it is far from elegant and is cumbersome to implement a number of times in a Mathematica script.

Reap[
  For[k = 1, k <= n, k++,
   p = Prime[k];
   If[PrimeQ[p + 2], Sow[p]]
  ]
 ][[-1, 1]]

An ideal solution would be a built-in function which allows an option similar to this.

Table[Prime[k], {k, n}, AddIf -> PrimeQ[# + 2] &]

Answer 1

You can perhaps try something like this:

Clear[f, primesList]
f = With[{p = Prime[#]},Piecewise[{{p, PrimeQ[p + 2]}}, {}] ] &;
primesList[k_] := Union@Flatten@(f /@ Range[k]);

---

If you want both the prime p and the prime p+2, then the solution is

Clear[f, primesList]
f = With[{p = Prime[#]},Piecewise[{{p, PrimeQ[p + 2]}}, {}] ] &;
primesList[k_] := 
  Module[{primes = f /@ Range[k]}, 
   Union@Flatten@{primes, primes + 2}];

Answer 2

I think the Reap/Sow approach is likely to be most efficient in terms of memory usage. Some alternatives might be:

DeleteCases[(With[{p=Prime[#]},If[PrimeQ[p+2],p,{}] ] ) & /@ Range[K]),_List]

Or (this one might need some sort of DeleteCases to eliminate Null results):

FoldList[[(With[{p=Prime[#2]},If[PrimeQ[p+2],p] ] )& ,1.,Range[2,K] ]

Both hold a big list of integers 1 to K in memory, but the Primes are scoped inside the With[] construct.

Answer 3

Well, someone has to allocate memory somewhere for the full table size, since it is not known before hand what the final size will be.

In the good old days before functional programming :), this sort of thing was solved by allocating the maximum array size, and then using a separate index to insert to it so no holes are made. Like this

x=Table[0,{100}];  (*allocate maximum possible*)
j=0;
Table[ If[PrimeQ[k+2], x[[++j]]=k],{k,100}];

x[[1;;j]]  (*the result is here *)

{1,3,5,9,11,15,17,21,27,29,35,39,41,45,51,57,59,65,69,71,77,81,87,95,99}

Answer 4

Yes, this is another answer. Another alternative that includes the flavour of the Reap/Sow approach and the FoldList approach would be to use Scan.

result = {1};
Scan[With[{p=Prime[#]},If[PrimeQ[p+2],result={result,p}]]&,Range[2,K] ];
Flatten[result]

Again, this involves a long list of integers, but the intermediate Prime results are not stored because they are in the local scope of With. Because p is a constant in the scope of the With function, you can use With rather than Module, and gain a bit of speed.

Answer 5

I will interpret this more as a question about automation and software engineering rather than about the specific problem at hand, and given a large number of solutions posted already. Reap and Sow are good means (possibly, the best in the symbolic setting) to collect intermediate results. Let us just make it general, to avoid code duplication.

What we need is to write a higher-order function. I will not do anything radically new, but will simply package your solution to make it more generally applicable:

Clear[tableGen];
tableGen[f_, iter : {i_Symbol, __}, addif : Except[_List] : (True &)] :=
 Module[{sowTag},   
  If[# === {}, #, First@#] &@
       Last@Reap[Do[If[addif[#], Sow[#,sowTag]] &[f[i]], iter],sowTag]];

The advantages of using Do over For are that the loop variable is localized dynamically (so, no global modifications for it outside the scope of Do), and also the iterator syntax of Do is closer to that of Table (Do is also slightly faster).

Now, here is the usage

In[56]:= tableGen[Prime, {i, 10}, PrimeQ[# + 2] &]

Out[56]= {3, 5, 11, 17, 29}

In[57]:= tableGen[Prime, {i, 3, 10}, PrimeQ[# + 1] &]

Out[57]= {}

In[58]:= tableGen[Prime, {i, 10}]

Out[58]= {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}

EDIT

This version is closer to the syntax you mentioned (it takes an expression rather than a function):

ClearAll[tableGenAlt];
SetAttributes[tableGenAlt, HoldAll];
tableGenAlt[expr_, iter_List, addif : Except[_List] : (True &)] :=
 Module[{sowTag}, 
  If[# === {}, #, First@#] &@
    Last@Reap[Do[If[addif[#], Sow[#,sowTag]] &[expr], iter],sowTag]];

It has an added advantage that you may even have iterator symbols defined globally, since they are passed unevaluated and dynamically localized. Examples of use:

In[65]:= tableGenAlt[Prime[i], {i, 10}, PrimeQ[# + 2] &]

Out[65]= {3, 5, 11, 17, 29}

In[68]:= tableGenAlt[Prime[i], {i, 10}]

Out[68]= {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}

Note that since the syntax is different now, we had to use the Hold-attribute to prevent the passed expression expr from premature evaluation.

EDIT 2

Per @Simon's request, here is the generalization for many dimensions:

ClearAll[tableGenAltMD];
SetAttributes[tableGenAltMD, HoldAll];
tableGenAltMD[expr_, iter__List, addif : Except[_List] : (True &)] :=
Module[{indices, indexedRes, sowTag},
  SetDelayed @@  Prepend[Thread[Map[Take[#, 1] &, List @@ Hold @@@ Hold[iter]], 
      Hold], indices];
  indexedRes = 
    If[# === {}, #, First@#] &@
      Last@Reap[Do[If[addif[#], Sow[{#, indices},sowTag]] &[expr], iter],sowTag];
  Map[
    First, 
    SplitBy[indexedRes , 
      Table[With[{i = i}, Function[Slot[1][[2, i]]]], {i,Length[Hold[iter]] - 1}]], 
    {-3}]];

It is considerably less trivial, since I had to Sow the indices together with the added values, and then split the resulting flat list according to the indices. Here is an example of use:

{i, j, k} = {1, 2, 3};
tableGenAltMD[i + j + k, {i, 1, 5}, {j, 1, 3}, {k, 1, 2}, # < 7 &]

{{{3, 4}, {4, 5}, {5, 6}}, {{4, 5}, {5, 6}, {6}}, {{5, 6}, {6}}, {{6}}}

I assigned the values to i,j,k iterator variables to illustrate that this function does localize the iterator variables and is insensitive to possible global values for them. To check the result, we may use Table and then delete the elements not satisfying the condition:

In[126]:= 
DeleteCases[Table[i + j + k, {i, 1, 5}, {j, 1, 3}, {k, 1, 2}], 
    x_Integer /; x >= 7, Infinity] //. {} :> Sequence[]

Out[126]= {{{3, 4}, {4, 5}, {5, 6}}, {{4, 5}, {5, 6}, {6}}, {{5, 6}, {6}}, {{6}}}

Note that I did not do extensive checks so the current version may contain bugs and needs some more testing.

EDIT 3 - BUG FIX

Note the important bug-fix: in all functions, I now use Sow with a custom unique tag, and Reap as well. Without this change, the functions would not work properly when expression they evaluate also uses Sow. This is a general situation with Reap-Sow, and resembles that for exceptions (Throw-Catch).

EDIT 4 - SyntaxInformation

Since this is such a potentially useful function, it is nice to make it behave more like a built-in function. First we add syntax highlighting and basic argument checking through

SyntaxInformation[tableGenAltMD] = {"ArgumentsPattern" -> {_, {_, _, _., _.}.., _.},
                                    "LocalVariables" -> {"Table", {2, -2}}};

Then, adding a usage message allows the menu item "Make Template" (Shift+Ctrl+k) to work:

tableGenAltMD::usage = "tableGenAltMD[expr,{i,imax},addif] will generate \
a list of values expr when i runs from 1 to imax, \
only including elements if addif[expr] returns true.
The default of addiff is True&."

A more complete and formatted usage message can be found in this gist.

Answer 6

Here's another couple of alternatives using NextPrime:

pairs1[pmax_] := Select[Range[pmax], PrimeQ[#] && NextPrime[#] == 2 + # &]

pairs2[pnum_] := Module[{p}, NestList[(p = NextPrime[#];
                      While[p + 2 != (p = NextPrime[p])]; 
                      p - 2) &, 3, pnum]] 

and a modification of your Reap/Sow solution that lets you specify the maximum prime:

pairs3[pmax_] := Module[{k,p},
                   Reap[For[k = 1, (p = Prime[k]) <= pmax, k++,
                        If[PrimeQ[p + 2], Sow[p]]]][[-1, 1]]]

The above are in order of increasing speed.

In[4]:= pairs2[10000]//Last//Timing
Out[4]= {3.48,1261079}
In[5]:= pairs1[1261079]//Last//Timing
Out[5]= {6.84,1261079}
In[6]:= pairs3[1261079]//Last//Timing
Out[7]= {0.58,1261079}