I don’t know of any libraries. I can see that on p438 (10.1.10) of Abrimowitz and Stegun there is a recursion formula for Spherical Bessel Functions in using trigonometric functions which seems to be very simple.
My thought is not so much that you use that to evaluate the function directly since there could be a serious performance hit doing that, rather that you use might be able to use it in a script to generate the Fortran code for each order n. Might be an interesting case for using parameter packs in C++
I have a code that computes j_n ( m chi ) and y_n ( chi ) for double precision chi and complex m. It’s based upon a method of Infeld and a routine DBMIE by J. V. Dave (1968) – used for Mie scattering – and adapted with extensions by C. D. Cantrell (1988). I upgraded it to Fortran 2000 in 2007. I made a simplistic (not comprehensive) comparison to reference codes for arbitrary order Bessel and Neumann functions. The maximum absolute error for j_n (x) with log10(x) = (-1(0.2)2), n = 1…100 and chi=1.0 (i.e., real j_n only) was 1.65e-15. The maximum relative error for y_n was 1.33e-12.
Wrt the library by Zhang and JIn: I had some problems in the past with evaluating the modified Bessel function K (never sure of the exact denomination, I use them too seldom). The routine did not converge properly for some range of the argument and came up with Nan, if I remember correctly. That does not mean that the library is no use, of course, far from that.
FWIW, I found code for the complex gamma function that seems fairly accurate. I checked a few data points using the relation gamma(x) = (x-1)*gamma(x-1). The relative errors were in the order of 1.0e-15.
(Not a function I really need, but I was curious as the Fortran standard does not define the gamma function for complex arguments)
j_n is unconditionally unstable for forward recursion. y_n is unconditionally unstable for backward recursion. So if you compute Hankel functions by forward (or reverse) recursion, you get (y_n, y_n) or (j_n,j_n).
No one that knows what they are doing will compute Hankel functions with strictly forward or backwards recursion. Fortunately, the Hankel function of first kind is h_n = cmplx(j_n,y_n), and so is trivial to compute.
Also note the claim of unconditionally unstable is a bit strong. For example, one can use forward recursion for j_n provided n < x. Why would one want to use forward recursion if possible, well speed of course.
integer, parameter :: n = 20, nmax = 10
real jfwd(0:n), jbwd(0:n), x, ix, jmx1, jmx, j0, j1
x = 10
ix = 1 / x
j0 = sin(x) * ix
j1 = (j0 - cos(x)) * ix
! Forward recursion.
jfwd(0) = j0
jfwd(1) = j1
do i = 1, n - 1
jfwd(i+1) = (2 * i + 1) * ix * jfwd(i) - jfwd(i-1)
! Backward recursion.
jmx1 = 0
jmx = 1
do i = n + nmax, n + 2, -1
jbwd(0) = (2 * i + 1) * (ix * jmx - jmx1 / (2 * i + 1))
jmx1 = jmx
jmx = jbwd(0)
do i = n + 1, 1, -1
jbwd(i - 1) = (2 * i + 1) * ix * jmx - jmx1
jmx1 = jmx
jmx = jbwd(i - 1)
if (abs(jmx) >= abs(jmx1)) then
j0 = j0 / jbwd(0)
j0 = j1 / jbwd(1)
jbwd = j0 * jbwd
! j(0) = -5.44021110889369813405E-2
! j(1) = 7.84669417987515470918E-2
! j(9) = 1.00096409548490576567E-2
! j(10) = 6.46051544925642642714E-2
! j(11) = 3.55744148858943784028E-2
! j(20) = 2.30837196131946871671E-6
do i = 0, n
write(*,'(I2, 2ES15.6)') i, jfwd(i), jbwd(i)
end program van
I leave it as an exercise to others to copy-n-paste to test.