Mofreh R. Zaghloul

Mofreh R Zaghloul

We present an efficient multi-accuracy algorithm for the computations of a set of special functions of a complex argument, z=x+iy. These functions include the complex probability function w(z), and closely related functions such as the error function erf(z), complementary error function erfc(z), imaginary error function erfi(z), scaled complementary error function, erfcx(z), the plasma dispersion function Z(z), Dawson s function Daw(z), and Fresnel integrals S(z) and C(z). Computational results from the present algorithm are compared with results from competitive algorithms and widely used software packages showing superior accuracy and efficiency of the present algorithm. In particular, the present results highlight concerns about the accuracy of evaluating such special functions using commercial packages like Mathematica and free/open source packages like the MIT-C++ package.

Mofreh R. Zaghloul, Ahmed N. Ali

We present a MATLAB function for the numerical evaluation of the Faddeyeva function w(z). The function is based on a newly developed accurate algorithm. In addition to its higher accuracy, the software provides a flexible accuracy vs efficiency trade-off through a controlling parameter that may be used to reduce accuracy and computational time and vice versa. Verification of the flexibility, reliability and superior accuracy of the algorithm is provided through comparison with standard algorithms available in other libraries and software packages.

Mofreh R. Zaghloul, Steven G. Johnson

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.

Mofreh R Zaghloul, Steven G. Johnson

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.