0
$\begingroup$
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.

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – look Nov 13 '18 at 23:28
  • $\begingroup$ Reference for the quotations? $\endgroup$ – Michael E2 Nov 14 '18 at 1:03
  • $\begingroup$ You have numbered equations (3), (14), etc. What is paper or book they are from? $\endgroup$ – Michael E2 Nov 14 '18 at 2:41
  • 1
    $\begingroup$ This is cross-posed to Wolfram Community (which should have been noted in both posts). $\endgroup$ – Daniel Lichtblau Nov 14 '18 at 14:52
  • $\begingroup$ Yes, Thank you for noticing@DanielLichtblau $\endgroup$ – look Nov 14 '18 at 23:37
2
$\begingroup$

Please double-check for correctness.

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – look Nov 14 '18 at 7:00
  • $\begingroup$ 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.) $\endgroup$ – Henrik Schumacher Nov 14 '18 at 7:31
  • $\begingroup$ @look Hmmm. I guess I found the issue. I forgot to make updates to Cold. Please have a look. $\endgroup$ – Henrik Schumacher Nov 14 '18 at 20:56
  • $\begingroup$ 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. $\endgroup$ – look Nov 14 '18 at 23:51
  • $\begingroup$ Please try this: plist = {0} xlist = RandomReal[{0, 1}, 20]; Max[Abs[Outer[re, plist, xlist] - Outer[R, plist, xlist]]] $\endgroup$ – look Nov 16 '18 at 15:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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