Double Precision Arithmetic Research Articles

Algorithms for numerically stable scattering amplitudes

The numerically stable evaluation of scattering matrix elements near the infrared limit of gauge theories is of great importance for the success of collider physics experiments. We present a novel algorithm that utilizes double-precision arithmetic and reaches higher precision than a naive quadruple-precision implementation at smaller computational cost. The method is based on physics-driven modifications to propagators, vertices, and external polarizations. Published by the American Physical Society 2024

Computing the Mittag-Leffler function of a matrix argument

It is well-known that the two-parameter Mittag-Leffler (ML) function plays a key role in Fractional Calculus. In this paper, we address the problem of computing this function, when its argument is a square matrix. Effective methods for solving this problem involve the computation of higher order derivatives or require the use of mixed precision arithmetic. In this paper, we provide an alternative method that is derivative-free and works entirely using IEEE standard double precision arithmetic. If certain conditions are satisfied, our method uses a Taylor series representation for the ML function; if not, it switches to a Schur-Parlett technique that will be combined with the Cauchy integral formula. A detailed discussion on the choice of a convenient contour is included. Theoretical and numerical issues regarding the performance of the proposed algorithm are discussed. A set of numerical experiments shows that our novel approach is competitive with the state-of-the-art method for IEEE double precision arithmetic, in terms of accuracy and CPU time. For matrices whose Schur decomposition has large blocks with clustered eigenvalues, our method far outperforms the other. Since our method does not require the efficient computation of higher order derivatives, it has the additional advantage of being easily extended to other matrix functions (e.g., special functions).

Nine-Stage Runge–Kutta–Nyström Pairs Sharing Orders Eight and Six

We explore second-order systems of non-stiff initial-value problems (IVPs), particularly those cases where the first derivatives are absent. These types of problems are of significant interest and have applications in various domains, such as astronomy and physics. Runge–Kutta–Nyström (RKN) pairs stand out as highly effective methods of addressing these IVPs. In order to create a pair with eighth and sixth orders, we need to address a certain known set of equations concerning the coefficients. When constructing such pairs for use in double-precision arithmetic, we often need to meet various conditions. Primarily, we aim to maintain small coefficient magnitudes to prevent a loss of accuracy. Nevertheless, in the context of quadruple precision, we can tolerate larger coefficients. This flexibility enables us to establish pairs with eighth and sixth orders that exhibit significantly reduced truncation errors. Traditionally, these pairs are constructed to go through eight stages per step. Here, we propose using nine stages per step. Then we have available more coefficients in order to further reduce truncation errors. As a result, we construct a novel pair that, as anticipated, achieves superior performance compared to equivalent-order pairs in various significant problem scenarios.

Calculating Correctly Rounded Exponential Function in Double-Precision Using Extended Double-Precision Arithmetic

The article presents an effective algorithm for calculating the correctly rounded exponent for any of the rounding modes specified by the IEEE 754 standard. The argument and the function value are double-precision numbers, but the algorithm itself uses extended double-precision arithmetic for calculations. At the argument reduction stage, we approximate the argument with numbers whose exponent value has a short mantissa, which makes calculations faster. The article gives a formal description of the algorithm and a proof of its correctness. The function has the shortest maximum execution time among the considered functions calculating correctly rounded exponent. The execution time slightly depends on the value of the argument, which may be important for critical applications.

Calculation of Fresnel integrals of real and complex arguments up to 28 significant digits

We present numerical algorithms for efficient and highly accurate computations of Fresnel’s sine and cosine integrals S(z) and C(z) for real (z = x) and complex (z = x + iy) arguments. The algorithms are based on series expansion for small values of |z|, expansion in Chebyshev subinterval polynomial approximations for intermediate real values (y = 0), toggled with sum of half-integer order Bessel functions for complex arguments, together with asymptotic series expansion for large values of |z|. The present algorithms, implemented in a Fortran elemental module, can be run using any of the single, double, or quadruple precision arithmetic. Results from the present code have been benchmarked versus comprehensive data tables, generated using Maple, Mathematica, and Matlab symbolic toolbox. Compared to a maximum of 16 significant digits in the literature, the present algorithms can calculate S(x) and C(x) up to 28 significant digits in the range xϵ[0, 106] and S(z) and C(z) to the same accuracy in the domain where |S(z)| and |C(z)| each is less than the largest finite floating point number in the precision under consideration with ∣argz∣<π2\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\mid \\mathit{\\arg}\\ z\\mid <\\frac{\\pi }{2}$$\\end{document} or more specifically 3×10-4<yx<2×103\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$3\ imes {10}^{-4}<\\left|\\frac{y}{x}\\right|<2\ imes {10}^3$$\\end{document} for S(z) and 4×10-4<yx<2.2×103\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$4\ imes {10}^{-4}<\\left|\\frac{y}{x}\\right|<2.2\ imes {10}^3$$\\end{document} for C(z).

