Hello, I'm looking into implementation of gsl_sf_bessel_Jn_array(). There stat_np1 is computed for nmax+1 and stat_n for nmax, so that later the backward recurrence relation is used to obtain values for nmin through nmax-1. It seems strange to me though that nmax+1 is used at all — can't the function just go down from nmax and nmax-1 instead of nmax+1 and nmax? It seems currently it does one extra computation only to effectively throw away its result. This also becomes a problem e.g. when x==0.7 and nmax==141: the underflow error makes it fail, although if the function started from nmax instead of nmax+1, it'd succeed as it currently does for nmax==140.
Could someone clarify this to me? Regards, Ruslan
