Efficient Multi-Precision Computation of Bessel Functions for Real Orders and Complex Arguments with a Fortran Implementation -- Part III: Regular Bessel Functions of the First and Second Kinds $J_ν(z)$ and $Y_ν(z)$
For computational scientists and engineers needing reliable Bessel function evaluations over a wide parameter range, this work provides a faster and more robust alternative to existing algorithms, particularly in quadruple precision.
This paper presents numerically stable and efficient multi-precision algorithms for computing Bessel functions J_ν(z) and Y_ν(z) with real orders and complex arguments, implemented in double and quadruple precision. The algorithms achieve execution times of 35-72% of Algorithm 644 in double precision and require only a few thousandths to hundredths of the time of Algorithm 912 while maintaining comparable or higher accuracy.
This paper is the final part in a series devoted to the development of numerically stable and efficient algorithms, together with multi-precision Fortran implementations, for computing Bessel functions with real orders and complex arguments. Parts~I and~II presented stable and efficient algorithms for the modified Bessel functions I_nu(z) and K_nu(z); here we treat Bessel functions J_nu(z) and Y_nu(z). The proposed algorithms support complex arguments for both positive and negative real orders and are implemented in native double and quadruple precision. Quadruple precision substantially increases dynamic range and accuracy (by roughly an order of magnitude in reliably computable |nu| and |z|), thereby extending applicability to problems requiring 20-30 digits. Comprehensive accuracy and performance comparisons are carried out against both the widely used Algorithm~644, restricted to double-precision arithmetic, and the more recent Algorithm~912, which supports both double- and quad-precision arithmetic as well as complex orders. In double precision, the present implementation consistently outperforms Algorithm~644, achieving execution times of approximately 35-67% for J_nu(z) and 44-72% for Y_nu(z), while also producing reliable results where Algorithm~644 fails. In comparison with Algorithm~912, the present algorithm achieves comparable accuracy in double-precision computations and significantly higher accuracy in quad-precision calculations. At the same time, it requires only a small fraction of the computational cost (from a few thousandths to a few hundredths) of the time taken by Algorithm~912, depending on the precision and parameter regime. Furthermore, unlike Algorithm~912, whose applicability is restricted to a limited region of the (Re(nu),z) plane, the present algorithm remains stable and accurate over the full tested domain for which reliable reference values are available.