Dynamic programming is a technique for avoiding the repeated computation of the same values in a recursive program. Each value computed is immediately stored. If the value is needed again, it is not computed but simply looked up in the table. (1)

I use orthogonal polynomials a fair bit in my work. Since Mathematica supports only the classical ones, I often have to write my own functions. For instance, the monic Charlier polynomials satisfy the three-term recurrence

$$C_{n+1}^{(a)}(x)=(x-a-n)C_n^{(a)}(x)-an C_{n-1}^{(a)}(x)$$

with $C_0^{(a)}(x)=1$ and $C_1^{(a)}(x)=x-a$.

If I want to be able to use monic Charlier polynomials in Mathematica, I can do this:

CharlierC[0, a_, x_] := 1;
CharlierC[1, a_, x_] := x - a;
CharlierC[n_Integer, a_, x_] := (x - a - n + 1) CharlierC[n - 1, a, x] -
                                 a (n - 1) CharlierC[n - 2, a, x]

The problem with this route, of course, is that the effort expended to generate, say, CharlierC[20, a, x] can't be used for evaluating CharlierC[50, a, x]. For a one-argument recursive function (e.g. Fibonacci), dynamic programming is fine and dandy for saving evaluation effort. For a multiple-argument function, imagine what would happen if one had used the definition CharlierC[n_Integer, a_, x_] := CharlierC[n, a, x] = (* stuff *) and then executed Plot[{CharlierC[5, 1, x], CharlierC[6, x, 2]}, {x, -1, 1}].

Is there a way to reap the benefits of dynamic programming on a multiple-argument function, while storing only results where the recursion variable (n in the Charlier example) changes?

  • $\begingroup$ As I am sure you know CharlierC[n, x, a] == (-1)^n HypergeometricU[-n, 1 - n + x, a], so there is no need to use the recursive definition for this one. $\endgroup$
    – Sasha
    Commented Jan 18, 2012 at 18:46
  • $\begingroup$ @Sasha: I certainly do, but I also deal with even more complicated polynomials than Charlier polynomials; I just wanted a simple example of an unsupported polynomial system. $\endgroup$ Commented Jan 18, 2012 at 21:47
  • $\begingroup$ Out of curiosity, is this related to probability anyhow? Like I just learned that Charlier polynomials are orthogonal with respect to Poisson measure, which makes them well suited for developing distribution-specific quadrature rules $\endgroup$
    – Sasha
    Commented Jan 18, 2012 at 21:53
  • $\begingroup$ @Sasha: I'm actually working with Askey-Wilson class polynomials, as well as $q$-analogs of classical orthogonal polynomials, for one of my personal projects. I know about being able to use QHypergeometricPFQ[], but it seems that even a recursive definition can sometimes go a bit faster than using that. But yeah, even the humble Charlier polynomial is fascinating! :) $\endgroup$ Commented Jan 18, 2012 at 22:00
  • $\begingroup$ Can you give an update of any of the functions available in versions up to 11.1 provide more orthogonal polynomial other that the classic ones? $\endgroup$ Commented Apr 27, 2017 at 13:16

5 Answers 5


Yes, there is, although the speed-up is not as dramatic as for 1D memoization:

