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
