Observed large error from spherical Bessel functions at high orders

[Ziqi Fan]

Dear developers and maintainers of GNU Scientific Library,

My name is Ziqi Fan. I am a PhD candidate from University of Florida. I
have been using GSL for at least 3 years for my numerical solver project.
Recently, I met a bug in my solver and the bug should not exist from a
theoretical perspective. So I checked very carefully my own implementation
and found that the bug originated from using the routine for spherical
bessel functions of the second kind: gsl_sf_bessel_yl.

For instance, I tested the function using a large order n = 45, and got the
following result: y_45(6.975948) = -726330209582507479571265224704.000000,
err = 19103761820604644.000000, where the first term is the output and the
second term is an estimated error provided by the routine
"gsl_sf_bessel_yl_e". The error is relatively small compared to the output,
but it can easily cause divergence of a numerical algorithm, considering
the magnitude of its value.

I once implemented spherical bessel functions myself using the numerical
recipe, and I understand that large error occurs when x << n, where x is a
positive input to the function, and n is the order of the function. I am
wondering if it is possible to resolve the large absolute error at high
order of the spherical bessel functions of the second kind. If I have an
opportunity to communicate with an engineer or mathematician from you, I
may be able to help contribute a routine with a higher precision. If so, a
new routine with a higher precision can really contribute to the invention
of many new algorithms for powerful numerical solvers.

I really appreciate your effort in providing and maintaining this great
software, and I look forward to hearing from you!

Best,
Ziqi