# 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]. Jan 28, 2012 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. Jan 28, 2012 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/…. Jan 28, 2012 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? Jan 28, 2012 at 17:38
• @LeonidShifrin, that comment seems worthy of a credited answer. You'd certainly get an upvote from at least me. Jan 28, 2012 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
][] /. {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 (* False *)
palindromicQ (* 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
][] /. {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. Jan 28, 2012 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 :) Jan 29, 2012 at 0:30
• +1 ! Did you have this on the shelf or did you write it for this question? Jan 29, 2012 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. Jan 29, 2012 at 5:24
• The syntax with sblom's package of WReach's solutions is Prime~Map~Lazy[Integers]~Select~palindromicQ~Take~30//List Apr 9, 2012 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:= Table[palPrimeClosure[], {100}]

Out= {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:= Table[palPrimeClosure[], {50}]

Out= {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:= palPrimeClosure500 = makePalindromePrimeC

Out= (r$10054 = NestWhile[NextPrime, p$10054,
With[{d = IntegerDigits[#1]}, d != Reverse[d]] &];
p$10054 = NextPrime[r$10054]; r$10054) & In:= Table[palPrimeClosure500[], {30}] Out= {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. Jan 28, 2012 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. Jan 28, 2012 at 7:14 • Wow! Nice. I should add that trick to my toolbox. Jan 28, 2012 at 10:38 • +1. I use similar techniques. I linked to some additional examples in the comment above. Jan 28, 2012 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, #[] + 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! Mar 18, 2013 at 12:56 • How about pumping with things like Subsets or Permutations? Apr 25, 2013 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. Jan 28, 2012 at 4:08 • @sblom Thank you ;-) Jan 28, 2012 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, 2012 at 4:52
• @acl There is infix all over. = === := + /; all infix. ;-) Jan 28, 2012 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. ;) Jan 28, 2012 at 3:33
• About twice as fast as Artes' solution, before my suggestion, but slower than both of Vitaliy's. Jan 28, 2012 at 3:40
• @rcollyer: there's always the tradeoff between space and time, you know. :) Jan 28, 2012 at 3:45
• You just have to have to make sure your have the optimal Lorentz contraction. Jan 28, 2012 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. :) Jan 28, 2012 at 4:06

I recently published a lazy list package on GitHub, so I might as well show it off here:

https://github.com/ssmit1986/lazyLists

Installation instructions can be found on the GitHub landing page above (or in the README.md file).

My code overloads the normal list processing functions like Map, MapIndexed, MapThread, FoldList, Pick, Cases, and Select. lazyRange[min, step] is the basic constructor for generating infinite ranges of integers. Take will return a lazyList in which the first argument is the extracted elements and the second is the tail of the list (which can be used to continue extracting more elements). I implemented it using ReplaceRepeated, since that seemed to be the fastest option I could find.

Example:

<< lazyLists;
Take[Select[Map[Prime, lazyRange[]], PalindromeQ], 400]

Out= lazyList[{<<primes>>}, <<tail>>]


Edit I'm still updating this repository with new functionalities and I'm open to suggestions. See the README file on the landing page for the update history.

• (+1) Please, describe how to install. Sep 19, 2018 at 23:50
• "Just download the repository (and unzip if necessary)." -- I did that. I was asking on behalf of the (GitHub) unwashed masses. Sep 20, 2018 at 14:29
• Ok, got it. The landing page on github should be fairly self-explanatory now, I hope :). Sep 20, 2018 at 15:41
• @AntonAntonov , nice tutorial here: Lazy lists in Mathematica ;-) Sep 20, 2018 at 16:11

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], 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], 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. Jan 28, 2012 at 3:43
• From Spartacus's answer, there is the alternative test IntegerString[k] === StringReverse[IntegerString[k]]. I wonder which test is quicker... Jan 28, 2012 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 Jan 28, 2012 at 4:58

I like WReach's answer because it shows how to compose an expression that is similar to the example in the question. I like Sal Mangano's answer because it is concise and shows how to make something that behaves like a C# enumerator. I hope to provide a little of each.

Enumerator[state_:0, increment_:(# + 1 &)] := Module[
{s = state},
(s = increment[s]) &
];

Where[predicate_][enumerator_] := Function[
NestWhile[enumerator, enumerator[], Not @* predicate]
];

TakeEnumerator[count_] := Table[#[], count] &;

PalindromicQ := PalindromeQ @* IntegerDigits;


Please note that Enumerator and Where both produce pure functions that behave like C# enumerators but TakeEnumerator does not.

With that the example in the question can be expressed as

Enumerator[0, NextPrime] // Where[PalindromicQ] // TakeEnumerator


Result:

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


Like in Sal's answer, you can assign to something that behaves as an enumerator.

palindrome = Enumerator[0, NextPrime] // Where[PalindromicQ];
palindrome // TakeEnumerator


Result:

{2, 3, 5, 7, 11, 101, 131, 151, 181, 191}


And repeat for the next 10:

palindrome // TakeEnumerator


Result:

{313, 353, 373, 383, 727, 757, 787, 797, 919, 929}


And you can quickly compose enumerators:

Enumerator[] /* (#^2 &) // TakeEnumerator


Result:

{1, 4, 9, 16, 25, 36, 49, 64, 81, 100}


or

Enumerator[] /* (#^2 &) /* Sqrt // TakeEnumerator


Result:

{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}