CharlierC[0, a_, x_] := 1;
CharlierC[1, a_, x_] := x - a;
CharlierC[n_Integer, a_, x_] :=
  Module[{al, xl},
    Set @@ Hold[CharlierC[n, al_, xl_],          
        Expand[(xl - al - n + 1) CharlierC[n - 1, al, xl] - 
             al (n - 1) CharlierC[n - 2, al, xl]
    CharlierC[n, a, x]

(Thanks to @Mike Bantegui for pointing out the wastefulness of Simplify, which has now been removed). What you memoize here are function definitions.Expand is used to not accumulate the complexity too fast. The idea is that I first create a new pattern-based definition, using a number of tricks to fool the scoping variable - renaming mechanism but localize pattern variables, and then evaluate this definition.

For example:

In[249]:= CharlierC[20,a,x];//Timing
Out[249]= {0.063,Null}

In[250]:= CharlierC[25,a,x];//Timing
Out[250]= {0.078,Null}

While with clear definitions:

In[260]:= CharlierC[25,a,x];//Timing
Out[260]= {0.094,Null}

Here are a first few generated definitions:

In[262]:= Take[DownValues[CharlierC],4]
      al$4106^2-xl$4106-2 al$4106 xl$4106+xl$4106^2,
   -al$4105^3+2 xl$4105+3 al$4105 xl$4105+3 al$4105^2 xl$4105
       -3 xl$4105^2-3 al$4105 xl$4105^2+xl$4105^3}
  • $\begingroup$ I timed both versions, and yours is significantly slower on my machine. 0.016 seconds (mine) vs 1.185 seconds (yours) for n = 20. Still, very cool idea you have. $\endgroup$ Commented Jan 17, 2012 at 23:59
  • 1
    $\begingroup$ @Mike My code creates definitions for absolutely arbitrary a. Also, some time is spent on Expand and Simplify. It is slow the first run, but instant all the following runs. If your task is to vary a (say, plot when a is in some range), I only need to compute stuff once, and it will work for every value of a. If there is a single numerical a, then your code wins, but then arguably a is a constant rather than a true function argument, so the memoization is effectively 1D. $\endgroup$ Commented Jan 18, 2012 at 0:04
  • $\begingroup$ @LeonidShifrin: If you want to vary a and x, you can simply evaluate it and attach it to a variable: poly = CharlierC[100, a, x] and then evaluate it for any arbitrary {a, x}. You just have to be careful to only evaluate CharlierC[n, a, x] for sybmolic a and x for my version, and always attach it to something beforehand. $\endgroup$ Commented Jan 18, 2012 at 0:05
  • $\begingroup$ I also deal with polynomials with more than three arguments, like the Hahn polynomials, and I expect to be varying indices and parameters a lot. So, Leonid's proposal seems to be most expedient for me. $\endgroup$ Commented Jan 18, 2012 at 0:14
  • $\begingroup$ @Mike Good point, I will edit. Simplify was taking the most time, I should have tried Expand alone. As to your first point: since in my case you compute all the previous ones in the process, you need almost no time to compute say CharlierC[101, a, x] fully symbolically then. $\endgroup$ Commented Jan 18, 2012 at 0:19

Nice question. This is my suggested implementation. Evaluate all code at once.

Clear[CharlierC, "CharlierC`*"]

CharlierC (* create symbol in current context *)


implementation[0] := 1;
implementation[1] := x - a;
implementation[n_Integer] := 
 implementation[n] = Expand[(x - a - n + 1) implementation[n - 1] - 
                            a (n - 1) implementation[n - 2] ]

CharlierC[n_Integer ? NonNegative, ai_, xi_] := 
 Block[{a, x}, implementation[n] /. {a -> ai, x -> xi}]


The Block isn't strictly necessary as a and x are already created in the CharlierC` context where they should be safe.

This function memoizes the symbolic representation of CharlierC only, and only for each n (not a and x). Then substitutes in whatever variables or numbers you pass in.

  • $\begingroup$ Now let me look at the other answers which I didn't check on purpose before posting (so the fun won't be spoiled!) $\endgroup$
    – Szabolcs
    Commented Jan 18, 2012 at 16:49
  • 2
    $\begingroup$ +1. This is simpler than what I suggested. Andy Ross did the first step just like you, but did not do replacements and deleted his answer, which is a pity. I only later realized that his (and now yours) code is no less general than mine, and perhaps also faster. $\endgroup$ Commented Jan 18, 2012 at 16:51
  • $\begingroup$ @Leonid My code is much slower than yours, but it is not completely clear to me why. There was also a new, pure function based solution posted today. Now it seems so obvious that pure functions are the simplest solution here ... $\endgroup$
    – Szabolcs
    Commented Feb 6, 2012 at 10:13
  • 2
    $\begingroup$ @Leonid My code was so slow because of the wrong placing of Expand, corrected now. $\endgroup$
    – Szabolcs
    Commented Feb 6, 2012 at 10:15
  • $\begingroup$ I see. That new solution is nice indeed. It did not cross my mind to use pure functions here. $\endgroup$ Commented Feb 6, 2012 at 11:58

I would also suggest to use pure functions here:

CharlierC[0] = 1 &;
CharlierC[1] = #2 - #1 &;
CharlierC[n_Integer] := (CharlierC[n] = 
      Expand[(#2 - #1 - n + 1) CharlierC[n - 1][#1, #2] - #1 (n - 
           1) CharlierC[n - 2][#1, #2]]] &);
CharlierC[20][a, x] // AbsoluteTiming

(* ==> {0.0312414, 
 a^20 - 121645100408832000 x - 128047474114560000 a x - 
  67580611338240000 a^2 x - 23851980472320000 a^3 x - 
  6335682312960000 a^4 x - 1351612226764800 a^5 x - 
  241359326208000 a^6 x - 37132204032000 a^7 x - 
  5028319296000 a^8 x - 609493248000 a^9 x - 67044257280 a^10 x - 
  6772147200 a^11 x - 634888800 a^12 x - ...

I haven't made comparisons, but this seems to be faster.

  • $\begingroup$ The disadvantage with this solution is that it memoizes for the second argument as well (a). If we evaluate CharlierC[5,a], the cached intermediate results can't be used by CharlierC[5,b]. I think it should be possible to fix this. Can you update the solution? $\endgroup$
    – Szabolcs
    Commented Feb 6, 2012 at 9:59
  • $\begingroup$ Done :-) Is it okay now? $\endgroup$
    – faleichik
    Commented Feb 6, 2012 at 10:08
  • $\begingroup$ +1, great solution! But to nitpick a bit more, Evaluate is not needed in CharlierC[1] = Evaluate[#2 - #1] & I placed Expand at the wrong position in my solution which caused a serious performance problem. I corrected it now. $\endgroup$
    – Szabolcs
    Commented Feb 6, 2012 at 10:18
  • $\begingroup$ Yes, I've forgotten to remove this from the previous version. $\endgroup$
    – faleichik
    Commented Feb 6, 2012 at 10:23
  • $\begingroup$ Nice code - +1. $\endgroup$ Commented Feb 6, 2012 at 11:57

You can utilize Dynamic Programming in the following way:

CharlierC[0, a_, x_] := 1
CharlierC[1, a_, x_] := x - a
CharlierC[n_Integer, a_, x_] := CharlierC[n, a, x] = Expand[Expand[(x - a - n + 1) CharlierC[n - 1, a, x]] - Expand[a (n - 1) CharlieC[n - 2, a, x]]]

It basically creates an in-memory store of previous evaluated function values instead of evaluating it every time you call it.

Mind you, if you actually tried evaluating CharlierC[50, a, x] it would still take a very long time. That's a very long polynomial you're trying to evaluate.


The addition of the Expand[..] forces Mathematica to expand and combine like terms. It's actually possible to evaluate CharlierC[50, a, x].

If you wanted to really use that expression, here's how I would do it:

poly = CharlierC[50, a, x];

a = 0.01;
Plot[poly, {x, ...}];

Here's some timings. In between each run I did a ClearAll and evaluated the definition:

AbsoluteTiming[CharlierC[20, a, x];] (* 0.0250014 *)
AbsoluteTiming[CharlierC[40, a, x];] (* 0.1510087 *)
AbsoluteTiming[CharlierC[60, a, x];] (* 0.5150295 *)

What I would recommend is combining the ideas that Leonid and I have:

CharlierC[0, a_, x_] := 1;
CharlierC[1, a_, x_] := x - a;
CharlierC[n_Integer, a_, x_] := Module[{al, xl},
   Set @@ Hold[CharlierC[n, al_, xl_],
      Expand[(xl - al - n + 1) CharlierC[n - 1, al, xl]] - 
       Expand[al (n - 1) CharlierC[n - 2, al, xl]]
   CharlierC[n, a, x]];

The use of Expand works much faster than Simplify.

  • $\begingroup$ Okay, try executing CharlierC[10, 3, 1] and subsequently ?CharlierC. Messy, no? Using that version of CharlierC[] within something like Plot[] will make things worse. $\endgroup$ Commented Jan 17, 2012 at 23:50
  • 1
    $\begingroup$ I did actually use Expand in my initial answer as well - it was Simplify which was slow. And, you only need a single top-level one I think - I also just made an edit. But, again, thanks for your point, and +1. $\endgroup$ Commented Jan 18, 2012 at 0:30

Decided to add one more version: what happened is that my question How to memoize with patterns? has been marked as a duplicate of this one, which I agree to. However I believe there are more general situations where it should work. I believe my version covers them.

toPatterns[x_] := 
 x /. Map[# -> Pattern[#, _] &, 
   Union@Select[Level[x, {-1}], Head[#] === Symbol &]]
myCharlierC[0, a_, x_] := 1;
myCharlierC[1, a_, x_] := x - a;
myCharlierC[n_Integer, a_, x_] := 
  myCharlierC[n, toPatterns[a], toPatterns[x]] = 
   Expand[(x - a - n + 1) myCharlierC[n - 1, a, 
       x] - (n - 1) a  myCharlierC[n - 2, a, x]];

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.