# how to improve the performance of this recursion code in mathematica

Clear[r, re, p, pmax, delta, imagesize, delta]
ClearSystemCache[]
re[0, r_] := Sqrt[8/Pi]*((1 - r)/r)^(1/4)*1;
re[1, r_] := Sqrt[8/Pi]*((1 - r)/r)^(1/4)*-1*2*(1 - 2*r);
re[p_, r_] := re[p, r] = Sqrt[8/Pi]*((1 - r)/r)^(1/4)*(-1)^p*(re[1, r]*re[p - 1, r] - re[p - 2, r]);
imagesize = 32;
pmax = 10;
delta = 2/imagesize;
Table[r = Sqrt[x^2 + y^2]; re[pmax, r], {x, -1 + delta/2, 1 - delta/2, delta}, {y, 1 - delta/2, -1 + delta/2, -delta}];


when the imagesize and pmax increase, the time will become unacceptable. So, I would ask if we can use compile of other methods to speed up,like: for imagesize is 256 and pmax is 120, the time will be about 10 seconds. In my code, I also use the memoization to store the value during the evaluation which I will use in the future.

• You might try this Simplify[RSolve[{re[0] == c, re[1] == c*-1*2*(1 - 2*r), re[p] == c*(-1)^p*((c*-1*2*(1 - 2*r))*re[p - 1] - re[p - 2])}, re[p], p]] /. c -> Sqrt[8/Pi]*((1 - r)/r)^(1/4) There are lots of common subexpressions in that result which you could calculate once, store in a variable and the replace all those subexpressions with simple variables. If you do all that carefully and simplify the result you might end up with a reasonable sized expression with a few variables and no recursion. – @Bill
– look
Nov 13, 2018 at 23:28
• Reference for the quotations? Nov 14, 2018 at 1:03
• You have numbered equations (3), (14), etc. What is paper or book they are from? Nov 14, 2018 at 2:41
• This is cross-posed to Wolfram Community (which should have been noted in both posts). Nov 14, 2018 at 14:52
• Yes, Thank you for noticing@DanielLichtblau
– look
Nov 14, 2018 at 23:37

re = Compile[{{pmax, _Integer}, {r, _Real}},
Block[{C1, Cold, Cp, buffer},
If[0. <= r <= 1.,
If[pmax == 0, Cp = 1., C1 = 2. (1. - 2. r);
Cold = 1.;
Cp = C1;
Do[buffer = Cp;
Cp = C1 Cp - Cold;
Cold = buffer, {p, 2, pmax}];
];
Sqrt[8./Pi] ((1 - r)/r)^(1/4) (-1.)^pmax Cp,
0.
]
],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True]


The sum expression for testing purposes:

R[p_, r_] := Sqrt[8./Pi] ((1 - r)/r)^(1/4) Sum[
(-1)^k (p - k)!/k!/(p - 2 k)! (2 (2 r - 1))^(p - 2 k),
{k, 0, Quotient[p, 2]}]


Usage example:

imagesize = 256;
delta = 2./imagesize;
x = Subdivide[-1. + delta/2, 1. - delta/2, imagesize - 1];
y = Subdivide[-1. + delta/2, 1. - delta/2, imagesize - 1];
result = re[120, Sqrt[0.5 Outer[Plus, x^2, y^2]]]; // AbsoluteTiming // First


0.006015

Notice that the function re will throw errors if you try to feed it with second arguments from outside the unit interval.

Check for correctness:

plist = RandomInteger[{1, 20}, 20];
xlist = RandomReal[{0, 1}, 20];
Max[Abs[Outer[re, plist, xlist] - Outer[R, plist, xlist]]]


1.2568*10^-9

• it seems the result is not right, I checked the code:"Sqrt[8./Pi] ((1 - r)/r)^(1/4) (-1.)^pmax Cp". According to the equation (14), the p will change, not the only pmax. And thank you very much.
– look
Nov 14, 2018 at 7:00
• Hm. I don't think that this makes sense. For given $p_\mathrm{max}$, we want to compute $R_{p_\mathrm{max}}(r)$; therefore, we need $C_{p_\mathrm{max}}(r)$. And the latter is computed recursively. (Notice that I did it iteratively; this is in general faster than doing it recursively.) Nov 14, 2018 at 7:31
• @look Hmmm. I guess I found the issue. I forgot to make updates to Cold. Please have a look. Nov 14, 2018 at 20:56
• Thanks for updating, It seems the result is different from another answer. according to my result from my own method, I think his method is right so far. But your performance is much better.
– look
Nov 14, 2018 at 23:51
• Please try this: plist = {0} xlist = RandomReal[{0, 1}, 20]; Max[Abs[Outer[re, plist, xlist] - Outer[R, plist, xlist]]]
– look
Nov 16, 2018 at 15:54