Accuracy problem with LaguerreA function

Hello,

I need to use associated Laguerre polynomials of large n (but small k only) to solve some physics problems. It looks like there might be some numerical issues with LaguerreA Igor function when n gets relatively large (>50). The generated function get very "noisy". To confirm this, I have compared what gives LaguerreA with what gives the direct definition of LaguerreA (i.e. differentiating the Laguerre function). It looks like the "error" can be very large. An example is attached for n=55. It gets worse for larger n. The code I used is below.
Since I need to use fonctions with very large n, it is a problem if I use the LaguerreA Igor function. I can, of course, generate the fonction with the given code but having an accurate LaguerreA fonction would be very convenient.
I am not sure if this is a known behavior of LaguerreA or if I am doing something wrong.
Best,



function generate_laguerre(xrange,xpts,nrange) //this function generate a wave "laguerrewave" that contains Laguerre polynomials of degree up to nrange and x value up to xrange (with xpts points)
variable xrange,xpts,nrange
variable xpts_int=trunc(xpts)
variable n=trunc(nrange)

make /o /n=(xpts_int,n+1) laguerrewave=nan
setscale /i x,0,xrange,laguerrewave
make /o /n=(xpts_int) laguerre_xwave=nan
laguerre_xwave[]=p*xrange/(xpts_int-1)
laguerrewave[][]=laguerre(q,laguerre_xwave[p])


end


function generate_laguerrea_1_2()  //generate waves based on the generated "laguerrewave":
//"laguerrewave_1" is the associated Laguerre polynomial with k=1 get from differentiating "laguerrewave",
//"laguerrewave_2" is the associated Laguerre polynomial with k=2,
//"laguerrea_1" is the associated Laguerre polynomial with k=1 get with LaguerreA,
//"laguerrea_2" is the associated Laguerre polynomial with k=2,
//"laguerrea_1res" is the difference between the associated Laguerre polynomial with k=1 that is given by Laguerrea Igor function and "laguerrewave_1"
//"laguerrea_2res" is the difference between the associated Laguerre polynomial with k=2 that is given by Laguerrea Igor function and "laguerrewave_2"

wave laguerrewave=root:laguerrewave
wave laguerre_xwave=root:laguerre_xwave
variable nsize=dimsize(laguerrewave,1)
duplicate /o laguerrewave,laguerrewave_1,laguerrewave_1temp,laguerrewave_2,laguerrewave_2temp,laguerrea_1,laguerrea_2,laguerrea_1res,laguerrea_2res
redimension /n=(-1,(nsize-1)) laguerrewave_1,laguerrea_1,laguerrea_1res
redimension /n=(-1,(nsize-2)) laguerrewave_2,laguerrea_2,laguerrea_2res
differentiate /dim=0 laguerrewave /x=laguerre_xwave /D=laguerrewave_1temp //First differrentiate to get Laguerre(n,1,x)
laguerrewave_1[][]=(-1)^1*laguerrewave_1temp[p][q+1] //Definition of Laguerrea
laguerrea_1[][]=laguerrea(q,1,laguerre_xwave[p]) //Generate directly from Laguerrea Igor function for comparison
differentiate /dim=0 laguerrewave_1temp /x=laguerre_xwave /D=laguerrewave_2temp //First differrentiate to get Laguerre(n,2,x)
laguerrewave_2[][]=(-1)^2*laguerrewave_2temp[p][q+2] //Definition of Laguerrea
laguerrea_2[][]=laguerrea(q,2,laguerre_xwave[p]) //Generate directly from Laguerrea Igor function for comparison
laguerrea_1res=laguerrea_1-laguerrewave_1 //wave difference for comparison
laguerrea_2res=laguerrea_2-laguerrewave_2 //wave difference for comparison


end


When n gets large it is not ideal to compute LaguerreA() because the computation involves sums of terms that consist of factorials of n and binomials of k+n. On the other hand, Laguerre() function is computed from a recurrence relation which does not have numerical stability problems at small x. For large x and large n you need to be careful with both approaches.