Laguerre Polynomial overflow...

["Robert G. Brown" <rgb at phy dot duke dot edu>]

Dear GSL listpersons:

I'm using gsl-1.0-1 prebuilt in Red Hat 7.2 (and have been using various
gsl routines from roughly 0.7-something on).  Currently I'm testing an
ODE-based solution to the time-independent Schrodinger equation by
solving the O(3) simple harmonic oscillator problem, the radial
eigensolution to which is basically a laguerre polynomial times a few
other simple functions.

To compare the ODE solution to the exact solution, I've been evaluating

 ...
 status = gsl_sf_laguerre_n_e(nsho, lsho+0.5 , rsho*rsho, &lgp);
 if(status != GSL_SUCCESS) return(0.0);
 ret = sqrt(2.0*fac(nsho)/fac(lsho + nsho + 1.0))*pow(rsho,(double)lsho)*
         lgp.val*exp(-rsho*rsho*0.5)/norm;

 fprintf(stderr,"sho3d( %2d, %2d, %2d, %.5e) = %.5e\n",nsho,lsho,msho,rsho,ret!
 fprintf(stderr,"gsl_sf_laguerre_n( %2d, %2d, %.5e) = %.5e\n",nsho,lsho,rsho,lgp.val);
 return(ret);

which works fine for nsho = 0 or 1, but when I hit nsho > 2, the routine
runs for some values of lsho and rsho but then dies on an overflow
beyond a certain point:

sho3d(  2,  0,  0, 2.21000e+00) = 7.38523e-02
gsl_sf_laguerre_n(  2,  0, 4.88410e+00) = 1.59197e+00
sho3d(  2,  0,  0, 2.22000e+00) = 7.70712e-02
gsl_sf_laguerre_n(  2,  0, 4.92840e+00) = 1.69856e+00
sho3d(  2,  0,  0, 2.23000e+00) = 8.02146e-02
gsl_sf_laguerre_n(  2,  0, 4.97290e+00) = 1.80762e+00
gsl: laguerre.c:96: ERROR: overflow
Abort

The interesting thing is that (as you can see) the solution is extremely
pedestrian here -- n is two, \ell is 0, and the argument is around 5.
The returned value is order unity.  I can see no good reason for an
overflow to be generated externally, and presume that it is occuring
during the operation of e.g. the recursion relation that generates the
result internally.  I also seem to generate no GSL error condition --
the code crashes even though I try to trap the overflow and continue.

For what it is worth, the solution is dead on up to that point,
validated somewhat backwards by the nearly exact coincidence of this
solution with the one I'm obtaining by explicit ODE solution with e.g.
rk45 (also from gsl elsewhere in the code).

If nobody on this list has encountered this or can suggest any fix
beyond going into the source and figuring out where and why the overflow
occurs and fixing it, I'm happy enough to do the latter, but given the
hassle of building from source and the possibility that it has already
been fixed in a more recent version, I thought I'd ask first.

   rgb

Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu