Improving the Accuracy of the Fast Inverse Square Root by Modifying Newton–Raphson Corrections

Abstract

Direct computation of functions using low-complexity algorithms can be applied both for hardware constraints and in systems where storage capacity is a challenge for processing a large volume of data. We present improved algorithms for fast calculation of the inverse square root function for single-precision and double-precision floating-point numbers. Higher precision is also discussed. Our approach consists in minimizing maximal errors by finding optimal magic constants and modifying the Newton–Raphson coefficients. The obtained algorithms are much more accurate than the original fast inverse square root algorithm and have similar very low computational costs.

1. Introduction

Efficient performance of algebraic operations in the framework of floating-point arithmetic is a subject of considerable importance [1,2,3,4,5,6]. Approximations of elementary functions are crucial in scientific computing, computer graphics, signal processing, and other fields of engineering and science [7,8,9,10]. Our aim is to compute elementary functions at a very low computational cost without using memory resources. Direct evaluation of functions could be of interest in any systems where storage capabilities challenge the processing of a large volume of data. This problem is crucial, for instance, in high-energy physics experiments [11,12,13].

In this paper, we consider approximation and fast computation of the inverse square root function, which has numerous applications (see [8,10,14,15,16,17]), especially in 3D computer graphics, where it is needed for normalization of vectors [4,18,19]. The proposed algorithms are aimed primarily at floating-point platforms with limited hardware resources, such as microcontrollers, some field-programmable gate arrays (FPGAs), and graphics processing units (GPUs) that cannot use fast look-up table (LUT)-based hardware instructions, such as SSE (i.e., Streaming SIMD (single instruction, multiple data) Extensions) or Advanced Vector Extensions (AVX). We mean here devices and chips containing floating-point multipliers, adders–subtractors, and fused-multiply adders. Therefore, our algorithms can easily be implemented on such a platform. We also offer them as an alternative to library functions that provide full precision, but are very time consuming. This was the motivation for considering the cases of higher precision in Section 3.2. By selecting the precision and number of iterations, the desired accuracy can be obtained. We propose the use of our codes as direct insertions into more general algorithms without referring to the corresponding library of mathematical functions. In the double-precision mode, most modern processors do not have SSE instructions like rsqrt (such instructions appeared only in AVX-512, which is supported only by the latest processor models). In such cases, one can use our algorithms (with the appropriate number of iterations) as a fast alternative to the library function 1/sqrt(x).

In most cases, the initial seed needed to start the approximation is taken from a memory-consuming look-up table (LUT), although the so-called “bipartite table methods” (actually used on many current processors) make it possible to considerably lower the table sizes [20,21]. The “fast inverse square root” code works in a different way. It produces the initial seed in a cheap and effective way using the so-called magic constant [4,19,22,23,24,25]. We point out that this algorithm is still useful in numerous software applications and hardware implementations (see, e.g., [17,26,27,28,29,30]). Recently, we presented a new approach to the fast inverse square root code InvSqrt, presenting a rigorous derivation of the well-known code [31]. Then, this approach was used to construct a more accurate modification (called InvSqrt1) of the fast inverse square root (see [32]). It will be developed and generalized in the next sections, where we will show how to increase the accuracy of the InvSqrt code without losing its advantages, including the low computational cost. We will construct and test two new algorithms, InvSqrt2 and InvSqrt3.

The main idea of the algorithm InvSqrt consists in interpreting bits of the input floating-point number as an integer [31]. In this paper, we consider positive floating-point normal numbers

$$x=(1+m_x)2^{e_x},m_x\in [0,1),e_x\in \mathbb{Z},$$

(1)

and, in Section 3.1, we also consider subnormal numbers. We use the standard IEEE-754, where single-precision floating-point numbers are encoded with 32 bits. For positive numbers, the first bit is zero. The next eight bits encode $e_x$, and the remaining 23 bits represent the mantissa $m_x$. The same 32 bits can be treated as an integer $I_x$:

$$I_x=N_m(B+e_x+m_x),$$

(2)

where $N_m=2^{23}$ and $B=127$. In this case $B+e_x$ is a natural number not exceeding 254. The case of higher precision is analogous (see Section 3.2).

The crucial step of the algorithm InvSqrt consists in shifting all bits to the right by one bit and subtracting the result of this operation from a “magic constant” R (and the optimum value of R has to be guessed or determined). In other words,

$$I_{y_0}=R-\lfloor I_x/2\rfloor .$$

(3)

Originally, R was proposed as $0x5F3759DF$ (see [19,23]). Interpreted in terms of floating-point numbers, $I_{y_0}$ approximates the inverse square root function surprisingly well ($y_0\approx y=1/\sqrt{x}$). This trick works because (3) is close to dividing the floating-point exponent by $-2$. The number R is needed because the floating-point exponents are biased (see (2)).

The magic constant R is usually given as a hexadecimal integer. The same bits encode the floating-point number $R_f$ with an exponent $e_R$ and mantissa $m_R$. According to (1), $R_f=(1+m_R)2^{e_R}$. In [31], we have shown that if $e_R=\frac{1}{2}\left(B-1\right)$ (e.g., $e_R=63$ in the 32-bit case), then the function (3) (defined on integers) is equivalent to the following piece-wise linear function (when interpreted in terms of corresponding floating-point numbers):

$${\tilde{y}}_0(\tilde{x},t)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{8}\left(6-2\tilde{x}+t\right)}\hfill & \mathrm{for}\tilde{x}\in [1,2)\hfill \\ {\displaystyle \frac{1}{8}\left(4+t-\tilde{x}\right)}\hfill & \mathrm{for}\tilde{x}\in [2,t)\hfill \\ {\displaystyle \frac{1}{16}\left(8+t-\tilde{x}\right)}\hfill & \mathrm{for}\tilde{x}\in [t,4)\hfill \end{array}\right.$$

(4)

where

$$t=2+4m_R+2{\mu }_{\tilde{x}}{N}_m^{-1},$$

(5)

$m_R$ is the mantissa of R (i.e., $m_R:=N_m^{-1}R-\lfloor N_m^{-1}R\rfloor$), and, finally, ${\mu }_{\tilde{x}}=0$ or ${\mu }_{\tilde{x}}=1$ depending on the parity of the last digit of the mantissa of $\tilde{x}$.

