Efficient Multi-Precision Computation of Bessel Functions for Real Orders and Complex Arguments with Fortran Implementation -- Part I: The Modified Bessel Function of the First Kind, $I_ν(z)$

Abstract

This paper is the first in a series devoted to the development of efficient and highly accurate algorithms, with multiprecision \texttt{Fortran} implementations, for the computation of Bessel functions. In this first part, we present a \emph{novel, self-contained, efficient,} and \emph{multiprecision} algorithm for evaluating the modified Bessel function of the first kind, $I_ν(z)$. The method integrates several analytic representations of $I_ν(z)$, carefully selected to ensure both high accuracy and suitability for high-precision computation, together with optimally determined transition boundaries between computational regions. This design achieves high efficiency while fully preserving numerical accuracy. Unlike other widely used algorithms and libraries, such as AMOS, Boost, and GSL, which either reject negative orders $ν$ or rely on special-case symmetries valid only for integer orders, the present algorithm provides a stable approach for evaluating $I_ν(z)$ for arbitrary real orders, including $ν< 0$, and complex arguments $z$. The developed robust \texttt{Fortran} implementation provides support for both double and native quadruple-precision arithmetic. The availability of quadruple precision further enhances numerical stability, extends the reliable computational domain in $(ν, |z|)$ by approximately an order of magnitude in each direction, and enables accuracies exceeding 26 significant digits. This advancement substantially broadens the applicability of the method to demanding high-precision problems in science and engineering. Compared to AMOS (Algorithm~644), which is restricted to double precision, the present algorithm exhibits superior accuracy and efficiency, with benchmark tests demonstrating execution times reduced to 38--71\% of those of AMOS in double precision.