# Functional style using lazy lists?

Let's say I want to answer the question "what are the first 400 palindromic prime numbers?"

The first approach that comes to my mind from the set of languages that I know is to use some sort of lazy list materialization, a la IEnumerable (and yield) in C#, generators in Python, or sequence blocks in F#.

For example, in C#:

PrimesEnumerator().Where(n => n.ToString() == n.ToString().Reverse()).Take(400);


This would cause the PrimesGenerator to be pumped for primes long enough for the Where() clause to find enough numbers that meet the requirement for Take() to meet its quota.

The best I've come up with in Mathematica is something like:

i = 1; results = List[];
While[Length[results] != 400,
If[IntegerDigits[Prime[i]] == Reverse[IntegerDigits[Prime[i]]],
results = Append[results,Prime[i]]];
i = i + 1]


It surprises me that I end up writing in such an imperative style in Mathematica. Am I missing something that would let me write this entirely functionally? Maybe even with lazy lists?

Update: I took inspiration from WReach's work of art answer, and made a package that took his ideas and expanded them into a broad, general solution for lazy data in Mathematica. I describe its usage in an answer below.

-
The code could be streamlined a bit: For[i = 1; k = 1; results = {}, k <= 300, i++, If[IntegerDigits[Prime[i]] == Reverse[IntegerDigits[Prime[i]]], results = {results, Prime[i]}; k++]]; Flatten[results]. –  Guess who it is. Jan 28 '12 at 3:46
Seems like most solutions to a "find the first N" type of problem suffers from "guess how big of a Range[] you need to iterate through". It seems, at least in the general case, that there's no great way to avoid that. Maybe I need to write a "Lazy Generators" package for Mathematica. –  sblom Jan 28 '12 at 4:04
Roman Maeder gave an implementation of lazy lists in one of his books (Mathematica Programmer). I used a technique described by @Sal, for a similar purposes: iterators - stackoverflow.com/questions/8596134/…, mathematica.stackexchange.com/questions/559/… (Module -> Advanced uses), more examples and links in my third post here: groups.google.com/group/comp.soft-sys.math.mathematica/…. –  Leonid Shifrin Jan 28 '12 at 11:11
@J.M., why results = {results, Prime[i]}; Flatten[results] instead of results = Append[results,Prime[i]]? I'm guessing it's to defer reallocation costs for results? –  sblom Jan 28 '12 at 17:38
@LeonidShifrin, that comment seems worthy of a credited answer. You'd certainly get an upvote from at least me. –  sblom Jan 28 '12 at 17:39

A "lazy list", "functional style" solution to this problem might look something like this:

sIntegers[] ~sMap~ Prime ~sFilter~ palindromicQ ~sTake~ 400 // sList


No such notation is built into Mathematica. However, creating such notations is Mathematica's strong suit. Let's do it.

First, we need to define the notion of a "stream". Streams are inherently lazy, so let's use HoldAll:

SetAttributes[stream, {HoldAll}]


A stream can be empty:

sEmptyQ[stream[]] := True


... or it can be non-empty, having two elements:

sEmptyQ[stream[_, _]] = False;


The first element of the stream is called the "head":

sHead[stream[h_, _]] := h


The remaining elements of the stream are called the "tail":

sTail[stream[_, t_]] := t


Armed with these definitions, we can now express an infinite stream of integers thus:

sIntegers[n_:1] :=
With[{nn = n+1}, stream[n, sIntegers[nn]]]

sIntegers[] // sEmptyQ                 (* False *)
sIntegers[] // sHead                   (* 1 *)
sIntegers[] // sTail // sHead          (* 2 *)
sIntegers[] // sTail // sTail // sHead (* 3 *)


Infinite streams are difficult to display in a notebook. Let's introduce sTake which truncates a stream to a fixed length:

sTake[s_stream, 0] := stream[]
sTake[s_stream, n_] /; n > 0 :=
With[{nn = n-1}, stream[sHead[s], sTake[sTail[s], nn]]]


Let's also introduce sList, which converts a (finite) stream into a list:

sList[s_stream] :=
Module[{tag}
, Reap[
, tag
][[2]] /. {l_} :> l
]


Now we can inspect an integer stream directly:

sIntegers[] ~sTake~ 10 // sList
(* {1, 2, 3, 4, 5, 6, 7, 8, 9, 10} *)


sMap applies a function to every element of a stream:

sMap[stream[], _] := stream[]
sMap[s_stream, fn_] := stream[fn[sHead[s]], sMap[sTail[s], fn]]