The function ${\mu }_{\tilde{x}}$ is two-valued, so a given parameter t may correspond to either two values of R or one value of R (when the term containing ${\mu }_{\tilde{x}}$ has no influence on the bits of the mantissa $m_R$). The function $y=1/\sqrt{x}$, the function (3), and all Newton–Raphson corrections considered below are invariant under the scaling $\tilde{x}=2^{-2n}x$ and $\tilde{y}=2^{n}y$ for any integer n. Therefore, we can confine ourselves to numbers from the interval $[1,4)$. Here and in the following, the tilde always denotes quantities defined on the interval $[1,4)$.

In this paper, we focus on the Newton–Raphson corrections, which form the second part of the InvSqrt code. Following and developing ideas presented in our recent papers [31,32], we propose modifications of the Newton–Raphson formulas, which result in algorithms that have the same or similar computational cost as/to InvSqrt, but improve the accuracy of the original code, even by several times. The modifications consist in changing both the Newton–Raphson coefficients and the magic constant. Moreover, we extend our approach to subnormal numbers and to higher-precision cases.

2. Modified Newton–Raphson Formulas

The standard Newton–Raphson corrections ${\tilde{y}}_1$ and ${\tilde{y}}_2$ for the zeroth approximation ${\tilde{y}}_0$ given by (4) are given by the following formulas:

$$\begin{array}{c}{\displaystyle {\tilde{y}}_1(\tilde{x},t)=\frac{3}{2}{\tilde{y}}_0(\tilde{x},t)-\frac{1}{2}\tilde{x}{\tilde{y}}_0^3(\tilde{x},t),}\\ {\displaystyle {\tilde{y}}_2(\tilde{x},t)=\frac{3}{2}{\tilde{y}}_1(\tilde{x},t)-\frac{1}{2}\tilde{x}{\tilde{y}}_1^3(\tilde{x},t),}\end{array}$$

(6)

(analogous formulas hold for the next corrections as well; see [31]). The relative error functions ${\tilde{\delta }}_j(\tilde{x},t)$ (where $j=0,1,2,\dots$) can be expressed as:

$${\tilde{\delta }}_j(\tilde{x},t)=\sqrt{\tilde{x}}{\tilde{y}}_j(\tilde{x},t)-1.$$

(7)

The function ${\tilde{\delta }}_0(\tilde{x},t)$, which is very important for the further analysis, is thoroughly described and discussed in [31]. Using (7), we substitute ${\tilde{y}}_j=(1+{\tilde{\delta }}_j)/\sqrt{\tilde{x}}$ (for $j=0,1,2,\dots$) into (6), $\tilde{x}$ cancels out, and the formulas (6) assume the following form:

$${\tilde{\delta }}_j=-\frac{1}{2}{\tilde{\delta }}_{j-1}^2(3+{\tilde{\delta }}_{j-1}),(j=1,2,\dots ),$$

(8)

where ${\tilde{\delta }}_j={\tilde{\delta }}_j(\tilde{x},t)$. We immediately see that every correction increases the accuracy, even by several orders of magnitude (due to the factor ${\tilde{\delta }}_{j-1}^2$). Thus, a very small number of corrections is sufficient to reach the machine precision (see the end of Section 4).

The above approximations depend on the parameter t (which can be expressed by the magic constant R, see (5)). The best approximation is obtained for $t=t_k$ minimizing $\left|\right|{\tilde{\delta }}_k\left(t\right)\left|\right|$, i.e.,

$$\left|\right|{\tilde{\delta }}_k\left(t_k\right)\left|\right|=\underset{t\in (2,4)}{\inf }\left|\right|{\tilde{\delta }}_k\left(t\right)\left|\right|\equiv \underset{t\in (2,4)}{\inf }\left(\underset{\tilde{x}\in [1,4)}{\sup }\left|{\tilde{\delta }}_k(\tilde{x},t)\right|\right).$$

(9)

In this paper we confine ourselves to the case $t=t_1$ (i.e., we assume $t_2=t_1$) because the more general case (where the magic constant is also optimized with respect to the assumed number of iterations) is much more cumbersome, and the related increase in accuracy is negligible. Then, we get

$$t_1^{\left(0\right)}\approx 3.7298003,R^{\left(0\right)}=0x5F375A86,$$

(10)

for details, see [31]. The theoretical relative errors are given by

$${\Delta }_{1\max }^{\left(0\right)}\equiv \left|\right|{\tilde{\delta }}_1\left(t_1\right)|\approx 1.75118·{10}^{-3},{\Delta }_{2\max }^{\left(0\right)}\equiv \left|\right|{\tilde{\delta }}_2\left(t_1\right)|\approx 4.60·{10}^{-6}.$$

(11)

The superscript $\left(0\right)$ indicates values corresponding to the algorithm InvSqrt (other superscripts will denote modifications of this algorithm).

The idea of increasing the accuracy by a modification of the Newton–Raphson formulas is motivated by the fact that ${\tilde{\delta }}_k(\tilde{x},t)⩽0$ for any $\tilde{x}$ (see [31,32]). Therefore, we can try to shift the graph of ${\tilde{\delta }}_1$ upwards (making it more symmetric with respect to the horizontal axis). Then, the errors of the first correction are expected to decrease twice and the errors of the second correction are expected to decrease by about eight times (for more details, see [32]). Indeed, according to (8), reducing the first correction by a factor of 2 will reduce the second correction by a factor of 4. The second correction is also non-positive, so we may shift the graph of ${\tilde{\delta }}_2$, once more improving the accuracy by the factor of 2. This procedure can be formalized by postulating the following modification of the Newton–Raphson formulas (6):

$$\begin{array}{c}{\tilde{y}}_1=\frac{1}{2}{\tilde{y}}_0(3-{\tilde{y}}_0^2\tilde{x})+\frac{1}{2}{d}_1\left(a_1{\tilde{y}}_0+b_1{\tilde{y}}_1\right),\hfill \\ {\tilde{y}}_2=\frac{1}{2}{\tilde{y}}_1(3-{\tilde{y}}_1^2\tilde{x})+\frac{1}{2}{d}_2\left(a_2{\tilde{y}}_1+b_2{\tilde{y}}_2\right),\hfill \end{array}$$

