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]


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

This is another answer from stackexchange.com

I am afraid I cannot agree with you. It will spend more time

May I ask what is the further step to optimize this code so that the runtime will be less than 1 sec?

Your original code never completes (I killed it after 5 minutes). I had to protect the evaluation of re[] for values of r that make it complex by adding an if statement.

If[r <= 1, re[pmax, r], 0]

If you make the table limits floating point, (delta/2. vs delta/2) it takes about 24.5 seconds.

You have several issues.

1. You recompute Sqrt[8/Pi] many times -- do it once. 2.This code takes 4.9 seconds or a speedup of 5 times faster.

   Clear[r,re,p,pmax,delta,imagesize,delta]
   ClearSystemCache[];
   sqrt8 = Sqrt[8/Pi];
   cf1=Compile[{{rr,_Real}},sqrt8*((1-rr)/rr)^(1/4)*1];
   re[0,r_]:=cf1[r];
   cf2=Compile[{{rr,_Real}},sqrt8*((1-rr)/rr)^(1/4)*-1*2*(1-2*rr)];
   re[1,r_]:=cf2[r];
   cf3 = Compile[{{pp,_Integer},{rr,_Real}},sqrt8*((1-rr)/rr)^(1/4)*(-1)^pp];
   re[p_,r_]:=re[p,r]=cf3[p,r]*(re[1,r]*re[p-1,r]-re[p-2,r]);
   imagesize=256;
   pmax=120;
   delta=2/imagesize;
   {foo,bar}=Timing[Table[r=Sqrt[x^2+y^2];If[r<=1,re[pmax,r],0],{x,-1+delta/2.,1-delta/2.,delta},{y,1-delta/2.,-1+delta/2.,-delta}]];
   

There are some further speedups like maybe not looping through values of r that make no sense but this is a start.

Regards,

Neil

You have a bug. r is defined to be Sqrt[x^2+y^2] but also 0.1? I think you need to fix some notation. I also would not use r to be the pattern variable in your function definitions. (Three uses of the same variable name can cause strange results). Can you fix the code and repost it in complete form?

Regards,

Neil

Your first problem is that you are not computing the expressions numerically. They have exact integer values so MMA is keeping these huge expressions in memory. Add N[] around your expressions so they become numerical. This will give a huge speed up. Second, your results are coming out complex. Is this correct? I have some other suggestions but I wanted to verify if the expressions are correct before diving deeper.

Regards,

Neil