Efficient multiple-precision computation of the scaled complementary error function and the Dawson integral

AbstractWe present algorithms to approximate the scaled complementary error function, $$\mathit{exp}\left({x}^{2}\right)erfc(x)$$ exp x 2 e r f c ( x ) , and the Dawson integral, $${e}^{{-x}^{2}}\underset{0}{\overset{x}{\int }}{e}^{{t}^{2}}dt$$ e - x 2 ∫ x 0 e t 2 d t , to the best accuracy in the standard single, double, and quadruple precision arithmetic. The algorithms are based on expansion in Chebyshev subinterval polynomial approximations together with expansion in terms of Taylor series and/or Laplace continued fraction. The present algorithms, implemented as Fortran elemental modules, have been benchmarked versus competitive algorithms available in the literature and versus functions built-in in modern Fortran compilers, in addition to comprehensive tables generated with variable precision computations using the Matlab™ symbolic toolbox. The present algorithm for calculating the scaled complementary error function showed an overall significant efficiency improvement (factors between 1.3 and 20 depending on the compiler and tested dataset) compared to the built-in function “Erfc_Scaled” in modern Fortran compilers, whereas the algorithm for calculating the Dawson integral is exceptional in calculating the function to 32 significant digits (compared to 19 significant digits reported in the literature) while being more efficient than competitive algorithms as well.

Mixed precision support in HPC applications: What about reliability?

In their quest for exascale and beyond, High-Performance Computing (HPC) systems continue becoming ever larger and more complex. Application developers, on the other hand, leverage novel methods to improve the efficiency of their own codes: a recent trend is the use of floating-point mixed precision, or the careful interlocking of single- and double-precision arithmetic, as a tool to improve performance as well as reduce network and memory boundedness. However, while it is known that modern HPC systems suffer hardware faults at daily rates, the impact of reduced precision on application reliability is yet to be explored. In this work we aim to fill this gap: first, we propose a qualitative survey to identify the branches of HPC where mixed precision is most popular. Second, we show the results of instruction-level fault injection experiments on a variety of representative HPC workloads, comparing vulnerability to Silent Data Errors (SDEs) under different numerical configurations. Our experiments indicate that use of single and mixed precision leads to comparatively more frequent and more severe SDEs, with concerning implications regarding their use on extreme-scale, fault-prone HPC platforms.

OrbitN: A Symplectic Integrator for Planetary Systems Dominated by a Central Mass—Insight into Long-term Solar System Chaos

Reliable studies of the long-term dynamics of planetary systems require numerical integrators that are accurate and fast. The challenge is often formidable because the chaotic nature of many systems requires relative numerical error bounds at or close to machine precision (∼10−16, double-precision arithmetic); otherwise, numerical chaos may dominate over physical chaos. Currently, the speed/accuracy demands are usually only met by symplectic integrators. For example, the most up-to-date long-term astronomical solutions for the solar system in the past (widely used in, e.g., astrochronology and high-precision geological dating) have been obtained using symplectic integrators. However, the source codes of these integrators are unavailable. Here I present the symplectic integrator orbitN (lean version 1.0) with the primary goal of generating accurate and reproducible long-term orbital solutions for near-Keplerian planetary systems (here the solar system) with a dominant mass M 0. Among other features, orbitN-1.0 includes M 0’s quadrupole moment, a lunar contribution, and post-Newtonian corrections (1PN) due to M 0 (fast symplectic implementation). To reduce numerical round-off errors, Kahan compensated summation was implemented. I use orbitN to provide insight into the effect of various processes on the long-term chaos in the solar system. Notably, 1PN corrections have the opposite effect on chaoticity/stability on a 100 Myr versus Gyr timescale. For the current application, orbitN is about as fast as or faster (factor 1.15–2.6) than comparable integrators, depending on hardware. 1 1 The orbitN source code (C) is available at http://github.com/rezeebe/orbitN.

Dynamic Precision Approach for Accelerating Large-Scale Eigenvalue Solvers in Electronic Structure Calculations on Graphics Processing Units.