(12)

where $a_k+b_k=1$ for $k=1,2$. Thus, we have four independent parameters ($d_1$, $d_2$, $a_1$, and $a_2$) to be determined. In other words,

$$\begin{array}{c}{\displaystyle {\tilde{y}}_1=c_{11}{\tilde{y}}_0-c_{21}\tilde{x}{\tilde{y}}_0^3,}\hfill \\ {\displaystyle {\tilde{y}}_2=c_{12}{\tilde{y}}_1-c_{22}\tilde{x}{\tilde{y}}_1^3,}\hfill \end{array}$$

(13)

where four coefficients $c_{jk}$ can be expressed by the four coefficients $a_k$ and $d_k$:

$$c_{1k}=\frac{3+a_k{d}_k}{2-(1-a_k)d_k},c_{2k}=\frac{1}{2-(1-a_k)d_k}(k=1,2).$$

(14)

We point out that the Newton–Raphson corrections and any of their modifications of the form (13) are obviously invariant with respect to the scaling mentioned at the end of Section 1. Therefore, we can continue to confine our analysis to the interval $[1,4)$.

Below, we present three different algorithms (InvSqrt1, InvSqrt2, InvSqrt3) constructed along the above principles (the last two of them are first introduced in this paper). They will be denoted by superscripts in parentheses, e.g., ${\tilde{y}}_k^{\left(N\right)}$ means the kth modified Newton–Raphson correction to the algorithm InvSqrt N. We always assume that the zeroth approximation is given by (4), i.e.,

$${\tilde{y}}_0^{\left(N\right)}={\tilde{y}}_0(N=1,2,3),$$

(15)

and relative error functions, ${\Delta }_j^{\left(N\right)}$, are expressed as

$${\Delta }_j^{\left(N\right)}(\tilde{x},t)=\sqrt{\tilde{x}}{\tilde{y}}_j^{\left(N\right)}(\tilde{x},t)-1.$$

(16)

We point out that the coefficients of our algorithms are obtained without taking rounding errors into account. This issue will be shortly discussed at the end of Section 4.

2.1. Algorithm InvSqrt1

Assuming $a_1=a_2=0$ and $b_1=b_2=1$, we transform (12) into

$$\begin{array}{c}{\displaystyle {\tilde{y}}_1^{\left(1\right)}=\frac{1}{2}{\tilde{y}}_0^{\left(1\right)}\left(3-{\left({\tilde{y}}_0^{\left(1\right)}\right)}^2\tilde{x}\right)+\frac{1}{2}{d}_1{\tilde{y}}_1^{\left(1\right)},}\hfill \\ {\displaystyle {\tilde{y}}_2^{\left(1\right)}=\frac{1}{2}{\tilde{y}}_1^{\left(1\right)}\left(3-{\left({\tilde{y}}_1^{\left(1\right)}\right)}^2\tilde{x}\right)+\frac{1}{2}{d}_2{\tilde{y}}_2^{\left(1\right)}.}\hfill \end{array}$$

(17)

Therefore, ${\tilde{y}}_1^{\left(1\right)}$ and ${\tilde{y}}_2^{\left(1\right)}$ depend on $\tilde{x}$, t, $d_1$, and $d_2$. Parameters $t=t_1^{\left(1\right)}$ and $d_1=d_1^{\left(1\right)}$ are determined by minimization of $\left|\right|{\Delta }_1^{\left(1\right)}(\tilde{x},t)\left|\right|$. Then, the parameter $d_2=d_2^{\left(1\right)}$ is determined by minimization of $\left|\right|{\Delta }_2^{\left(1\right)}(\tilde{x},t_1^{\left(1\right)})\left|\right|$ (for details, see [32]). As a result, we get:

$$d_1^{\left(1\right)}\approx 1.75118·{10}^{-3},d_2^{\left(1\right)}\approx 1.15234·{10}^{-6},$$

(18)

and $t_1^{\left(1\right)}=t_1^{\left(0\right)}$ (see (10)). Therefore, $R^{\left(1\right)}=R^{\left(0\right)}$, i.e., InvSqrt1 has the same magic constant as InvSqrt. The theoretical relative errors are given by

$$\begin{array}{c}{\Delta }_{1\max }^{\left(1\right)}\equiv \left|\right|{\Delta }_1^{\left(1\right)}(\tilde{x},t_1^{\left(1\right)})\left|\right|\approx 0.87636·{10}^{-3},\hfill \\ {\Delta }_{2\max }^{\left(1\right)}\equiv \left|\right|{\Delta }_2^{\left(1\right)}(\tilde{x},t_1^{\left(1\right)})\left|\right|\approx 5.76·{10}^{-7}.\hfill \end{array}$$

(19)

The algorithm (17) can be written in the form (13), where:

$$c_{1k}^{\left(1\right)}=\frac{3}{2-d_k^{\left(1\right)}},c_{2k}^{\left(1\right)}=\frac{1}{2-d_k^{\left(1\right)}}(k=1,2).$$

(20)

Taking into account numerical values for $d_1^{\left(1\right)}$ and $d_2^{\left(1\right)}$, we obtain the following values of the parameters $c_{jk}^{\left(1\right)}$:

$$\begin{array}{c}{c}_{11}^{\left(1\right)}\approx 1.5013145387528147176730252470373223,\hfill \\ c_{21}^{\left(1\right)}\approx 0.50043817958427157255767508234577407,\hfill \\ c_{12}^{\left(1\right)}\approx 1.5000008642589575005473878767725752,\hfill \\ c_{22}^{\left(1\right)}\approx c_{21}^{\left(1\right)}·0.99912498383253616899527502360939620.\hfill \end{array}$$

(21)

This large number of digits, which is much higher than that needed for the single-precision computations, will be useful later in the case of higher precision.

