Efficient Calculation of Modified Bessel Functions of the First Kind, $I_ν (z)$, for Real Orders and Complex Arguments: Fortran Implementation with Double and Quadruple Precision

TL;DR

This work addresses the challenge of efficiently and accurately computing the modified Bessel function $I_{\nu}(z)$ for real $\nu$ and complex $z$ by presenting a self-contained Fortran algorithm that supports double- and quadruple-precision arithmetic. The method combines a small-$z$ power-series expansion, a large-$z$ asymptotic expansion, and a carefully designed intermediate-region seed-recurrence strategy that avoids dependence on $K_{\nu}(z)$, with empirically tuned region boundaries. Accuracy is validated against Maple references across extensive test grids, showing substantial speedups over Algorithm 644 in double precision and enabling high-precision results with quad precision (relative errors on the order of $10^{-26}$ to $10^{-27}$). The approach broadens the practical domain for high-precision Bessel evaluations and establishes a framework transferable to other Bessel families, with broad impact for physics and engineering applications.

Abstract

We present an efficient self-contained algorithm for computing the modified Bessel function of the first kind $I_ν(z)$, implemented in a robust Fortran code supporting double and quadruple (quad) precision. The algorithm overcomes the limitations of Algorithm 644, which is restricted to double precision and applies overly conservative underflow and overflow thresholds, leading to failures in large parameter regions. Accuracy is validated against high-precision Maple calculations, and benchmarking shows execution time reductions to 54%-80% of Algorithm 644 (in double precision). Quad precision enhances numerical stability and broadens the domain of computations, making the implementation well suited for high-precision applications in physics and engineering. This work also provides a foundation for the development of efficient algorithms for other Bessel functions.

Figures (6)

• Figure 1: Schematic computational regions and corresponding methods used in the present algorithm. Specific numerical values of the points $S_{1}$, $S_{2}$, and $S_{3}$ are summarized in Table A in the appendix.
• Figure 2: The grid of tested points using double-precision arithmetic, along with additional points superimposed on or near the boundaries of the regions defined in the present algorithm, as shown in Fig. 1. Maple calculations for these input points fall within the range of representable floating-point real numbers in double-precision arithmetic.
• Figure 3: Colormap plots of the base-10 logarithm of the ratio of the relative error in calculating the real part (a) and imaginary part (b) of $I_{\nu}(z)$ of the dataset points tested using the present algorithm relative to Algorithm-644, with Maple calculations as the reference.
• Figure 4: The grid of tested points using quad-precision arithmetic, along with superimposed points on or near the borders of the regions in the present algorithm, as shown in Fig. 1. Maple calculations for these input points fall within the range of the minimum and maximum floating-point real numbers in quad-precision arithmetic.
• Figure 5: Colormap plots of the relative error in calculating the real part (a) and imaginary part (b) of $I_{\nu} (z)$ for the dataset points tested using the present algorithm, with Maple calculations as the reference.