sIntegers[] ~sMap~ Prime ~sTake~ 10 // sList
(* {2, 3, 5, 7, 11, 13, 17, 19, 23, 29} *)


sFilter selects elements from a stream that satisfy a given filter predicate:

sFilter[s_, pred_] :=
NestWhile[sTail, s, (!sEmptyQ[#] && !pred[sHead[#]])&] /.
stream[h_, t_] :> stream[h, sFilter[t, pred]]

sIntegers[] ~sFilter~ OddQ ~sTake~ 15 // sList
(* {1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29} *)


We now have almost all of the pieces in place to address the original problem. All that is missing is a predicate that detects palindromic numbers:

palindromicQ[n_] := IntegerDigits[n] /. d_ :> d === Reverse[d]

palindromicQ[123] (* False *)
palindromicQ[121] (* True *)


Now, we can solve the problem:

sIntegers[] ~sMap~ Prime ~sFilter~ palindromicQ ~sTake~ 400 // sList

(* {2,3,5,7,11,101, ... ,3528253,3541453,3553553,3558553,3563653,3569653} *)


The stream facility we have defined here is very basic. It lacks error checking, and further consideration should be given to optimization. However, it demonstrates the power of Mathematica's symbolic programming paradigm.

The following listing gives the complete set of definitions:

ClearAll[stream]
SetAttributes[stream, {HoldAll, Protected}]

sEmptyError[] := (Message[stream::empty]; Abort[])
stream::empty = "Attempt to access beyond the end of a stream.";

ClearAll[sEmptyQ, sHead, sTail, sTake, sList, sMap, sFilter, sIntegers]

sEmptyQ[stream[]] := True
sEmptyQ[stream[_, _]] = False;

sTail[stream[]] := sEmptyError[]
sTail[stream[_, t_]] := t

sTake[s_stream, 0] := stream[]
sTake[s_stream, n_] /; n > 0 :=
With[{nn = n-1}, stream[sHead[s], sTake[sTail[s], nn]]]

sList[s_stream] :=
Module[{tag}
, Reap[
, tag
][[2]] /. {l_} :> l
]

sMap[stream[], _] := stream[]
sMap[s_stream, fn_] := stream[fn[sHead[s]], sMap[sTail[s], fn]]

sFilter[s_, pred_] :=
NestWhile[sTail, s, (!sEmptyQ[#] && !pred[sHead[#]])&] /.
stream[h_, t_] :> stream[h, sFilter[t, pred]]

sIntegers[n_:1] :=
With[{nn = n+1}, stream[n, sIntegers[nn]]]

palindromicQ[n_] := IntegerDigits[n] /. d_ :> d === Reverse[d]

-
I hoped you would answer. Your approach is, in fact, very similar to the approach taken by Roman Maeder in his implementation of lazy lists and streams in his "Mathematica programmer" book, except that he overloaded built-ins like Map etc with UpValues, while you created new heads. Very nice discussion and code - +1. –  Leonid Shifrin Jan 28 '12 at 23:37
@Leonid Thanks for the +1. I was surprised that I got through writing this one without you getting there first :) When I started out, I had Common Lisp series in mind. But it ended up more like a toy Haskell :) –  WReach Jan 29 '12 at 0:30
+1 ! Did you have this on the shelf or did you write it for this question? –  Mr.Wizard Jan 29 '12 at 5:14
Color me impressed. I tried to implement something like that years ago and failed. I might have been able to do it now but not that quickly or cleanly. –  Mr.Wizard Jan 29 '12 at 5:24
The syntax with sblom's package of WReach's solutions is Prime~Map~Lazy[Integers]~Select~palindromicQ~Take~30//List –  faysou Apr 9 '12 at 20:20

One way to get the lazy aspect is to use a closure, or the closest way for Mathematica to fake a closure.

This is the closures constructor:

makePalindromePrimeC[start_: 1] := Module[{p = Prime[start], r},
((r = NestWhile[NextPrime, p,
With[{d = IntegerDigits[#]}, d != Reverse[d]] &]);
p = NextPrime[r]; r) &]


This creates one:

palPrimeClosure = makePalindromePrimeC[]


Now you use it to generate some:

In[259]:= Table[palPrimeClosure[], {100}]

Out[259]= {2, 3, 5, 7, 11, 101, 131, 151, 181, 191, 313, 353, 373, \
383, 727, 757, 787, 797, 919, 929, 10301, 10501, 10601, 11311, 11411, \
12421, 12721, 12821, 13331, 13831, 13931, 14341, 14741, 15451, 15551, \
16061, 16361, 16561, 16661, 17471, 17971, 18181, 18481, 19391, 19891, \
19991, 30103, 30203, 30403, 30703, 30803, 31013, 31513, 32323, 32423, \
33533, 34543, 34843, 35053, 35153, 35353, 35753, 36263, 36563, 37273, \
37573, 38083, 38183, 38783, 39293, 70207, 70507, 70607, 71317, 71917, \
72227, 72727, 73037, 73237, 73637, 74047, 74747, 75557, 76367, 76667, \
77377, 77477, 77977, 78487, 78787, 78887, 79397, 79697, 79997, 90709, \
91019, 93139, 93239, 93739, 94049}


Generate some more:

In[260]:= Table[palPrimeClosure[], {50}]

Out[260]= {94349, 94649, 94849, 94949, 95959, 96269, 96469, 96769, \
97379, 97579, 97879, 98389, 98689, 1003001, 1008001, 1022201, \
1028201, 1035301, 1043401, 1055501, 1062601, 1065601, 1074701, \
1082801, 1085801, 1092901, 1093901, 1114111, 1117111, 1120211, \
1123211, 1126211, 1129211, 1134311, 1145411, 1150511, 1153511, \
1160611, 1163611, 1175711, 1177711, 1178711, 1180811, 1183811, \
1186811, 1190911, 1193911, 1196911, 1201021, 1208021}


Now create an entirely independent instance that starts searching at the 500th prime:

In[261]:= palPrimeClosure500 = makePalindromePrimeC[500]

Out[261]= (r$10054 = NestWhile[NextPrime, p$10054,
With[{d = IntegerDigits[#1]}, d != Reverse[d]] &];
p$10054 = NextPrime[r$10054]; r$10054) & In[262]:= Table[palPrimeClosure500[], {30}] Out[262]= {10301, 10501, 10601, 11311, 11411, 12421, 12721, 12821, \ 13331, 13831, 13931, 14341, 14741, 15451, 15551, 16061, 16361, 16561, \ 16661, 17471, 17971, 18181, 18481, 19391, 19891, 19991, 30103, 30203, \ 30403, 30703}  - It is indeed much faster then a blunt approach scanning all the numbers. – Vitaliy Kaurov Jan 28 '12 at 6:47 That's the closest thing I've seen so far. And I can even imagine getting rid of the NextPrime part by using another closure-as-a-generator for the prime numbers within the palindrome filter. – sblom Jan 28 '12 at 7:14 Wow! Nice. I should add that trick to my toolbox. – kkm Jan 28 '12 at 10:38 +1. I use similar techniques. I linked to some additional examples in the comment above. – Leonid Shifrin Jan 28 '12 at 11:18 I took inspiration from WReach's work of art answer, and made a package that took his ideas and expanded them into a broad, general solution for lazy data in Mathematica. You can find my implemenation on github. To use the package to answer the original question from this post, you'd do something like: palindromicQ[n_] := IntegerDigits[n] /. d_ :> d === Reverse[d] Needs["Lazy"] Lazy[Primes] ~Select~ palindromicQ ~Take~ 400 // List  It can also do things like triangular numbers (from Project Euler #12): divisorsLength[n_] := Apply[Times, #[[2]] + 1 & /@ FactorInteger[n]]; Needs["Lazy"] triangles = FoldList[Plus, 0, Lazy[Integers]]; triangles ~Select~ (divisorsLength[#] > 500 &) // First  Or some of the even kookier Project Euler questions: Needs["Lazy"] Rest[Lazy[Integers]]~Take~ 9999 ~Select~ ((Total@Most@Divisors@Total@Most@Divisors[#] === #) &) ~Select~ ((Total@Most@Divisors[#] =!= #) &) // Total  - I hadn't seen this post until today. I haven't tested it yet but I commend the effort! – Mr.Wizard Mar 18 '13 at 12:56 How about pumping with things like Subsets or Permutations? – Sjoerd C. de Vries Apr 25 '13 at 5:26 Not the "lazy" but shortest: Select[ToString /@ Prime[Range[10^4]], # == StringReverse[#] &]  Faster: Select[Prime[Range[10^4]], IntegerDigits[#] == Reverse[IntegerDigits[#]] &]  Faster: Pick[#, (# == Reverse[#]) & /@ IntegerDigits /@ #, True] &@ Prime[Range[10^4]]  Fastest: (twice faster than the rest, but slower than "lazy") Select[Pick[#,(#==Reverse[#])&/@ IntegerDigits/@ #, True]&@ Range[10^6], PrimeQ]  Result: {2, 3, 5, 7, 11, 101, 131, 151, 181, 191, 313, 353, 373, 383, 727, 757, 787, 797, 919, 929, 10301, 10501, 10601, 11311, 11411, 12421, 12721, 12821, 13331, 13831, 13931, 14341, 14741, 15451, 15551, 16061, 16361, 16561, 16661, 17471, 17971, 18181, 18481, 19391, 19891, 19991, 30103, 30203, 30403, 30703, 30803, 31013, 31513, 32323, 32423, 33533, 34543, 34843, 35053, 35153, 35353, 35753, 36263, 36563, 37273, 37573, 38083, 38183, 38783, 39293, 70207, 70507, 70607, 71317, 71917, 72227, 72727, 73037, 73237, 73637, 74047, 74747, 75557, 76367, 76667, 77377, 77477, 77977, 78487, 78787, 78887, 79397, 79697, 79997, 90709, 91019, 93139, 93239, 93739, 94049, 94349, 94649, 94849, 94949, 95959, 96269, 96469, 96769, 97379, 97579, 97879, 98389, 98689} - I really like that pipeline of functions you strung together in the Pick[] version. – sblom Jan 28 '12 at 4:08 @sblom Thank you ;-) – Vitaliy Kaurov Jan 28 '12 at 4:54 Here is an attempt to do this in a more "functional" way, without thought of efficiency: pQ = # === StringReverse@# & @ IntegerString@# &; ppf[x_?PrimeQ] /; pQ[x] := x ppf[x_] := ppf @ NextPrime @ x$IterationLimit = 1*^6;

NestList[ppf[# + 1] &, 2, 399]

-
where's the infix? –  acl Jan 28 '12 at 4:52
@acl There is infix all over. = === := + /; all infix. ;-) –  Mr.Wizard Jan 28 '12 at 4:56

A bit elaborate, but this works nicely for generating the first four hundred palindromic primes:

NestList[
Function[p, NestWhile[NextPrime, p,
(IntegerDigits[#] != Reverse[IntegerDigits[#]]) &, {2, 1}]],
2, 399]


An alternative route:

NestList[
Function[p, FixedPoint[NextPrime, p,
SameTest -> (IntegerDigits[#2] == Reverse[IntegerDigits[#2]] &)]],
2, 399]


This alternate is a bit faster than the previous snippet on my system, but you should do your own tests.

-
I'll only note that the advantage of this approach is that it does not have to internally generate a list of length greater than 400. ;) –  Guess who it is. Jan 28 '12 at 3:33
About twice as fast as Artes' solution, before my suggestion, but slower than both of Vitaliy's. –  rcollyer Jan 28 '12 at 3:40
@rcollyer: there's always the tradeoff between space and time, you know. :) –  Guess who it is. Jan 28 '12 at 3:45
You just have to have to make sure your have the optimal Lorentz contraction. –  rcollyer Jan 28 '12 at 3:51
I definitely like that there's no "guess the Range[]" going on here. If only NextPrime[] could be generalized to other sequences. :) –  sblom Jan 28 '12 at 4:06

There are many possible ways. Here is one in a functional style:

PalindromicPrimeQ[k_Integer] :=
IntegerDigits[k] == Reverse[IntegerDigits[k]] && PrimeQ[k]
Take[Select[Range[4 10^6], PalindromicPrimeQ], 400]


returns a list of length 400 like this:

{2,3,...101,131,...3558553, 3563653, 3569653}


Updating after rcollyer's suggestions I would rather follow this way :

PalindromicQ[k_Integer] :=
IntegerDigits[k] == Reverse[IntegerDigits[k]]
Take[Select[Prime[Range[300000]], PalindromicQ], 400];


Edit

It is worty to point out that Spartacus' test pQ seems to be a bit faster than PalindromicQ:

pQ = # === StringReverse@# &@IntegerString@# &;
Take[Select[Prime[Range[300000]], pQ], 400]; // Timing


Its timing : 2.152 versus 2.465 in case of PalindromicQ

-
Using Select[Range[10^5], PalindromicPrimeQ] and dropping Take gives the same results as Vitaliy's, but about 10 times as slow. But, if you select from Prime[Range[10^4]] (with Take still dropped), you get something a little slower than his Select method. –  rcollyer Jan 28 '12 at 3:43
From Spartacus's answer, there is the alternative test IntegerString[k] === StringReverse[IntegerString[k]]`. I wonder which test is quicker... –  Guess who it is. Jan 28 '12 at 4:24
I checked the both tests for my second approach and and Mr.Spartacus' test is faster, timings were respectively 2.153 to 2.48 –  Artes Jan 28 '12 at 4:58