Thus, finally, we obtained a new algorithm InvSqrt1 that has the same structure as InvSqrt, but with different values of numerical coefficients (see [32]). In the case of two iterations, the code $InvSqrt1$ has more algebraic operations (one additional multiplication) in comparison to InvSqrt.

2.2. InvSqrt2 Algorithm

Assuming $a_1=a_2=1$ and $b_1=b_2=0$, we transform (12) into

$$\begin{array}{c}{\displaystyle {\tilde{y}}_1^{\left(2\right)}=\frac{1}{2}{\tilde{y}}_0^{\left(2\right)}\left(3-{\left({\tilde{y}}_0^{\left(2\right)}\right)}^2\tilde{x}\right)+\frac{1}{2}{d}_1{\tilde{y}}_0^{\left(2\right)},}\hfill \\ {\displaystyle {\tilde{y}}_2^{\left(2\right)}=\frac{1}{2}{\tilde{y}}_1^{\left(2\right)}\left(3-{\left({\tilde{y}}_1^{\left(2\right)}\right)}^2\tilde{x}\right)+\frac{1}{2}{d}_2{\tilde{y}}_1^{\left(2\right)},}\hfill \end{array}$$

(22)

where ${\tilde{y}}_1^{\left(2\right)}$ and ${\tilde{y}}_2^{\left(2\right)}$ depend on $\tilde{x}$, t, $d_1$, and $d_2$.

Parameters $t=t_1^{\left(2\right)}$ and $d_1=d_1^{\left(2\right)}$ are determined by the minimization of $\left|\right|{\Delta }_1^{\left(2\right)}(\tilde{x},t)\left|\right|$. Then, the parameter $d_2=d_2^{\left(2\right)}$ is determined by the minimization of $\left|\right|{\Delta }_2^{\left(2\right)}(\tilde{x},t_1^{\left(2\right)})\left|\right|$ (see Appendix A.1 for details). As a result, we get:

$$d_1^{\left(2\right)}=1.75791023259·{10}^{-3},d_2^{\left(2\right)}\simeq 1.159352515·{10}^{-6},$$

(23)

and

$$t_1^{\left(2\right)}\equiv t^{\left(2\right)}\approx 3.73157124016,R^{\left(2\right)}=1597466888=0x5F376908.$$

(24)

The theoretical relative errors are given by

$$\begin{array}{c}{\Delta }_{1\max }^{\left(2\right)}\equiv \left|\right|{\Delta }_1^{\left(2\right)}(\tilde{x},t_1^{\left(2\right)})\left|\right|\approx 0.87908·{10}^{-3},\hfill \\ {\Delta }_{2\max }^{\left(2\right)}\equiv \left|\right|{\Delta }_2^{\left(2\right)}(\tilde{x},t_1^{\left(2\right)})\left|\right|\approx 5.80·{10}^{-7}.\hfill \end{array}$$

(25)

The coefficients in (13) are given by

$$c_{1k}^{\left(2\right)}=\frac{3+d_k^{\left(2\right)}}{2},c_{2k}^{\left(2\right)}=\frac{1}{2}.$$

(26)

Taking into account the numerical values for $d_1^{\left(2\right)}$ and $d_2^{\left(2\right)}$ (see (23)), we obtain the following values of the parameters $c_{jk}^{\left(2\right)}$:

$$\begin{array}{c}{c}_{11}^{\left(2\right)}\approx 1.5008789551163345746409291568502392,\hfill \\ c_{12}^{\left(2\right)}\approx 1.5000005796762576644996810350809289,\hfill \\ c_{21}^{\left(2\right)}=c_{22}^{\left(2\right)}=0.5,\hfill \end{array}$$

(27)

where the large number of digits will be useful later in the case of higher precision. Thus, we completed the derivation of the code $InvSqrt2$:

1. floatInvSqrt2(float x){

2.   float halfx = 0.5f*x;

3.   int i = *(int*) &x;

4.   i = 0x5F376908 - (i$>>$1);

5.   float y = *(float*) &i;

6.   y* = 1.50087896f - halfx*y*y;

7.   y* = 1.50000057f - halfx*y*y;

8.   return y;

9. }

The code InvSqrt2 contains a new magic constant ($R^{\left(2\right)}$) and has two lines (6 and 7) that were modified in comparison with the code $InvSqrt$. We point out that $InvSqrt2$ has the same number of algebraic operations as $InvSqrt$.

2.3. InvSqrt3 Algorithm

Now, we consider the algorithm (13) in its most general form:

$$\begin{array}{c}{\tilde{y}}_1^{\left(3\right)}=k_1{\tilde{y}}_0^{\left(3\right)}\left(k_2-\tilde{x}{\left({\tilde{y}}_0^{\left(3\right)}\right)}^2\right),\hfill \\ {\tilde{y}}_2^{\left(3\right)}=k_3{\tilde{y}}_1^{\left(3\right)}\left(k_4-\tilde{x}{\left({\tilde{y}}_1^{\left(3\right)}\right)}^2\right),\hfill \end{array}$$

(28)

where $k_j$, $k_2$, $k_3$, and $k_4$ are constant. In Appendix A.2, we determine parameters $t=t_1^{\left(3\right)}$, $k_1$, and $k_2$ by minimization of $\left|\right|{\Delta }_1^{\left(3\right)}(\tilde{x},t)\left|\right|$. Then, the parameters $k_3$ and $k_4$ are determined by minimization of $\left|\right|{\Delta }_2^{\left(3\right)}(\tilde{x},t_1^{\left(3\right)})\left|\right|$. As a result, we get:

$$\begin{array}{c}{k}_1\approx 0.70395201,k_2\approx 2.3892451,\hfill \\ k_3\approx 0.50000005,k_4\approx 3.0000004,\hfill \end{array}$$

(29)

and

$$t_1^{\left(3\right)}\equiv t^{\left(3\right)}=3,R^{\left(3\right)}=1595932672=0x5F200000.$$

(30)

The theoretical relative errors are given by

