Functional style using lazy lists?

111

73

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.

sblom

Posted 2012-01-28T02:33:27.567

Reputation: 6 273

@sblom your package is an interesting work. Some usage examples would further increase its value in my opinion. – faysou – 2012-04-09T11:54:05.993

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]. – J. M.'s ennui – 2012-01-28T03:46:15.190

1Seems 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 – 2012-01-28T04:04:35.570

3

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 - http://stackoverflow.com/questions/8596134/efficient-alternative-to-outer-on-sparse-arrays-in-mathematica/8609071#8609071, http://mathematica.stackexchange.com/questions/559/what-are-the-use-cases-for-different-scoping-constructs/569#569 (Module -> Advanced uses), more examples and links in my third post here: http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/5eff6213a0ab5f51.

– Leonid Shifrin – 2012-01-28T11:11:33.823

@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 – 2012-01-28T17:38:32.057

@LeonidShifrin, that comment seems worthy of a credited answer. You'd certainly get an upvote from at least me. – sblom – 2012-01-28T17:39:52.333

The technique behind it was described in the aswer by @Sal already, so I did not think my additions would merit a separate answer. Thanks anyways, I am happy if that was helpful. – Leonid Shifrin – 2012-01-28T17:51:51.747

What is used in this question is very similar. Just pointing to a practical use on this very site.

– Szabolcs – 2012-01-28T17:53:16.720

@sblom: it's a bit faster than Append[]/AppendTo[] when I checked. You might want to experiment on your setup to compare. – J. M.'s ennui – 2012-01-28T18:09:14.600

Answers

98

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[
      NestWhile[(Sow[sHead[#], tag]; sTail[#])&, s, !sEmptyQ[#]&]
    , 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;

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

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[
      NestWhile[(Sow[sHead[#], tag]; sTail[#])&, s, !sEmptyQ[#]&]
    , 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]

WReach

Posted 2012-01-28T02:33:27.567

Reputation: 62 787

3The syntax with sblom's package of WReach's solutions is Prime~Map~Lazy[Integers]~Select~palindromicQ~Take~30//List – faysou – 2012-04-09T20:20:39.270

9I 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 – 2012-01-28T23:37:12.090

3

@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 – 2012-01-29T00:30:28.237

2+1 ! Did you have this on the shelf or did you write it for this question? – Mr.Wizard – 2012-01-29T05:14:03.413

@Spartacus I wrote it for the question, although it borrows concepts liberally from work I have done before in other languages. – WReach – 2012-01-29T05:21:34.723

3Color 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 – 2012-01-29T05:24:21.337

I am total newbie in this area but I love to program with Pythonic list comprehensions. Is it easy to remix this style to create a style for functional Python programmers to get more involved in Mathematica? I moved my curiosity here and I would very thankful if you overviewed it.

– hhh – 2015-06-27T22:22:13.580

43

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}

Sal Mangano

Posted 2012-01-28T02:33:27.567

Reputation: 1 325

It is indeed much faster then a blunt approach scanning all the numbers. – Vitaliy Kaurov – 2012-01-28T06:47:10.250

2That'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 – 2012-01-28T07:14:51.787

Wow! Nice. I should add that trick to my toolbox. – kkm – 2012-01-28T10:38:04.350

+1. I use similar techniques. I linked to some additional examples in the comment above. – Leonid Shifrin – 2012-01-28T11:18:08.870

31

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

sblom

Posted 2012-01-28T02:33:27.567

Reputation: 6 273

I hadn't seen this post until today. I haven't tested it yet but I commend the effort! – Mr.Wizard – 2013-03-18T12:56:16.013

How about pumping with things like Subsets or Permutations? – Sjoerd C. de Vries – 2013-04-25T05:26:26.707

22

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}

Vitaliy Kaurov

Posted 2012-01-28T02:33:27.567

Reputation: 66 672

4I really like that pipeline of functions you strung together in the Pick[] version. – sblom – 2012-01-28T04:08:04.957

@sblom Thank you ;-) – Vitaliy Kaurov – 2012-01-28T04:54:20.807

20

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]

Mr.Wizard

Posted 2012-01-28T02:33:27.567

Reputation: 259 163

where's the infix? – acl – 2012-01-28T04:52:50.307

1@acl There is infix all over. = === := + /; all infix. ;-) – Mr.Wizard – 2012-01-28T04:56:26.717

13

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.

J. M.'s ennui

Posted 2012-01-28T02:33:27.567

Reputation: 115 520

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. ;) – J. M.'s ennui – 2012-01-28T03:33:57.750

About twice as fast as Artes' solution, before my suggestion, but slower than both of Vitaliy's.

– rcollyer – 2012-01-28T03:40:21.910

@rcollyer: there's always the tradeoff between space and time, you know. :) – J. M.'s ennui – 2012-01-28T03:45:40.387

You just have to have to make sure your have the optimal Lorentz contraction. – rcollyer – 2012-01-28T03:51:41.740

@rcollyer Note that Vitaliy's both approaches are slower than mine if they have to return the first 400 palindromic primes. J.M's method as I pointed out is faster (since doesn't have to know in advance the length of range of integers to select) but NestWhile makes resemblance of procedural approach. So one has to decide what he prefers – Artes – 2012-01-28T03:57:03.997

@Artes, I had missed the point about it being the first 400, and hadn't tested it out that far. But, for the shorter list (113 primes), Vitaliy's was faster. The addition of Prime supercharged your method and put it in the same league. But, required some knowledge of where to put the cutoff. – rcollyer – 2012-01-28T04:06:07.730

3I definitely like that there's no "guess the Range[]" going on here. If only NextPrime[] could be generalized to other sequences. :) – sblom – 2012-01-28T04:06:44.210

I'm not sure if Vitaliy's approach is slower, but his method does't return exactly the first 400 palindromic primes. – Artes – 2012-01-28T04:24:18.037

@Artes to be clear: slower with a shorter list. Yours has a much slower growth curve and for 400 palindromic primes, yours beats him. – rcollyer – 2012-01-28T04:29:02.557

@rcollyer Yes, I had to make a correction, but I believe it is worthy if one makes extensively such seeking tasks. J.M.'s aproach and Mr.Spartacus' one are slower, although elegant. – Artes – 2012-01-28T04:43:54.940

12

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[138]= 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.

Sjoerd Smit

Posted 2012-01-28T02:33:27.567

Reputation: 15 418

(+1) Please, describe how to install. – Anton Antonov – 2018-09-19T23:50:05.190

"Just download the repository (and unzip if necessary)." -- I did that. I was asking on behalf of the (GitHub) unwashed masses. – Anton Antonov – 2018-09-20T14:29:52.270

2Ok, got it. The landing page on github should be fairly self-explanatory now, I hope :). – Sjoerd Smit – 2018-09-20T15:41:38.533

2

@AntonAntonov , nice tutorial here: Lazy lists in Mathematica ;-)

– Vitaliy Kaurov – 2018-09-20T16:11:09.897

9

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

Artes

Posted 2012-01-28T02:33:27.567

Reputation: 51 831

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 – 2012-01-28T03:43:28.850

From Spartacus's answer, there is the alternative test IntegerString[k] === StringReverse[IntegerString[k]]. I wonder which test is quicker... – J. M.'s ennui – 2012-01-28T04:24:43.793

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 – 2012-01-28T04:58:34.207

5

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[400]

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[10]

Result:

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

And repeat for the next 10:

palindrome // TakeEnumerator[10]

Result:

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

And you can quickly compose enumerators:

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

Result:

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

or

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

Result:

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

Rik Renich

Posted 2012-01-28T02:33:27.567

Reputation: 151