Single precision (SP) arithmetic can be greatly accelerated as compared to double precision (DP) arithmetic on graphics processing units (GPUs). However, the use of SP in the whole process of electronic structure calculations is inappropriate for the required accuracy. We propose a 3-fold dynamic precision approach for accelerated calculations but still with the accuracy of DP. Here, SP, DP, and mixed precision are dynamically switched during an iterative diagonalization process. We applied this approach to the locally optimal block preconditioned conjugate gradient method to accelerate a large-scale eigenvalue solver for the Kohn-Sham equation. We determined a proper threshold for switching each precision scheme by examining the convergence pattern on the eigenvalue solver only with the kinetic energy operator of the Kohn-Sham Hamiltonian. As a result, we achieved up to 8.53× and 6.60× speedups for band structure and self-consistent field calculations, respectively, for test systems under various boundary conditions on NVIDIA GPUs.

Runge–Kutta–Nyström Pairs of Orders 8(6) for Use in Quadruple Precision Computations

The second-order system of non-stiff Initial Value Problems (IVP) is considered and, in particular, the case where the first derivatives are absent. This kind of problem is interesting since since it arises in many significant problems, e.g., in Celestial mechanics. Runge–Kutta–Nyström (RKN) pairs are perhaps the most successful approaches for solving such type of IVPs. To achieve a pair attaining orders eight and six, we have to solve a well-defined set of equations with respect to the coefficients. Here, we provide a simplified form of these equations in a robust algorithm. When creating such pairings for use in double precision arithmetic, numerous conditions are often satisfied. First and foremost, we strive to keep the coefficients’ magnitudes small to prevent accuracy loss. We may, however, allow greater coefficients when working with quadruple precision. Then, we may build pairs of orders eight and six with significantly smaller truncation errors. In this paper, a novel pair is generated that, as predicted, outperforms state-of-the-art pairs of the same orders in a collection of important problems.

Quasi-Fermi-Based Charge Transport Scheme for Device Simulation in Cryogenic, Wide Bandgap, and High-Voltage Applications

We present a novel approach to solving the transport problem in semiconductors. We reformulate the drift-diffusion (DD) equations in terms of the quasi-Fermi-energies as solution variables; a drastic increase in numerical stability is achieved, which permits the simulation of devices at cryogenic temperatures as well as wide bandgap devices using double precision arithmetic, instead of extended precision arithmetic which would otherwise be required to solve these applications using regular DD.

Superancillary Equations for Nonpolar Pure Fluids Modeled with the PC-SAFT Equation of State

Superancillary equations are presented for the PC-SAFT equation of state of Gross and Sadowski for nonpolar pure fluids. The equations cover the range of a number of segments from 1 to 64. These equations, formed as nested sets of Chebyshev expansions, represent the densities for vapor–liquid equilibria better than can be achieved by iterative calculations in double precision arithmetic. Furthermore, the expressions are more than 16 times faster to evaluate than the iterative calculations from the fastest thermodynamic property library. By their construction, the superancillary equations are guaranteed to be reliable over the entire temperature range. For user-friendliness, the functions have been packaged into a Python package available on PYPI.

Uplink Signal Detection via Look-Up Table-Based AMP for Massive MIMO Systems

We propose a novel look-up table (LUT)-based multi-user detection (MUD) algorithm to imitate approximate message passing (AMP) for uplink signal detection in massive multi-user multi-input multi-output (MU-MIMO) systems. When AMP is implemented with double-precision arithmetic, the dramatic increase in power consumption and process latency becomes a major barrier to practical application, as the system scale expands. To circumvent this issue, we design a novel referential AMP detector composed of many small LUTs cascaded hierarchically, where only informative integer-valued messages are exchanged on a factor graph (FG). Our method tracks the discrete distribution of the LUT outputs according to the multi-layer structure, and successively determines non-uniform quantization thresholds of the LUTs, using an unsupervised learning approach. Specifically, the thresholds at each level of the hierarchy, <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">i.e</i> ., layer, are optimized by clustering with the Lloyd-Max algorithm with initial values given by the k-means++ method, in order to minimize the performance degradation caused by quantization errors. The efficacy of the proposed referential AMP detector is confirmed by numerical simulations, which show that the proposed method significantly outperforms the state-of-the-art (SotA) in terms of the bit error rate (BER) performance in various channel models. The results also indicate the referential AMP is robust to changes in wireless channels, and then clarify the reason in terms of the algorithmic structure.

Singular value decomposition in extended double precision arithmetic