$$\begin{array}{c}{\Delta }_{1\max }^{\left(3\right)}\equiv \left|\right|{\Delta }_1^{\left(3\right)}(\tilde{x},t_1^{\left(3\right)})\left|\right|\approx 0.65007·{10}^{-3},\hfill \\ {\Delta }_{2\max }^{\left(3\right)}\equiv \left|\right|{\Delta }_2^{\left(3\right)}(\tilde{x},t_1^{\left(3\right)})\left|\right|\approx 3.17·{10}^{-7}.\hfill \end{array}$$

(31)

They are significantly smaller (by $26\%$ and $45\%$, respectively) than the analogous errors for InvSqrt1 and InvSqrt2 (see (19) and (25)). The comparison of error functions for InvSqrt2 and InvSqrt3 (in the case of one correction) is presented in Figure 1.

The numerical values of coefficients $c_{ij}^{\left(3\right)}$ (compare with (13)) are given by:

$$\begin{array}{c}{c}_{11}^{\left(3\right)}=k_1{k}_2\approx 1.68191390868723079,\hfill \\ c_{21}^{\left(3\right)}=k_1\approx 0.703952009104829370,\hfill \\ c_{12}^{\left(3\right)}=k_3{k}_4\approx 1.50000036976749938,\hfill \\ c_{22}^{\left(3\right)}=k_3\approx 0.500000052823927419.\hfill \end{array}$$

(32)

Thus, we obtained the following code, called InvSqrt3:

1. floatInvSqrt3(float x){

2.   int i = *(int*) &x;

3.   i = 0x5F200000 - (i$>>$1);

4.   float y = *(float*) &i;

5.   y* = 1.68191391f - 0.703952009f*x*y*y;

6.   y* = 1.50000036f - 0.500000053f*x*y*y;

7.   return y;

8. }

The code InvSqrt3 has the same number of multiplications as InvSqrt1, which means that it is slightly more expensive than InvSqrt and InvSqrt2.

3. Generalizations

The codes presented in Section 2 can only be applied to normal numbers (1) of the type float. In this section, we show how to extend these results to subnormal numbers and to higher-precision formats.

3.1. Subnormal Numbers

Subnormal numbers are smaller than any normal number of the form of (1). In the single-precision case, positive subnormals can be represented as $m_x·2^{-126}$, where $m_x\in (0,1)$. They can also be characterized by nine first bits equal to zero (which also includes the case where $x=0$). In order to identify subnormals, we will make a bitwise conjunction (AND) of a given number with the integer $0x7f800000$, which has all eight exponent bits equal to 1 and all 23 mantissa bits equal to 0. This bitwise conjunction is zero if and only if the given number is subnormal (including 0).

In the case of the single precision, the multiplication by $2^{24}$ transforms any subnormal number into a normal number. Therefore, we make this transformation; then, we apply one of our algorithms and, finally, make the inverse transformation (i.e., multiplying the result by $2^{-12}$). Thus, we get an approximated value of the inverse square root of the subnormal number. Note that $2^{24}$ is the smallest power of 2 with an even exponent that transforms all subnormals into normal numbers.

In the case of InvSqrt3, the procedure described above can be written in the form of the following code.

1. floatInvSqrt3s(float x){

2.   int i = *(int*) &x;

3.   int k = i & 0x7f800000;

4.   if (k==0) {

5.    x = 16777216.f*x; //16777216.f=pow(2.0f, 24)

6.    i = *(int*) &x;

7.   }

8.   i = 0x5F200000 - (i$>>$1);

9.   float y = *(float*) &i;

10.   y* = 1.68191391f - 0.703952009f*x*y*y;

11.   y* = 1.50000036f - 0.500000053f*x*y*y;

12.   if (k==0) return 4096.f*y; //4096.f=pow(2.0f, 12)

13.   return y;

14. }

The maximum relative errors for this code are presented in Section 4 (see

Table 1. Relative numerical errors for the first and second corrections in the case of the type float (compiler 32-bit) for subnormal numbers.

Algorithm	${\mathbf{\Delta }}_{1,\mathit{N}\min }^{\left(\mathit{i}\right)}$	${\mathbf{\Delta }}_{1,\mathit{N}\max }^{\left(\mathit{i}\right)}$	${\mathbf{\Delta }}_{2,\mathit{N}\min }^{\left(\mathit{i}\right)}$	${\mathbf{\Delta }}_{2,\mathit{N}\max }^{\left(\mathit{i}\right)}$
InvSqrt1	$-0.87642\times {10}^{-3}$	$0.87644\times {10}^{-3}$	$-0.66221\times {10}^{-6}$	$0.63442\times {10}^{-6}$
InvSqrt2	$-0.87916\times {10}^{-3}$	$0.87911\times {10}^{-3}$	$-0.62060\times {10}^{-6}$	$0.65285\times {10}^{-6}$
InvSqrt3	$-0.65016\times {10}^{-3}$	$0.65006\times {10}^{-3}$	$-0.38701\times {10}^{-6}$	$0.35196\times {10}^{-6}$

).

3.2. Higher Precision

The above analysis was confined to the single-precision floating-point format. This is sufficient for many applications (especially microcontrollers), although the double-precision format is more popular. A trade-off between accuracy, computational cost, and memory usage is welcome [33]. In this subsection, we extend our analysis to double- and higher-precision formats. The calculations are almost the same. We just have to compute all involved constants with an appropriate accuracy. Low-bit arithmetic cases could be treated in exactly the same way. In this paper, however, we are focused on increasing the accuracy and on possible applications in distributed systems, so only the cases of higher precision are explicitly presented.

We present detailed results for double precision and some results (magic constants) for quadruple precision. Performing computations in C, we use the GCC Quad-Precision Math Library (working with numbers of type _float128). The crucial point is to express the magic constant R through the corresponding parameter t, which can be done with the formula:

$$R=N_m(3B-1)/2+\lfloor N_m(t-2)/4-{\mu }_{\tilde{x}}/2\rceil$$

(33)

where ${\mu }_{\tilde{x}}\in \{0,1\}$, t depends on the considered algorithm and $N_m$ and B depend on the precision format used. Namely,

