Very Slow LaguerreL Function (mupadmex (MEX-file))

I have a function that generates a Laguerre-Gaussian profile, using the built in laguerreL on a meshgrid of X-Y values. I have noticed that this function is very slow and when timing it the element that seems to be taking the most time is listed as "mupadmex (MEX-file)" which is called a very large number of times (far more than once per element of my arrays).

It would appear that this mupad function is related to the symbolic math toolbox that I know the laguerreL function is designed to work with. Would anyone be able to enlighten me on what this function is doing exactly and if there is any way to bypass it when using built in functions like laguerreL for a faster raw computation?

I have attached my function as well just in case there is something else wrong with my style that is the result of such a slow computation! I need it running for arrays bigger than 1024x1024 which already take several hours to compute!

Any advice would be greatly appreciated!

function lg = LG(z,r,z0,w0,k0,phi,n,m)
%%Creating Variables: 
% z: distance from waist
% r: square array of rt(X^2 + Y^2)
% z0,w0,k0: constants
% phi: array of atan2(Y,X)
% m,n: Laguerre indices

w = w0*sqrt(1+(z/z0)^2);
R = z*(1+(z0/z)^2);
Phase = -1j*m*phi;
Gouy = 1j*(2*n+m+1)*atan(z/z0);
Arg = ((2*r.^2)/(w^2));
Norm = sqrt((2*factorial(n))/(pi*factorial((abs(m)+n))));

%%Building Laguerre-Gaussain:

A = w0/w;
B = ((sqrt(2)*r)/(w)).^abs(m);
C = exp(-1*(r.^2)/(w^2));
P = -1j*k0*((r.^2)/(2*R));

lg = Norm*A*B.*C...
    .*laguerreL(n,abs(m),Arg)...
    .*exp(P+Phase+Gouy);
end

Answer 1

[Walter Roberson]

laguerreL is part of the Symbolic Toolbox, the implementation of which is the library of code accessed by calling mupadmex(). We are not told whether it is implemented as symbolic MuPAD routines, or if there is a special case for Digits at most 16 that allows it to call pure numeric routines: that is an implementation detail. You should not understand it as being a pure numeric routine: you should understand it as a software mathematics routine that uses whatever software an internal precision is necessary to achieve the mathematics goals.

Depending on the arguments you are passing, you might be able to convert it into a combination of binomial() and hypergeom(). Both of those are implemented in MuPAD, but for at least some of the cases, they could be implemented numerically. For example see http://www.mathworks.com/matlabcentral/fileexchange/5616-generalized-hypergeometric-function

LaguerreL(a, b, z) = binomial(a + b, a) hypergeom([-a], [b + 1], z)
              provided b ~= -1 and b + a ~= -1