A well-known and successful algorithm to compute the singular value decomposition (SVD) of a matrix was published by Golub and Reinsch (Numer. Math. 14:403–420, 1970), together with an implementation in Algol. We give an updated implementation in extended double precision arithmetic in the C programming language. Extended double precision is native for Intel x86 processors and provides improved accuracy at full hardware speed. The complete program for computing the SVD is listed. Additionally, a comprehensive explanation of the original algorithm of Golub and Reinsch (Numer. Math. 14:403–420, 1970) is given at an elementary level without referring to the more general results of Francis (Comput. J. 4:265–271, 1961, 1962).

Stable factorization for phase factors of quantum signal processing

This paper proposes a new factorization algorithm for computing the phase factors of quantum signal processing. The proposed algorithm avoids root finding of high degree polynomials by using a key step of Prony&amp;apos;s method and is numerically stable in the double precision arithmetics. Experimental results are reported for Hamiltonian simulation, eigenstate filtering, matrix inversion, and Fermi-Dirac operator.

Runge–Kutta Embedded Methods of Orders 8(7) for Use in Quadruple Precision Computations

High algebraic order Runge–Kutta embedded methods are commonly used when stringent tolerances are demanded. Traditionally, various criteria are satisfied while constructing these methods for application in double precision arithmetic. Firstly we try to keep the magnitude of the coefficients low, otherwise we may experience loss of accuracy; however, when working in quadruple precision we may admit larger coefficients. Then we are able to construct embedded methods of orders eight and seven (i.e., pairs of methods) with even smaller truncation errors. A new derived pair, as expected, is performing better than state-of-the-art pairs in a set of relevant problems.

Compact formulae for three-center nuclear attraction integrals over exponential type functions

In any ab initio molecular orbital calculations, the major task involves the computation of the so-called molecular multi-center integrals. Multi-center integral calculations is a very challenging mathematical problem in nature. Quantum mechanics only determines which integrals we evaluate, but the techniques employed for their evaluations are entirely mathematical. The three-center nuclear attraction integrals occur in a very large number even for small molecules and are among of the most difficult molecular integrals to compute efficiently. In the present contribution, we report analytical expressions for the three-center nuclear attraction integrals over exponential type functions. We describe how to compute the formula to obtain an efficient evaluation in double precision arithmetic. This requires the rational minimax approximants that minimize the maximum error on the interval of evaluation.

On the Highly Accurate Evaluation of the Voigt/Complex Error Function with Small Imaginary Argument

A rapidly convergent series, based on Taylor expansion of the imaginary part of the complex error function, is presented for highly accurate approximation of the Voigt/complex error function with small imaginary argument y ≤ 0.1. Error analysis and run-time tests in double-precision arithmetic reveals that in the real and imaginary parts, the proposed algorithm provides an average accuracy exceeding 10−15 and 10−16, respectively, and the calculation speed is as fast as that reported in recent publications. An optimized MATLAB code providing rapid computation with high accuracy is presented.

4倍精度演算による複素数階微分の2階微分に関する基礎的検討

Methods for improving the accuracy of second order derivative calculations using the complex first-order derivative approximation (CSDA method) is proposed. In this paper, the following two approaches are investigated: (1) perturbation direction on the complex number plane, (2) use of double and quadruple precision floating point arithmetics. The above approaches were applied to several numerical examples to discuss the accuracy and computational costs. As the results, the perturbation of CSDA on θ = 45° had the best accuracy and stability throughout all the examples. Moreover, it was also found that when combined with perturbation on θ = 45° and quadruple precision floating point arithmetics, the second-order derivatives could be computed with more than accuracy of the double precision arithmetic in the examples.

Evolution of the Graphics Processing Unit (GPU)

Graphics processing units (GPUs) power today’s fastest supercomputers, are the dominant platform for deep learning, and provide the intelligence for devices ranging from self-driving cars to robots and smart cameras. They also generate compelling photorealistic images at real-time frame rates. GPUs have evolved by adding features to support new use cases. NVIDIA’s GeForce 256, the first GPU, was a dedicated processor for real-time graphics, an application that demands large amounts of floating-point arithmetic for vertex and fragment shading computations and high memory bandwidth. As real-time graphics advanced, GPUs became programmable. The combination of programmability and floating-point performance made GPUs attractive for running scientific applications. Scientists found ways to use early programmable GPUs by casting their calculations as vertex and fragment shaders. GPUs evolved to meet the needs of scientific users by adding hardware for simpler programming, double-precision floating-point arithmetic, and resilience.