$$\begin{array}{ccc}\mathrm{Sin}\mathrm{gle}\mathrm{precision}\left(32\text{-bit}\right):\hfill & N_m=2^{23},\hfill & B=2^7-1,\hfill \\ \mathrm{Double}\mathrm{precision}\left(64\text{-bit}\right):\hfill & N_m=2^{52},\hfill & B=2^{10}-1,\hfill \\ \mathrm{Quadruple}\mathrm{precision}\left(128\text{-bit}\right):\hfill & N_m=2^{112},\hfill & B=2^{14}-1.\hfill \end{array}$$

(34)

In the case of the zeroth approximation (without Newton–Raphson corrections), the parameter t is given by:

$$t_0=3.7309795598377727818740863479840422,$$

(35)

which can be compared with [31]. The corresponding magic constants computed from the formula (33) read:

$$\begin{array}{c}32\mathrm{-bit}:R^{\left(0\right)}=0x5F37642F\hfill \\ 64\mathrm{-bit}:R^{\left(0D\right)}=0x5FE6EC85E7DE30DA\hfill \\ 128\mathrm{-bit}:R^{\left(0Q\right)}=0x5FFE6EC85E7DE30DAABC602711840B0F.\hfill \end{array}$$

(36)

In this paper, we focus on the case of Newton–Raphson corrections, where the value of the parameter t may depend on the algorithm. For InvSqrt and InvSqrt1, we have:

$$t_1^{\left(0\right)}=t_1^{\left(1\right)}=3.7298003391605705687151317499871860,$$

(37)

(see (Section 2.1); compare with [31,32]). Then, (33) yields the following magic constants:

$$\begin{array}{c}32\mathrm{-bit}:R^{\left(1\right)}=0x5F375A86\hfill \\ 64\mathrm{-bit}:R^{\left(1D\right)}=0x5FE6EB50C7B537A9\hfill \\ 128\mathrm{-bit}:R^{\left(1Q\right)}=0x5FFE6EB50C7B537A9CD9F02E504FCFC0.\hfill \end{array}$$

(38)

Actually, the above value of R in the 64-bit case (i.e., $R^{\left(1D\right)}$) corresponds to ${\mu }_{\tilde{x}}=1$ (the same value of R was obtained by Robertson for InvSqrt[24] with a different method). For ${\mu }_{\tilde{x}}=0$, we got an R greater by 1 (other results reported in this section do not depend on ${\mu }_{\tilde{x}}$). In the 128-bit case, Robertson obtained an R that was 1 less than our value (i.e., $R^{\left(1Q\right)}$).

In the case of InvSqrt2, we have

$$t^{\left(2\right)}=3.7315712401613957182292407381942955$$

(39)

(compare with (A11)), which yields:

$$\begin{array}{c}32\mathrm{-bit}:R^{\left(2\right)}=0x5F376908\hfill \\ 64\mathrm{-bit}:R^{\left(2D\right)}=0x5FE6ED2102DCBFDA\hfill \\ 128\mathrm{-bit}:R^{\left(2Q\right)}=0x5FFE6ED2102DCBFDA59415059AC483B5.\hfill \end{array}$$

(40)

Finally, for InvSqrt3, we obtained:

$$t^{\left(3\right)}=3$$

(41)

(see (A29)). The corresponding magic constants are given by:

$$\begin{array}{c}32\mathrm{-bit}:R^{\left(3\right)}=0x5F200000\hfill \\ 64\mathrm{-bit}:R^{\left(3D\right)}=0x5FE400000000000C\hfill \\ 128\mathrm{-bit}:R^{\left(3Q\right)}=0x5FFE4000000000000000000000000000.\hfill \end{array}$$

(42)

The parameters of the modified Newton–Raphson corrections for the higher-precision codes can be computed from the theoretical formulas used in the single-precision cases, taking into account an appropriate number of significant digits. In numerical experiments, we tested the algorithms InvSqrt1D, InvSqrt2D, and InvSqrt3D with the magic constants $R^{\left(1D\right)}$, $R^{\left(2D\right)}$, and $R^{\left(3D\right)}$, respectively, and the following coefficients in the modified Newton–Raphson iterations (compare with (21), (27), and (32), respectively):

$$\begin{array}{c}{c}_{11}^{\left(1D\right)}=1.50131453875281472,\hfill \\ c_{21}^{\left(1D\right)}=0.500438179584271573,\hfill \\ c_{12}^{\left(1D\right)}=1.50000086425895750,\hfill \\ c_{22}^{\left(1D\right)}=c_{21}^{\left(1D\right)}·0.999124983832536169.\hfill \end{array}$$

(43)

$$\begin{array}{c}{c}_{11}^{\left(2D\right)}=1.50087895511633457,\hfill \\ c_{12}^{\left(2D\right)}=1.50000057967625766,\hfill \\ c_{21}^{\left(2D\right)}=c_{22}^{\left(2D\right)}=0.5,\hfill \end{array}$$

(44)

$$\begin{array}{c}{c}_{11}^{\left(3D\right)}=1.68191390868723079,\hfill \\ c_{21}^{\left(3D\right)}=0.703952009104829370,\hfill \\ c_{12}^{\left(3D\right)}=1.50000036976749938,\hfill \\ c_{22}^{\left(3D\right)}=0.500000052823927419.\hfill \end{array}$$

(45)

The algorithm InvSqrt and its improved versions are usually implemented in the single-precision case with no more than two Newton–Raphson corrections. However, in the case of higher precision, higher accuracy of the result is welcome. Then, a higher number of modified Newton–Raphson iterations could be considered. As an example, we present the algorithm InvSqrt2D with four iterations:

1. doubleInvSqrt2D(double x){

2.   double halfx=0.5*x;

3.   long long i=*(long long*) &x;

4.   i=0x5FE6ED2102DCBFDA - (i$>>$1);

5.   double y =*(double*) &i;

6.   y* = 1.50087895511633457 - halfx*y*y;

7.   y* = 1.50000057967625766 - halfx*y*y;

8.   y* = 1.5000000000002520 - halfx*y*y;

9.   y* = 1.5000000000000000 - halfx*y*y;

10.   return y;

11. }

By removing Line 9, we obtain the code InvSqrt2D with three iterations, and by also removing Line 8, we get the code defined by (44). The maximum relative errors for this code are presented in Section 4 (see (52)).

4. Numerical Experiments

The numerical tests for the codes derived and presented in this paper were performed on an Intel Core i5-3470 processor using the TDM-GCC 4.9.2 32-bit compiler (when repeating these tests on the Intel i7-5700 processor, we obtained the same results, and comparisons with some other processors and compilers are given in Appendix B). In this section, we discuss round-off errors for the algorithms InvSqrt2 and InvSqrt3 (the case of single precision and two Newton–Raphson iterations) and then present the final results of analogous analysis for other codes described in this paper.

Applying algorithms InvSqrt2 and InvSqrt3, we obtain relative errors that differ slightly, due to round-off errors, from their analytical values (see Figure 2 and Figure 3; compare with [32] for an analogous discussion concerning InvSqrt1). Although we present only randomly chosen values in the figures, calculations were done for all float numbers x such that $e_x\in [-126,128)$.

The errors of numerical values returned by $InvSqrt2$

$${\Delta }_{2;N}^{\left(2\right)}\left(x\right)=\mathrm{sqrt}\left(x\right)\ast \mathit{InvSqrt}\mathit{2}\left(x\right)-1$$

(46)

belong (for $e_x\ne -126$) to the interval $(-6.21·{10}^{-7},6.53·{10}^{-7})$. For $e_x=-126$, we get a wider interval: $[-6.46·{10}^{-7},6.84·{10}^{-7}]$. These errors differ from the errors of ${\tilde{y}}_2^{\left(3\right)}(\tilde{x},t^{\left(2\right)})$, which were determined analytically (compare (25)). We define:

$${\epsilon }^{\left(2\right)}\left(\tilde{x}\right)=\frac{\mathit{InvSqrt}\mathit{2}\left(\mathit{x}\right)-{\tilde{y}}_2^{\left(2\right)}(\tilde{x},t^{\left(2\right)})}{{\tilde{y}}_2^{\left(2\right)}(\tilde{x},t^{\left(2\right)})}.$$

(47)

This function, representing the observed blur of the float approximation of the InvSqrt2 output, is symmetric with respect to its mean value

$$\langle {\epsilon }^{\left(2\right)}\rangle =2^{-1}{N}_m^{-1}\sum _{x\in [1,4)}{\epsilon }^{\left(2\right)}\left(\tilde{x}\right)=1.636·{10}^{-8}$$

(48)

(see the right part of Figure 2), and covers the following range of values:

$${\epsilon }^{\left(2\right)}\left(\tilde{x}\right)\in [-4.333·{10}^{-8},7.596·{10}^{-8}].$$

(49)

Analogous results for the code InvSqrt3 read:

$$\langle {\epsilon }^{\left(3\right)}\rangle =2^{-1}{N}_m^{-1}\sum _{x\in [1,4)}{\epsilon }^{\left(3\right)}\left(\tilde{x}\right)=-1.890·{10}^{-8}$$

(50)

$${\epsilon }^{\left(3\right)}\left(\tilde{x}\right)\in [-7.850·{10}^{-8},4.067·{10}^{-8}].$$

(51)

The results produced by the same hardware with a 64-bit compiler have a greater amplitude of the error oscillations as compared with the 32-bit case (also compare Appendix B).

The maximum errors for the code InvSqrt and all codes presented in the previous sections are given in

Table 2. Relative numerical errors for the first correction in the case of type float (compiler 32-bit). In the case of type double, the errors are equal to theoretical errors $\pm {\Delta }_{1\max }^{\left(i\right)}$ up to the accuracy given in the table.

Algorithm	i	${\mathbf{\Delta }}_{1\max }^{\left(\mathit{i}\right)}$	${\mathbf{\Delta }}_{1,\mathit{N}\min }^{\left(\mathit{i}\right)}$	${\mathbf{\Delta }}_{1,\mathit{N}\max }^{\left(\mathit{i}\right)}$	${\mathbf{\Delta }}_{1,\mathit{N}}^{\left(\mathit{i}\right)}$
InvSqrt	0	$1.75118\times {10}^{-3}$	$-1.75124\times {10}^{-3}$	$0.00008\times {10}^{-3}$	$1.75124\times {10}^{-3}$
InvSqrt1	1	$0.87636\times {10}^{-3}$	$-0.87642\times {10}^{-3}$	$0.87645\times {10}^{-3}$	$0.87645\times {10}^{-3}$
InvSqrt2	2	$0.87908\times {10}^{-3}$	$-0.87916\times {10}^{-3}$	$0.87914\times {10}^{-3}$	$0.87916\times {10}^{-3}$
InvSqrt3	3	$0.65007\times {10}^{-3}$	$-0.65017\times {10}^{-3}$	$0.65006·{10}^{-3}$	$0.65017\times {10}^{-3}$

(for codes with just one Newton–Raphson iteration) and

Table 3. Relative numerical errors for the second correction in the case of type float (compiler 32-bit). In the case of type double, the errors are equal to theoretical errors $\pm {\Delta }_{2\max }^{\left(i\right)}$ up to the accuracy given in the table.

Algorithm	i	${\mathbf{\Delta }}_{2\max }^{\left(\mathit{i}\right)}$	${\mathbf{\Delta }}_{2,\mathit{N}\min }^{\left(\mathit{i}\right)}$	${\mathbf{\Delta }}_{2,\mathit{N}\max }^{\left(\mathit{i}\right)}$	${\mathbf{\Delta }}_{2,\mathit{N}}^{\left(\mathit{i}\right)}$
InvSqrt	0	$4.59728\times {10}^{-6}$	$-4.65441\times {10}^{-6}$	$0.08336\times {10}^{-6}$	$4.65441\times {10}^{-6}$
InvSqrt1	1	$0.57617\times {10}^{-6}$	$-0.67207\times {10}^{-6}$	$0.64871\times {10}^{-6}$	$0.67207\times {10}^{-6}$
InvSqrt2	2	$0.57968\times {10}^{-6}$	$-0.64591\times {10}^{-6}$	$0.68363\times {10}^{-6}$	$0.68363\times {10}^{-6}$
InvSqrt3	3	$0.31694\times {10}^{-6}$	$-0.38701\times {10}^{-6}$	$0.35198\times {10}^{-6}$	$0.38701\times {10}^{-6}$

(the same codes but with two iterations).

Looking at the last column of Table 2 (this is the case of one iteration), we see that the code InvSqrt1 is slightly more accurate than InvSqrt2, and both are roughly almost two times more accurate than InvSqrt. However, it is the code InvSqrt3 that has the best accuracy. The computational costs of all these codes are practically the same (four multiplications in every case).

In the case of two iterations (Table 3), the code InvSqrt3 is the most accurate as well. Compared with InvSqrt, its accuracy is 12 times higher for single precision and 14.5 times higher for double precision. However, the computational costs of InvSqrt1 and InvSqrt3 (eight multiplications) are higher than the cost of InvSqrt (seven multiplications). Therefore, the code InvSqrt2 has some advantage, as it is less accurate than InvSqrt3 but cheaper. In the single-precision case the code InvSqrt2 is 6.8 times more accurate than InvSqrt.

We point out that the round-off errors in the single-precision case significantly decrease the gain of the accuracy of the new algorithms as compared with the theoretical values, especially in the case of two Newton–Raphson corrections (compare the third and the last column of Table 3).

The range of errors in the case of subnormal numbers (using the codes described in Section 3.1) is shown in Table 1. One can easily see that the relative errors are similar—in fact, even slightly lower—than in the case of normal numbers.

Although the original InvSqrt code used only one Newton–Raphson iteration, and in this paper, we focus mostly on two iterations, it is worthwhile to also briefly consider the case of more iterations. Then, the increased computational cost is accompanied by increased accuracy. We confine ourselves to the code InvSqrt2 (see the end of Section 3.2), which is less expensive than InvSqrt3 (and the advantage of InvSqrt2 increases with the number of iterations). In the double-precision case, the maximum error for three Newton–Raphson corrections is much lower, and the fourth correction yields the best possible accuracy.

$$\begin{array}{c}{\Delta }_{1D,N}^{\left(2\right)}=0.87908\times {10}^{-3},\hfill \\ {\Delta }_{2D,N}^{\left(2\right)}=0.57968\times {10}^{-6},\hfill \\ {\Delta }_{3D,N}^{\left(2\right)}=2.5213\times {10}^{-13},\hfill \\ {\Delta }_{4D,N}^{\left(2\right)}=1.1103\times {10}^{-16}.\hfill \end{array}$$

(52)

In the case of single precision, we already get the best possible accuracy for the third correction, given by adding the line y* = 1.5f - halfx*y*y as Line 8 in the code InvSqrt2 (see Section 2.2).

$$\begin{array}{c}{\Delta }_{1,N}^{\left(2\right)}=0.87916\times {10}^{-3},\hfill \\ {\Delta }_{2,N}^{\left(2\right)}=0.68363\times {10}^{-6},\hfill \\ {\Delta }_{3,N}^{\left(2\right)}=0.89367\times {10}^{-7}\hfill \end{array}$$

(53)

The derivation of all numerical codes presented in this paper did not take rounding errors into account. Therefore, the best floating-point parameters can be slightly different from the rounding of the best real parameters, all the more so since the distribution of the errors is still not exactly symmetric (compare fourth and fifth columns in Table 2 and Table 3). The full analysis of this problem is much more difficult than the analogous analysis for the original InvSqrt code because we now have several parameters to be optimized instead of a single magic constant. At the same time, the increase in accuracy is negligible. Actually, much greater differences in the accuracy appear in numerical experiments as a result of using different devices (see Appendix B).

As an example, we present the results of an experimental search in the case of the code InvSqrt3 with one Newton–Raphson correction (three parameters to be optimized). The modified Newton–Raphson coefficients are found to be

$$c_{11num}^{\left(3\right)}=1.681911588f,c_{21num}^{\left(3\right)}=k_1=0.7039490938f.$$

(54)

Figure 4 summarizes the last step of this analysis. The dependence of maximum errors on R shows clearly that the optimum value for the magic constant is slightly shifted as compared to the theoretical (real) value:

$$R_{num}^{\left(3\right)}=17+0x5f200000=0x5f200011.$$

(55)

The corresponding errors given by

$${\Delta }_{1,N\max }^{\left(3\right)}=6.50112284·{10}^{-4},{\Delta }_{1,N\min }^{\left(3\right)}=-6.501092575·{10}^{-4}$$

(56)

are nearly symmetric. They are smaller than the maximum error ${\Delta }_{1,N}^{\left(3\right)}$ corresponding to our theoretical values, but only by about $0.001\%$ (see Table 2).

5. Conclusions

We presented two new modifications (InvSqrt2 and $\mathit{InvSqrt}\mathit{3}$) of the fast inverse square root code in single-, double-, and higher-precision versions. Each code has its own magic constant. All new algorithms are much more accurate than the original code InvSqrt. One of the new algorithms, InvSqrt2, has the same computational cost as InvSqrt in the case of any precision. Another code, InvSqrt3, has the best accuracy, but is more expensive if the number of Newton–Raphson corrections is greater than 1. However, its gain in accuracy is very high, even by more than 12 times for two iterations (see Table 3 in Section 4).

Our approach was to modify the Newton–Raphson method by introducing arbitrary parameters, which are then determined by minimizing the maximum relative error. It is expected that such modifications will provide a significant increase in accuracy, especially in the case of asymmetric error distribution for Newton–Raphson corrections (and this is the case with the inverse square root function when these corrections are non-positive). One has to remember that due to rounding errors, our theoretical results may differ from the best floating-point parameters, but the difference is negligible (see the end of Section 4). In fact, parameters (magic constants and modified Newton–Raphson coefficients) from a certain range near the values obtained in this article seem equally good for all practical purposes.

Concerning potential applications, we have to acknowledge that for general-purpose computing, the SSE and AVX reciprocal square root instructions are faster and more accurate. We hope, however, that the proposed algorithms can be applied in embedded systems and microcontrollers without a hardware floating-point divider, and potentially in FPGAs. Moreover, in contrast to the SSE and AVX instructions, our approach can be easily extended to computational platforms of high precision, like 256-bit or 512-bit platforms.