Menu
Home
People
Places
Arts
History
Plants & Animals
Science
Life & Culture
Technology
Reference.org
Fast inverse square root
open-in-new
<h2 id="history">History</h2> <p><a href="/facts/William_Kahan/J1GNQDyN">William Kahan</a> and K.C. Ng at Berkeley wrote an unpublished paper in May 1986 describing how to calculate the square root using bit-fiddling techniques followed by Newton iterations.<a class="footnote-ref" id="fnref:4" href="#fn:4"><sup>4</sup></a> In the late 1980s, <a href="/facts/Cleve_Moler/o0wxumjp">Cleve Moler</a> at <a href="/facts/Ardent_Computer/Rhd2MEo2">Ardent Computer</a> learned about this technique<a class="footnote-ref" id="fnref:5" href="#fn:5"><sup>5</sup></a> and passed it along to his coworker Greg Walsh. Greg Walsh devised the now-famous constant and fast inverse square root algorithm. Gary Tarolli was consulting for Kubota, the company funding Ardent at the time, and likely brought the algorithm to <a href="/facts/3dfx_Interactive/8t0LNDBb">3dfx Interactive</a> circa 1994.<a class="footnote-ref" id="fnref:6" href="#fn:6"><sup>6</sup></a><a class="footnote-ref" id="fnref:7" href="#fn:7"><sup>7</sup></a> </p><p><a href="/facts/Jim_Blinn/UxlnCaWD">Jim Blinn</a> demonstrated a simple approximation of the inverse square root in a 1997 column for <i><a href="/facts/List_of_Institute_of_Electrical_and_Electronics_Engineers_publications/WDE73uBS">IEEE Computer Graphics and Applications</a></i>.<a class="footnote-ref" id="fnref:8" href="#fn:8"><sup>8</sup></a> Reverse engineering of other contemporary 3D video games uncovered a variation of the algorithm in <a href="/facts/Activision/NmDDD1wW">Activision</a>'s 1997 <a href="/facts/Interstate_%252776/GR3WAsuS">Interstate '76</a>.<a class="footnote-ref" id="fnref:9" href="#fn:9"><sup>9</sup></a> </p><p><i><a href="/facts/Quake_III_Arena/7WgcvDVt">Quake III Arena</a></i>, a <a href="/facts/First-person_shooter/LQQX6yAw">first-person shooter</a> video game, was released in 1999 by <a href="/facts/Id_Software/NsbHCMsV">id Software</a> and used the algorithm. Brian Hook may have brought the algorithm from 3dfx to id Software.<a class="footnote-ref" id="fnref:10" href="#fn:10"><sup>10</sup></a> A discussion of the code appeared on the Chinese developer forum CSDN in 2000,<a class="footnote-ref" id="fnref:11" href="#fn:11"><sup>11</sup></a> and <a href="/facts/Usenet/lYJWFBR6">Usenet</a> and the gamedev.net forum spread the code widely in 2002 and 2003.<a class="footnote-ref" id="fnref:12" href="#fn:12"><sup>12</sup></a> Speculation arose as to who wrote the algorithm and how the constant was derived; some guessed <a href="/facts/John_D._Carmack/PKRfqBQn">John Carmack</a>.<a class="footnote-ref" id="fnref:13" href="#fn:13"><sup>13</sup></a> <i>Quake III</i>'s full source code was released at <a href="/facts/QuakeCon_2005/R4vS0a2R">QuakeCon 2005</a>, but provided no answers. The authorship question was resolved in 2006 when Greg Walsh, the original author, contacted Beyond3D after their speculation gained popularity on <a href="/facts/Slashdot/vVpMuby2">Slashdot</a>.<a class="footnote-ref" id="fnref:14" href="#fn:14"><sup>14</sup></a> </p><p>In 2007 the <a href="/facts/Algorithm/fnl5NmRt">algorithm</a> was implemented in some dedicated hardware <a href="/facts/Vertex_shader/qfhkKFwA">vertex shaders</a> using <a href="/facts/Field-programmable_gate_array/M8LSGcAU">field-programmable gate arrays</a> (FPGA).<a class="footnote-ref" id="fnref:15" href="#fn:15"><sup>15</sup></a><a class="footnote-ref" id="fnref:16" href="#fn:16"><sup>16</sup></a> </p> <h2 id="motivation">Motivation</h2> <p>The inverse square root of a floating point number is used in <a href="/facts/Digital_signal_processing/PEG7azr8">digital signal processing</a> to normalize a vector, scaling it to <a href="/facts/Length/EwFyAqGH">length</a> 1 to produce a <a href="/facts/Unit_vector/GyaXneDn">unit vector</a>.<a class="footnote-ref" id="fnref:17" href="#fn:17"><sup>17</sup></a> For example, <a href="/facts/Computer_graphics/lapBE8wv">computer graphics</a> programs use inverse <a href="/facts/Square_root/AfrzfBdQ">square roots</a> to compute <a href="/facts/Angle_of_incidence_(optics)/DTaCjhTU">angles of incidence</a> and <a href="/facts/Reflection_(computer_graphics)/ThoxJKIy">reflection</a> for <a href="/facts/Computer_graphics_lighting/kHshvyKY">lighting</a> and <a href="/facts/Shading/ewcDEfda">shading</a>. <a href="/facts/3D_graphics/zlWW6eaF">3D graphics</a> programs must perform millions of these calculations every second to simulate lighting. When the code was developed in the early 1990s, most floating point processing power lagged the speed of integer processing.<a class="footnote-ref" id="fnref:18" href="#fn:18"><sup>18</sup></a> This was troublesome for 3D graphics programs before the advent of specialized hardware to handle <a href="/facts/Transform%2c_clipping%2c_and_lighting/wnDkGH4W">transform and lighting</a>. Computation of square roots usually depends upon many division operations, which for floating point numbers are <a href="/facts/Computationally_expensive/pStoBgBn">computationally expensive</a>. The fast inverse square generates a good approximation with only one division step. </p><p>The length of the <a href="/facts/Euclidean_vector/bCLrdksy">vector</a> is determined by calculating its <a href="/facts/Euclidean_norm/R2UbzmzM">Euclidean norm</a>: the square root of the sum of squares of the <a href="/facts/Vector_components/bCLrdksy">vector components</a>. When each component of the vector is divided by that length, the new vector will be a unit vector pointing in the same direction. In a 3D graphics program, all vectors are in three-<a href="/facts/Dimension_(vector_space)/z8m02A3W">dimensional</a> space, so v {\displaystyle {\boldsymbol {v}}} would be a vector ( v 1 , v 2 , v 3 ) {\displaystyle (v_{1},v_{2},v_{3})} . Then, </p> ‖ v ‖ = v 1 2 + v 2 2 + v 3 2 {\displaystyle \|{\boldsymbol {v}}\|={\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}} <p>is the Euclidean norm of the vector, and the normalized (unit) vector is </p> v ^ = 1 ‖ v ‖ v = 1 v 1 2 + v 2 2 + v 3 2 v , {\displaystyle {\begin{aligned}{\boldsymbol {\hat {v}}}&={\frac {1}{\left\|{\boldsymbol {v}}\right\|}}{\boldsymbol {v}}\\&={\frac {1}{\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}}\,{\boldsymbol {v}},\end{aligned}}} <p>where the fraction term is the inverse square root of v 1 2 + v 2 2 + v 3 2 {\displaystyle v_{1}^{2}+v_{2}^{2}+v_{3}^{2}} . </p><p>At the time, floating-point division was generally expensive compared to multiplication; the fast inverse square root algorithm bypassed the division step, giving it its performance advantage. </p> <h2 id="overview-of-the-code">Overview of the code</h2> <p>The following <a href="/facts/C_(programming_language)/Ky2No763">C</a> code is the fast inverse square root implementation from <i><a href="/facts/Quake_III_Arena/7WgcvDVt">Quake III Arena</a></i>, stripped of <a href="/facts/C_preprocessor/KuhuxKYP">C preprocessor</a> directives, but including the exact original comment text:<a class="footnote-ref" id="fnref:19" href="#fn:19"><sup>19</sup></a> </p> float Q_rsqrt( float number ) { long i; float x2, y; const float threehalfs = 1.5F; x2 = number * 0.5F; y = number; i = * ( long * ) &y; // evil floating point bit level hacking i = 0x5f3759df - ( i >> 1 ); // what the fuck? y = * ( float * ) &i; y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration // y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed return y; } <p>At the time, the general method to compute the inverse square root was to calculate an approximation for 1 x {\textstyle {\frac {1}{\sqrt {x}}}} , then revise that approximation via another method until it came within an acceptable error range of the actual result. <a href="/facts/Methods_of_computing_square_roots/dRG2BtBp">Common software methods</a> in the early 1990s drew approximations from a <a href="/facts/Lookup_table/ULPEAGTJ">lookup table</a>.<a class="footnote-ref" id="fnref:20" href="#fn:20"><sup>20</sup></a> The key of the fast inverse square root was to directly compute an approximation by utilizing the structure of floating-point numbers, proving faster than table lookups. The algorithm was approximately four times faster than computing the square root with another method and calculating the reciprocal via floating-point division.<a class="footnote-ref" id="fnref:21" href="#fn:21"><sup>21</sup></a> The algorithm was designed with the <a href="/facts/IEEE_754-1985/5uWBQrc0">IEEE 754-1985</a> 32-bit floating-point specification in mind, but investigation from Chris Lomont showed that it could be implemented in other floating-point specifications.<a class="footnote-ref" id="fnref:22" href="#fn:22"><sup>22</sup></a> </p><p>The advantages in speed offered by the fast inverse square root trick came from treating the 32-bit floating-point word<a class="footnote-ref" id="fnref:23" href="#fn:23"><sup>23</sup></a> as an <a href="/facts/Integer_(computer_science)/zWabo5M9">integer</a>, then subtracting it from a "<a href="/facts/Magic_number_(programming)/Xg7mVI7V">magic</a>" constant, <a href="/facts/Hexadecimal/02ohQW22">0x</a>5F3759DF.<a class="footnote-ref" id="fnref:24" href="#fn:24"><sup>24</sup></a><a class="footnote-ref" id="fnref:25" href="#fn:25"><sup>25</sup></a><a class="footnote-ref" id="fnref:26" href="#fn:26"><sup>26</sup></a><a class="footnote-ref" id="fnref:27" href="#fn:27"><sup>27</sup></a> This integer subtraction and <a href="/facts/Bitwise_operation/0gkqnUda">bit shift</a> results in a bit pattern which, when <a href="/facts/Type_conversion/FuPD945T">re-defined</a> as a floating-point number, is a rough approximation for the inverse square root of the number. One iteration of Newton's method is performed to gain some accuracy, and the code is finished. The algorithm generates reasonably accurate results using a unique first approximation for <a href="/facts/Newton%2527s_method/VlRhI2FL">Newton's method</a>; however, it is much slower and less accurate than using the <a href="/facts/Streaming_SIMD_Extensions/mMa6EFyB">SSE</a> instruction rsqrtss on x86 processors also released in 1999.<a class="footnote-ref" id="fnref:28" href="#fn:28"><sup>28</sup></a><a class="footnote-ref" id="fnref:29" href="#fn:29"><sup>29</sup></a> </p> <h3>Worked example</h3> <p>As an example, the number x = 0.15625 {\displaystyle x=0.15625} can be used to calculate 1 x ≈ 2.52982 {\textstyle {\frac {1}{\sqrt {x}}}\approx 2.52982} . The first steps of the algorithm are illustrated below: </p> 0011_1110_0010_0000_0000_0000_0000_0000 Bit pattern of both x and i 0001_1111_0001_0000_0000_0000_0000_0000 Shift right one position: (i >> 1) 0101_1111_0011_0111_0101_1001_1101_1111 The magic number 0x5F3759DF 0100_0000_0010_0111_0101_1001_1101_1111 The result of 0x5F3759DF - (i >> 1) <p>Interpreting as <a href="/facts/IEEE_754/dD3e7zAl">IEEE</a> 32-bit representation: </p> 0_01111100_01000000000000000000000 1.25 × 2−3 0_00111110_00100000000000000000000 1.125 × 2−65 0_10111110_01101110101100111011111 1.432430... × 263 0_10000000_01001110101100111011111 1.307430... × 21 <p>Reinterpreting this last bit pattern as a floating point number gives the approximation y = 2.61486 {\displaystyle y=2.61486} , which has an error of about 3.4%. After one iteration of <a href="/facts/Newton%2527s_method/VlRhI2FL">Newton's method</a>, the final result is y = 2.52549 {\displaystyle y=2.52549} , an error of only 0.17%. </p> <h3>Avoiding undefined behavior</h3> <p class="note">Main article: <a href="/facts/Type_punning/xYtejQnM">Type punning</a></p> <p>According to the C standard, reinterpreting a floating point value as an integer by casting then dereferencing the pointer to it is not valid (<a href="/facts/Undefined_behavior/awUm2vGd">undefined behavior</a>).<a class="footnote-ref" id="fnref:30" href="#fn:30"><sup>30</sup></a> This can be avoided by using alternative <a href="/facts/Type_punning/xYtejQnM">type punning</a> techniques such as C's unions<a class="footnote-ref" id="fnref:31" href="#fn:31"><sup>31</sup></a> or <a href="/facts/C%252B%252B20/hwpB1Vzp">C++20</a>'s std::bit_cast.<a class="footnote-ref" id="fnref:32" href="#fn:32"><sup>32</sup></a> </p> <h2 id="algorithm">Algorithm</h2> <p>The algorithm computes 1 x {\textstyle {\frac {1}{\sqrt {x}}}} by performing the following steps: </p> <ol><li>Alias the argument x {\displaystyle x} to an integer as a way to compute an approximation of the <a href="/facts/Binary_logarithm/Uy7kAtRm">binary logarithm</a> log 2 ( x ) {\textstyle \log _{2}(x)} </li> <li>Use this approximation to compute an approximation of log 2 ( 1 x ) = − 1 2 log 2 ( x ) {\textstyle \log _{2}\left({\frac {1}{\sqrt {x}}}\right)=-{\frac {1}{2}}\log _{2}(x)} </li> <li>Alias back to a float, as a way to compute an approximation of the base-2 exponential</li> <li>Refine the approximation using a single iteration of <a href="/facts/Newton%2527s_method/VlRhI2FL">Newton's method</a>.</li></ol> <h3>Floating-point representation</h3> <p class="note">Main article: <a href="/facts/Single-precision_floating-point_format/e6Fez7BM">Single-precision floating-point format</a></p> <p>Since this <a href="/facts/Algorithm/fnl5NmRt">algorithm</a> relies heavily on the bit-level representation of single-precision floating-point numbers, a short overview of this representation is provided here. To encode a non-zero <a href="/facts/Real_number/R02gw5Pb">real number</a> x {\displaystyle x} as a single precision float, the first step is to write x {\displaystyle x} as a <a href="/facts/Normalized_number/auQAhM4D">normalized</a> binary number:<a class="footnote-ref" id="fnref:33" href="#fn:33"><sup>33</sup></a> </p> x = ± 1. b 1 b 2 b 3 … × 2 e x {\displaystyle {\begin{aligned}x&=\pm 1.b_{1}b_{2}b_{3}\ldots \times 2^{e_{x}}\end{aligned}}} <p>where the exponent e x {\textstyle e_{x}} is an integer, and 1. b 1 b 2 b 3 … {\textstyle 1.b_{1}b_{2}b_{3}\ldots } is the binary representation of the <a href="/facts/Significand/JFCVyTz4">significand</a>. Since the single bit before the point in the significand is always 1, it does not need be stored. The equation can be rewritten as: </p> x = ( − 1 ) S x ⋅ 2 e x ( 1 + m x ) {\displaystyle {\begin{aligned}x&=(-1)^{S_{x}}\cdot 2^{e_{x}}(1+m_{x})\end{aligned}}} <p>where m x {\textstyle m_{x}} means 0. b 1 b 2 b 3 … {\textstyle 0.b_{1}b_{2}b_{3}\ldots } , so m x ∈ [ 0 , 1 ) {\textstyle m_{x}\in [0,1)} . From this form, three unsigned integers are computed:<a class="footnote-ref" id="fnref:34" href="#fn:34"><sup>34</sup></a> </p> <ul><li> S x {\textstyle S_{x}} , the "sign bit", is 0 {\displaystyle 0} if x {\displaystyle x} is positive and 1 {\displaystyle 1} negative or zero (1 bit)</li> <li> E x = e x + B {\textstyle E_{x}=e_{x}+B} is the "biased exponent", where B = 127 {\displaystyle B=127} is the "<a href="/facts/Exponent_bias/u5MIKgKP">exponent bias</a>"<a class="footnote-ref" id="fnref:35" href="#fn:35"><sup>35</sup></a> (8 bits)</li> <li> M x = m x × L {\textstyle M_{x}=m_{x}\times L} , where L = 2 23 {\textstyle L=2^{23}} <a class="footnote-ref" id="fnref:36" href="#fn:36"><sup>36</sup></a> (23 bits)</li></ul> <p>Thus: e x = E x − B {\textstyle e_{x}=E_{x}-B} and m x = M x L {\textstyle m_{x}={\frac {M_{x}}{L}}} . </p><p>These fields are then packed, left to right, into a 32-bit container.<a class="footnote-ref" id="fnref:37" href="#fn:37"><sup>37</sup></a> </p><p>As an example, consider again the number x = 0.15625 = 0.00101 2 {\displaystyle x=0.15625=0.00101_{2}} . Normalizing x {\displaystyle x} yields: </p> x = ( − 1 ) 0 ⋅ 2 − 3 ( 1 + 0.25 ) = + 2 − 3 ( 1 + 0.25 ) {\textstyle x=(-1)^{0}\cdot 2^{-3}(1+0.25)=+2^{-3}(1+0.25)} <p>and thus, the three unsigned integer fields are: </p> <ul><li> S = 0 {\displaystyle S=0} </li> <li> E = − 3 + 127 = 124 = 0111 1100 2 {\displaystyle E=-3+127=124=0111\ 1100_{2}} </li> <li> M = 0.25 × 2 23 = 2 097 152 = 0010 0000 0000 0000 0000 0000 2 {\displaystyle M=0.25\times 2^{23}=2\ 097\ 152=0010\ 0000\ 0000\ 0000\ 0000\ 0000_{2}} </li></ul> <p>these fields are packed as shown in the figure below: </p> <p>The number is represented in <a href="/facts/Binary_number/a6JDSRPs">binary</a> as: I x = S x ⋅ 2 31 + E x L + M x {\displaystyle I_{x}=S_{x}\cdot 2^{31}+E_{x}L+M_{x}} </p><p>Also, since this algorithm works on real numbers, x {\textstyle {\sqrt {x}}} is only defined for x ≥ 0 {\textstyle x\geq 0} . The code thus assumes x ≥ 0 {\textstyle x\geq 0} and S x = 0 {\displaystyle S_{x}=0} . </p><p>The number, given to calculate the square root, could be rewritten as: </p> I x = E x L + M x {\textstyle I_{x}=E_{x}L+M_{x}} x = ( − 1 ) 0 ⋅ 2 e x ( 1 + m x ) = + 2 e x ( 1 + m x ) {\textstyle x=(-1)^{0}\cdot 2^{e_{x}}(1+m_{x})=+2^{e_{x}}(1+m_{x})} <h3>Aliasing to an integer as an approximate logarithm</h3> <p>If 1 x {\textstyle {\frac {1}{\sqrt {x}}}} were to be calculated without a computer or a calculator, a <a href="/facts/Table_of_logarithms/uCAtPzJo">table of logarithms</a> would be useful, together with the identity log b ( 1 x ) = log b ( x − 1 2 ) = − 1 2 log b ( x ) {\textstyle \log _{b}\left({\frac {1}{\sqrt {x}}}\right)=\log _{b}\left(x^{-{\frac {1}{2}}}\right)=-{\frac {1}{2}}\log _{b}(x)} , which is valid for every base b {\displaystyle b} . The fast inverse square root is based on this identity, and on the fact that aliasing a float32 to an integer gives a rough approximation of its logarithm. Here is how: </p><p>If x {\displaystyle x} is a positive <a href="/facts/Normal_number_(computing)/5pvlHVvg">normal number</a>: </p> x = 2 e x ( 1 + m x ) {\displaystyle x=2^{e_{x}}(1+m_{x})} <p>then </p> log 2 ( x ) = e x + log 2 ( 1 + m x ) {\displaystyle \log _{2}(x)=e_{x}+\log _{2}(1+m_{x})} <p>and since m x ∈ [ 0 , 1 ) {\displaystyle m_{x}\in [0,1)} , the logarithm on the right-hand side can be approximated by<a class="footnote-ref" id="fnref:38" href="#fn:38"><sup>38</sup></a> </p> log 2 ( 1 + m x ) ≈ m x + σ {\displaystyle \log _{2}(1+m_{x})\approx m_{x}+\sigma } <p>where σ {\displaystyle \sigma } is a free parameter used to tune the approximation. For example, σ = 0 {\displaystyle \sigma =0} yields exact results at both ends of the interval, while σ = 1 2 − 1 + ln ( ln ( 2 ) ) 2 ln ( 2 ) ≈ 0.0430357 {\textstyle \sigma ={\frac {1}{2}}-{\frac {1+\ln(\ln(2))}{2\ln(2)}}\approx 0.0430357} yields the <a href="/facts/Approximation_theory/87AaPXLp">optimal approximation</a> (the best in the sense of the <a href="/facts/Uniform_norm/5yIkjNE3">uniform norm</a> of the error). However, this value is not used by the algorithm as it does not take subsequent steps into account. </p> <p>Thus there is the approximation </p> log 2 ( x ) ≈ e x + m x + σ . {\displaystyle \log _{2}(x)\approx e_{x}+m_{x}+\sigma .} <p>Interpreting the floating-point bit-pattern of x {\displaystyle x} as an integer I x {\displaystyle I_{x}} yields<a class="footnote-ref" id="fnref:39" href="#fn:39"><sup>39</sup></a> </p> I x = E x L + M x = L ( e x + B + m x ) = L ( e x + m x + σ + B − σ ) ≈ L log 2 ( x ) + L ( B − σ ) . {\displaystyle {\begin{aligned}I_{x}&=E_{x}L+M_{x}\\&=L(e_{x}+B+m_{x})\\&=L(e_{x}+m_{x}+\sigma +B-\sigma )\\&\approx L\log _{2}(x)+L(B-\sigma ).\end{aligned}}} <p>It then appears that I x {\displaystyle I_{x}} is a scaled and shifted piecewise-linear approximation of log 2 ( x ) {\displaystyle \log _{2}(x)} , as illustrated in the figure on the right. In other words, log 2 ( x ) {\displaystyle \log _{2}(x)} is approximated by </p> log 2 ( x ) ≈ I x L − ( B − σ ) . {\displaystyle \log _{2}(x)\approx {\frac {I_{x}}{L}}-(B-\sigma ).} <h3>First approximation of the result</h3> <p>The calculation of y = 1 x {\textstyle y={\frac {1}{\sqrt {x}}}} is based on the identity </p> log 2 ( y ) = − 1 2 log 2 ( x ) {\displaystyle \log _{2}(y)=-{\tfrac {1}{2}}\log _{2}(x)} <p>Using the approximation of the logarithm above, applied to both x {\displaystyle x} and y {\displaystyle y} , the above equation gives: </p> I y L − ( B − σ ) ≈ − 1 2 ( I x L − ( B − σ ) ) {\displaystyle {\frac {I_{y}}{L}}-(B-\sigma )\approx -{\frac {1}{2}}\left({\frac {I_{x}}{L}}-(B-\sigma )\right)} <p>Thus, an approximation of I y {\displaystyle I_{y}} is: </p> I y ≈ 3 2 L ( B − σ ) − 1 2 I x {\displaystyle I_{y}\approx {\tfrac {3}{2}}L(B-\sigma )-{\tfrac {1}{2}}I_{x}} <p>which is written in the code as </p> i = 0x5f3759df - ( i >> 1 ); <p>The first term above is the magic number </p> 3 2 L ( B − σ ) = 0 x 5 F 3759 D F {\displaystyle {\tfrac {3}{2}}L(B-\sigma )={\mathtt {0x5F3759DF}}} <p>from which it can be inferred that σ ≈ 0.0450466 {\displaystyle \sigma \approx 0.0450466} . The second term, 1 2 I x {\displaystyle {\frac {1}{2}}I_{x}} , is calculated by shifting the bits of I x {\displaystyle I_{x}} one position to the right.<a class="footnote-ref" id="fnref:40" href="#fn:40"><sup>40</sup></a> </p> <h3>Newton's method</h3> <p class="note">Main article: <a href="/facts/Newton%2527s_method/VlRhI2FL">Newton's method</a></p> <p>The number y = 1 x {\displaystyle y={\tfrac {1}{\sqrt {x}}}} is a solution of the equation 1 y 2 − x = 0 {\displaystyle {\tfrac {1}{y^{2}}}-x=0} . The approximation yielded by the earlier steps can be refined by using a <a href="/facts/Root-finding_method/IxTDkhBq">root-finding method</a>, a method that finds the <a href="/facts/Zero_of_a_function/mlK3dhoT">zero of a function</a>. The algorithm uses <a href="/facts/Newton%2527s_method/VlRhI2FL">Newton's method</a>: if there is an approximation, y n {\displaystyle y_{n}} for y {\displaystyle y} , then a better approximation y n + 1 {\displaystyle y_{n+1}} can be calculated by taking y n − f ( y n ) f ′ ( y n ) {\displaystyle y_{n}-{\tfrac {f(y_{n})}{f'(y_{n})}}} , where f ′ ( y n ) {\displaystyle f'(y_{n})} is the <a href="/facts/Derivative/7Ito70Dh">derivative</a> of f ( y ) {\displaystyle f(y)} at y n {\displaystyle y_{n}} .<a class="footnote-ref" id="fnref:41" href="#fn:41"><sup>41</sup></a> Applied to the equation 1 y 2 − x = 0 {\displaystyle {\tfrac {1}{y^{2}}}-x=0} , Newton's method gives </p> f ( y ) = 1 y 2 − x f ′ ( y ) = − 2 y 3 y n + 1 = y n − f ( y n ) f ′ ( y n ) = y n + y n 3 2 ( 1 y n 2 − x ) = y n ( 3 2 − x 2 y n 2 ) , {\displaystyle {\begin{aligned}f(y)&={\frac {1}{y^{2}}}-x\\f'(y)&=-{\frac {2}{y^{3}}}\\y_{n+1}&=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}\\&=y_{n}+{\frac {y_{n}^{3}}{2}}\left({\frac {1}{y_{n}^{2}}}-x\right)\\&=y_{n}\left({\frac {3}{2}}-{\frac {x}{2}}y_{n}^{2}\right),\end{aligned}}} <p>which is written in the code as y = y * ( threehalfs - ( x2 * y * y ) );. </p><p>By repeating this step, using the output of the function ( y n + 1 {\displaystyle y_{n+1}} ) as the input of the next iteration, the algorithm causes y {\displaystyle y} to <a href="/facts/Rate_of_convergence/JVGuzPoS">converge</a> to the inverse square root.<a class="footnote-ref" id="fnref:42" href="#fn:42"><sup>42</sup></a> For the purposes of the <a href="/facts/Id_Tech_3/xLmklUPJ"><i>Quake III</i> engine</a>, only one iteration was used. A second iteration remained in the code but was <a href="/facts/Comment_out/Q8hHFq8h">commented out</a>.<a class="footnote-ref" id="fnref:43" href="#fn:43"><sup>43</sup></a> </p> <h3>Accuracy</h3> <p>As noted above, the approximation is very accurate. The single graph on the right plots the error of the function (that is, the error of the approximation after it has been improved by running one iteration of Newton's method), for inputs starting at 0.01, where the standard library gives 10.0 as a result, and InvSqrt() gives 9.982522, making the relative difference 0.0017478, or 0.175% of the true value, 10. The absolute error only drops from then on, and the relative error stays within the same bounds across all orders of magnitude. </p> <h2 id="subsequent-improvements">Subsequent improvements</h2> <h3>Magic number</h3> <p>It is not known precisely how the exact value for the magic number was determined. Chris Lomont developed a function to minimize <a href="/facts/Approximation_error/orSvvIMC">approximation error</a> by choosing the magic number R {\displaystyle R} over a range. He first computed the optimal constant for the linear approximation step as 0x5F37642F, close to 0x5F3759DF, but this new constant gave slightly less accuracy after one iteration of Newton's method.<a class="footnote-ref" id="fnref:44" href="#fn:44"><sup>44</sup></a> Lomont then searched for a constant optimal even after one and two Newton iterations and found 0x5F375A86, which is more accurate than the original at every iteration stage.<a class="footnote-ref" id="fnref:45" href="#fn:45"><sup>45</sup></a> He concluded by asking whether the exact value of the original constant was chosen through derivation or <a href="/facts/Trial_and_error/NGm2GoEd">trial and error</a>.<a class="footnote-ref" id="fnref:46" href="#fn:46"><sup>46</sup></a> Lomont said that the magic number for 64-bit IEEE754 size type double is 0x5FE6EC85E7DE30DA, but it was later shown by Matthew Robertson to be exactly 0x5FE6EB50C7B537A9.<a class="footnote-ref" id="fnref:47" href="#fn:47"><sup>47</sup></a> </p><p>Jan Kadlec reduced the relative error by a further factor of 2.7 by adjusting the constants in the single Newton's method iteration as well,<a class="footnote-ref" id="fnref:48" href="#fn:48"><sup>48</sup></a> arriving after an exhaustive search at </p> conv.i = 0x5F1FFFF9 - ( conv.i >> 1 ); conv.f *= 0.703952253f * ( 2.38924456f - x * conv.f * conv.f ); return conv.f; <p>A complete mathematical analysis for determining the magic number is now available for single-precision floating-point numbers.<a class="footnote-ref" id="fnref:49" href="#fn:49"><sup>49</sup></a><a class="footnote-ref" id="fnref:50" href="#fn:50"><sup>50</sup></a> </p> <h3>Zero finding</h3> <p>Intermediate to the use of one vs. two iterations of Newton's method in terms of speed and accuracy is a single iteration of <a href="/facts/Halley%2527s_method/t0tcGcHA">Halley's method</a>. In this case, Halley's method is equivalent to applying Newton's method with the starting formula f ( y ) = 1 y 1 / 2 − x y 3 / 2 = 0 {\displaystyle f(y)={\frac {1}{y^{1/2}}}-xy^{3/2}=0} . The update step is then </p> y n + 1 = y n − f ( y n ) f ′ ( y n ) = y n ( 3 + x y n 2 1 + 3 x y n 2 ) , {\displaystyle y_{n+1}=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}=y_{n}\left({\frac {3+xy_{n}^{2}}{1+3xy_{n}^{2}}}\right),} <p>where the implementation should calculate x y n 2 {\displaystyle xy_{n}^{2}} only once, via a temporary variable. </p> <h3>Obsolescence</h3> <p>Subsequent additions by hardware manufacturers have made this algorithm redundant for the most part. For example, on <a href="/facts/X86/Kypq5fgu">x86</a>, Intel introduced the <a href="/facts/Streaming_SIMD_Extensions/mMa6EFyB">SSE</a> instruction rsqrtss in 1999. In a 2009 benchmark on the <a href="/facts/Intel_Core_2/9gJP9o8p">Intel Core 2</a>, this instruction took 0.85ns per float compared to 3.54ns for the fast inverse square root algorithm, and had less error.<a class="footnote-ref" id="fnref:51" href="#fn:51"><sup>51</sup></a> </p><p>Some low-cost embedded systems do not have specialized square root instructions. However, manufacturers of these systems usually provide trigonometric and other math libraries, based on algorithms such as <a href="/facts/CORDIC/sQqTjkd2">CORDIC</a>. </p> <h2 id="see-also">See also</h2> <ul><li><a href="/facts/Methods_of_computing_square_roots/dRG2BtBp">Methods of computing square roots § Approximations that depend on the floating point representation</a></li> <li><a href="/facts/Magic_number_(programming)/Xg7mVI7V">Magic number</a></li></ul> <h2 id="notes">Notes</h2> <h3>Bibliography</h3> <ul><li>Blinn, Jim (July 1997). "Floating Point Tricks". <i>IEEE Computer Graphics & Applications</i>. 17 (4): 80. <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1109%2F38.595279">10.1109/38.595279</a>.</li> <li>Blinn, Jim (2003). <i>Jim Blinn's Corner: Notation, notation notation</i>. Morgan Kaufmann. <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 1-55860-860-5.</li> <li>Eberly, David (2001). <a href="https://archive.org/details/isbn_9781558605930"><i>3D Game Engine Design</i></a>. Morgan Kaufmann. <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 978-1-55860-593-0.</li> <li>Goldberg, David (1991). "What every computer scientist should know about floating-point arithmetic". <i>ACM Computing Surveys</i>. 23 (1): 5–48. <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1145%2F103162.103163">10.1145/103162.103163</a>. <a href="/facts/S2CID_(identifier)/ldJsHa2Y">S2CID</a> <a href="https://api.semanticscholar.org/CorpusID:222008826">222008826</a>.</li> <li>Hardy, Godfrey (1908). <a href="http://www.gutenberg.org/files/38769"><i>A Course of Pure Mathematics</i></a>. <a href="/facts/Cambridge/gNRdbmwK">Cambridge, UK</a>: <a href="/facts/Cambridge_University_Press/fgEBSSRq">Cambridge University Press</a>. Centenary Edition, 2008, <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 0-521-72055-9.</li> <li>Hennessey, John; Patterson, David A. (1998). <a href="https://archive.org/details/computerorganiz000henn"><i>Computer Organization and Design</i></a> (2nd ed.). San Francisco, CA: Morgan Kaufmann Publishers. <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 978-1-55860-491-9.</li> <li>Lomont, Chris (February 2003). <a href="http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf">"Fast Inverse Square Root"</a> (PDF). Retrieved 2009-02-13.</li> <li>McEniry, Charles (August 2007). <a href="https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf">"The Mathematics Behind the Fast Inverse Square Root Function Code"</a> (PDF). Archived from <a href="http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf">the original</a> (PDF) on 2015-05-11.</li> <li>Middendorf, Lars; Mühlbauer, Felix; Umlauf, George; Bodba, Christophe (June 1, 2007). <a href="https://link.springer.com/content/pdf/10.1007/978-0-387-72258-0_14.pdf">"Embedded Vertex Shader in FPGA"</a> (PDF). In Rettberg, Achin (ed.). <i>Embedded System Design: Topics, Techniques and Trends</i>. <a href="/facts/International_Federation_for_Information_Processing/RgTMBVN9">IFIP</a> TC10 Working Conference:International Embedded Systems Symposium (IESS). et al. <a href="/facts/Irvine%2c_California/Sa2RCCIU">Irvine, California</a>: Springer. <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1007%2F978-0-387-72258-0_14">10.1007/978-0-387-72258-0_14</a>. <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 978-0-387-72257-3. <a href="https://web.archive.org/web/20190501235316/https://link.springer.com/content/pdf/10.1007/978-0-387-72258-0_14.pdf">Archived</a> (PDF) from the original on 2019-05-01. <a href="/facts/Open_access/r7rLYN6C"></a></li> <li>Striegel, Jason (2008-12-04). <a href="https://web.archive.org/web/20090215020337/http://www.hackszine.com/blog/archive/2008/12/quakes_fast_inverse_square_roo.html">"Quake's fast inverse square root"</a>. <i>Hackszine</i>. <a href="/facts/O%2527Reilly_Media/Fr3Kv7sm">O'Reilly Media</a>. Archived from <a href="http://www.hackszine.com/blog/archive/2008/12/quakes_fast_inverse_square_roo.html">the original</a> on 2009-02-15. Retrieved 2013-01-07.</li> <li><a href="/facts/IEEE_Computer_Society/asUJfKHp">IEEE Computer Society</a> (1985). <a href="https://standards.ieee.org/ieee/754/993/">"754-1985 - IEEE Standard for Binary Floating-Point Arithmetic"</a>. <a href="/facts/Institute_of_Electrical_and_Electronics_Engineers/UddC4pUW">Institute of Electrical and Electronics Engineers</a>.</li> <li>Moroz, Leonid V.; Walczyk, Cezary J.; Hrynchyshyn, Andriy; Holimath, Vijay; Cieslinski, Jan L. (January 2018). "Fast calculation of inverse square root with the use of magic constant analytical approach". <i>Applied Mathematics and Computation</i>. 316 (C). Elsevier Science Inc.: 245–255. <a href="/facts/ArXiv_(identifier)/H6EtgnBe">arXiv</a>:<a href="https://arxiv.org/abs/1603.04483">1603.04483</a>. <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1016%2Fj.amc.2017.08.025">10.1016/j.amc.2017.08.025</a>. <a href="/facts/S2CID_(identifier)/ldJsHa2Y">S2CID</a> <a href="https://api.semanticscholar.org/CorpusID:7494112">7494112</a>.</li></ul> <h2 id="further-reading">Further reading</h2> <ul><li>Kushner, David (August 2002). "The wizardry of Id". <i>IEEE Spectrum</i>. 39 (8): 42–47. <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1109%2FMSPEC.2002.1021943">10.1109/MSPEC.2002.1021943</a>.</li></ul> <h2 id="external-links">External links</h2> <ul><li><a href="https://web.archive.org/web/20120920120948/http://blog.quenta.org/2012/09/0x5f3759df.html">0x5f3759df</a>, further investigations into accuracy and generalizability of the algorithm by Christian Plesner Hansen</li> <li><a href="https://www.beyond3d.com/content/articles/8/">Origin of Quake3's Fast InvSqrt()</a></li> <li><a href="https://github.com/id-Software/Quake-III-Arena">Quake III Arena source code</a> (<a href="https://archive.softwareheritage.org/swh:1:dir:c6f07c2173a458d098de45d4c459a8f1916d900f;origin=https://github.com/id-Software/Quake-III-Arena;visit=swh:1:snp:4ab9bcef131aaf449a7c01370aff8c91dcecbf5f;anchor=swh:1:rev:dbe4ddb10315479fc00086f08e25d968b4b43c49">archived in Software Heritage</a>)</li> <li><a href="https://www.desmos.com/calculator/yoz6n1wlvu">Implementation of InvSqrt</a> in <a href="http://www.desmos.com">DESMOS</a></li> <li><a href="https://www.youtube.com/watch?v=p8u_k2LIZyo">"Fast Inverse Square Root — A Quake III Algorithm"</a> (<a href="/facts/YouTube/db276jst">YouTube</a>)</li></ul> <h2 id="references">References</h2> <ol> <li id="fn:1"><p>Ruskin, Elan (2009-10-16). "Timing square root". Some Assembly Required. Archived from the original on 2021-02-08. Retrieved 2015-05-07. <a href="https://web.archive.org/web/20210208132927/http://assemblyrequired.crashworks.org/timing-square-root/" target="_blank">https://web.archive.org/web/20210208132927/http://assemblyrequired.crashworks.org/timing-square-root/</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></p></li> <li id="fn:2"><p>feilipu. "z88dk is a collection of software development tools that targets the 8080 and z80 computers". GitHub. <a href="https://github.com/z88dk/z88dk/tree/master/libsrc/_DEVELOPMENT/math/float/math32#sqrt-and-invsqrt" target="_blank">https://github.com/z88dk/z88dk/tree/master/libsrc/_DEVELOPMENT/math/float/math32#sqrt-and-invsqrt</a> <a href="#fnref:2" class="footnote-back-ref">↩</a></p></li> <li id="fn:3"><p>Munafo, Robert. "Notable Properties of Specific Numbers". mrob.com. Archived from the original on 16 November 2018. <a href="https://mrob.com/pub/math/numbers-16.html#le009_16" target="_blank">https://mrob.com/pub/math/numbers-16.html#le009_16</a> <a href="#fnref:3" class="footnote-back-ref">↩</a></p></li> <li id="fn:4"><p>"sqrt implementation in fdlibm - See W. Kahan and K.C. Ng's discussion in comments in lower half of this code". <a href="http://www.netlib.org/fdlibm/e_sqrt.c" target="_blank">http://www.netlib.org/fdlibm/e_sqrt.c</a> <a href="#fnref:4" class="footnote-back-ref">↩</a></p></li> <li id="fn:5"><p>Moler, Cleve (19 June 2012). "Symplectic Spacewar". MATLAB Central - Cleve's Corner. MATLAB. Retrieved 2014-07-21. <a href="http://blogs.mathworks.com/cleve/2012/06/19/symplectic-spacewar/#comment-13" target="_blank">http://blogs.mathworks.com/cleve/2012/06/19/symplectic-spacewar/#comment-13</a> <a href="#fnref:5" class="footnote-back-ref">↩</a></p></li> <li id="fn:6"><p>Sommefeldt, Rys (2006-12-19). "Origin of Quake3's Fast InvSqrt() - Part Two". Beyond3D. Retrieved 2008-04-19. <a href="http://www.beyond3d.com/content/articles/15/" target="_blank">http://www.beyond3d.com/content/articles/15/</a> <a href="#fnref:6" class="footnote-back-ref">↩</a></p></li> <li id="fn:7"><p>Sommefeldt, Rys (2006-11-29). "Origin of Quake3's Fast InvSqrt()". Beyond3D. Retrieved 2009-02-12. <a href="http://www.beyond3d.com/content/articles/8/" target="_blank">http://www.beyond3d.com/content/articles/8/</a> <a href="#fnref:7" class="footnote-back-ref">↩</a></p></li> <li id="fn:8"><p>Blinn 1997, pp. 80–84. - Blinn, Jim (July 1997). "Floating Point Tricks". IEEE Computer Graphics & Applications. 17 (4): 80. doi:10.1109/38.595279. <a href="https://doi.org/10.1109%2F38.595279" target="_blank">https://doi.org/10.1109%2F38.595279</a> <a href="#fnref:8" class="footnote-back-ref">↩</a></p></li> <li id="fn:9"><p>Peelar, Shane (1 June 2021). "Fast reciprocal square root... in 1997?!". <a href="https://inbetweennames.net/blog/2021-05-06-i76rsqrt/" target="_blank">https://inbetweennames.net/blog/2021-05-06-i76rsqrt/</a> <a href="#fnref:9" class="footnote-back-ref">↩</a></p></li> <li id="fn:10"><p>Sommefeldt, Rys (2006-12-19). "Origin of Quake3's Fast InvSqrt() - Part Two". Beyond3D. Retrieved 2008-04-19. <a href="http://www.beyond3d.com/content/articles/15/" target="_blank">http://www.beyond3d.com/content/articles/15/</a> <a href="#fnref:10" class="footnote-back-ref">↩</a></p></li> <li id="fn:11"><p>"Discussion on CSDN". Archived from the original on 2015-07-02. <a href="https://web.archive.org/web/20150702180504/http://bbs.csdn.net/topics/41888" target="_blank">https://web.archive.org/web/20150702180504/http://bbs.csdn.net/topics/41888</a> <a href="#fnref:11" class="footnote-back-ref">↩</a></p></li> <li id="fn:12"><p>Lomont 2003, p. 1-2. - Lomont, Chris (February 2003). "Fast Inverse Square Root" (PDF). Retrieved 2009-02-13. <a href="http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf" target="_blank">http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf</a> <a href="#fnref:12" class="footnote-back-ref">↩</a></p></li> <li id="fn:13"><p>Sommefeldt, Rys (2006-11-29). "Origin of Quake3's Fast InvSqrt()". Beyond3D. Retrieved 2009-02-12. <a href="http://www.beyond3d.com/content/articles/8/" target="_blank">http://www.beyond3d.com/content/articles/8/</a> <a href="#fnref:13" class="footnote-back-ref">↩</a></p></li> <li id="fn:14"><p>Sommefeldt, Rys (2006-12-19). "Origin of Quake3's Fast InvSqrt() - Part Two". Beyond3D. Retrieved 2008-04-19. <a href="http://www.beyond3d.com/content/articles/15/" target="_blank">http://www.beyond3d.com/content/articles/15/</a> <a href="#fnref:14" class="footnote-back-ref">↩</a></p></li> <li id="fn:15"><p>Zafar, Saad; Adapa, Raviteja (January 2014). "Hardware architecture design and mapping of 'Fast Inverse Square Root' algorithm". 2014 International Conference on Advances in Electrical Engineering (ICAEE). pp. 1–4. doi:10.1109/ICAEE.2014.6838433. ISBN 978-1-4799-3543-7. S2CID 2005623. <a href="978-1-4799-3543-7" target="_blank">978-1-4799-3543-7</a> <a href="#fnref:15" class="footnote-back-ref">↩</a></p></li> <li id="fn:16"><p>Middendorf 2007, pp. 155–164. - Middendorf, Lars; Mühlbauer, Felix; Umlauf, George; Bodba, Christophe (June 1, 2007). "Embedded Vertex Shader in FPGA" (PDF). In Rettberg, Achin (ed.). Embedded System Design: Topics, Techniques and Trends. IFIP TC10 Working Conference:International Embedded Systems Symposium (IESS). et al. Irvine, California: Springer. doi:10.1007/978-0-387-72258-0_14. ISBN 978-0-387-72257-3. Archived (PDF) from the original on 2019-05-01. <a href="https://link.springer.com/content/pdf/10.1007/978-0-387-72258-0_14.pdf" target="_blank">https://link.springer.com/content/pdf/10.1007/978-0-387-72258-0_14.pdf</a> <a href="#fnref:16" class="footnote-back-ref">↩</a></p></li> <li id="fn:17"><p>Blinn 2003, p. 130. - Blinn, Jim (2003). Jim Blinn's Corner: Notation, notation notation. Morgan Kaufmann. ISBN 1-55860-860-5. <a href="#fnref:17" class="footnote-back-ref">↩</a></p></li> <li id="fn:18"><p>Sommefeldt, Rys (2006-11-29). "Origin of Quake3's Fast InvSqrt()". Beyond3D. Retrieved 2009-02-12. <a href="http://www.beyond3d.com/content/articles/8/" target="_blank">http://www.beyond3d.com/content/articles/8/</a> <a href="#fnref:18" class="footnote-back-ref">↩</a></p></li> <li id="fn:19"><p>"quake3-1.32b/code/game/q_math.c". Quake III Arena. id Software. Archived from the original on 2017-07-29. Retrieved 2017-01-21.{{cite web}}: CS1 maint: bot: original URL status unknown (link) <a href="https://web.archive.org/web/20170729072505/https://github.com/id-Software/Quake-III-Arena/blob/master/code/game/q_math.c#L552" target="_blank">https://web.archive.org/web/20170729072505/https://github.com/id-Software/Quake-III-Arena/blob/master/code/game/q_math.c#L552</a> <a href="#fnref:19" class="footnote-back-ref">↩</a></p></li> <li id="fn:20"><p>Eberly 2001, p. 504. - Eberly, David (2001). 3D Game Engine Design. Morgan Kaufmann. ISBN 978-1-55860-593-0. <a href="https://archive.org/details/isbn_9781558605930" target="_blank">https://archive.org/details/isbn_9781558605930</a> <a href="#fnref:20" class="footnote-back-ref">↩</a></p></li> <li id="fn:21"><p>Lomont 2003, p. 1. - Lomont, Chris (February 2003). "Fast Inverse Square Root" (PDF). Retrieved 2009-02-13. <a href="http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf" target="_blank">http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf</a> <a href="#fnref:21" class="footnote-back-ref">↩</a></p></li> <li id="fn:22"><p>Lomont 2003. - Lomont, Chris (February 2003). "Fast Inverse Square Root" (PDF). Retrieved 2009-02-13. <a href="http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf" target="_blank">http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf</a> <a href="#fnref:22" class="footnote-back-ref">↩</a></p></li> <li id="fn:23"><p>Use of the type long reduces the portability of this code on modern systems. For the code to execute properly, sizeof(long) must be 4 bytes, otherwise negative outputs may result. Under many modern 64-bit systems, sizeof(long) is 8 bytes. The more portable replacement is int32_t. <a href="/wiki/Sizeof" target="_blank">/wiki/Sizeof</a> <a href="#fnref:23" class="footnote-back-ref">↩</a></p></li> <li id="fn:24"><p>Sommefeldt, Rys (2006-11-29). "Origin of Quake3's Fast InvSqrt()". Beyond3D. Retrieved 2009-02-12. <a href="http://www.beyond3d.com/content/articles/8/" target="_blank">http://www.beyond3d.com/content/articles/8/</a> <a href="#fnref:24" class="footnote-back-ref">↩</a></p></li> <li id="fn:25"><p>Lomont 2003, p. 3. - Lomont, Chris (February 2003). "Fast Inverse Square Root" (PDF). Retrieved 2009-02-13. <a href="http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf" target="_blank">http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf</a> <a href="#fnref:25" class="footnote-back-ref">↩</a></p></li> <li id="fn:26"><p>McEniry 2007, p. 2, 16. - McEniry, Charles (August 2007). "The Mathematics Behind the Fast Inverse Square Root Function Code" (PDF). Archived from the original (PDF) on 2015-05-11. <a href="https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf" target="_blank">https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf</a> <a href="#fnref:26" class="footnote-back-ref">↩</a></p></li> <li id="fn:27"><p>Eberly 2001, p. 2. - Eberly, David (2001). 3D Game Engine Design. Morgan Kaufmann. ISBN 978-1-55860-593-0. <a href="https://archive.org/details/isbn_9781558605930" target="_blank">https://archive.org/details/isbn_9781558605930</a> <a href="#fnref:27" class="footnote-back-ref">↩</a></p></li> <li id="fn:28"><p>Ruskin, Elan (2009-10-16). "Timing square root". Some Assembly Required. Archived from the original on 2021-02-08. Retrieved 2015-05-07. <a href="https://web.archive.org/web/20210208132927/http://assemblyrequired.crashworks.org/timing-square-root/" target="_blank">https://web.archive.org/web/20210208132927/http://assemblyrequired.crashworks.org/timing-square-root/</a> <a href="#fnref:28" class="footnote-back-ref">↩</a></p></li> <li id="fn:29"><p>Fog, Agner. "Lists of instruction latencies, throughputs and micro-operation breakdowns for Intel, AMD and VIA CPUs" (PDF). Retrieved 2017-09-08. <a href="http://www.agner.org/optimize/instruction_tables.pdf" target="_blank">http://www.agner.org/optimize/instruction_tables.pdf</a> <a href="#fnref:29" class="footnote-back-ref">↩</a></p></li> <li id="fn:30"><p>"§ 6.5/7" (PDF), ISO/IEC 9899:2018, 2018, p. 55, archived from the original (PDF) on 2018-12-30, An object shall have its stored value accessed only by an lvalue expression that has one of the following types: [...] <a href="https://web.archive.org/web/20181230041359/http://www.open-std.org/jtc1/sc22/wg14/www/abq/c17_updated_proposed_fdis.pdf" target="_blank">https://web.archive.org/web/20181230041359/http://www.open-std.org/jtc1/sc22/wg14/www/abq/c17_updated_proposed_fdis.pdf</a> <a href="#fnref:30" class="footnote-back-ref">↩</a></p></li> <li id="fn:31"><p>"§ 6.5.2.3/3, footnote 97" (PDF), ISO/IEC 9899:2018, 2018, p. 59, archived from the original (PDF) on 2018-12-30, If the member used to read the contents of a union object is not the same as the member last used to store a value in the object, the appropriate part of the object representation of the value is reinterpreted as an object representation in the new type as described in 6.2.6 (a process sometimes called "type punning"). <a href="https://web.archive.org/web/20181230041359/http://www.open-std.org/jtc1/sc22/wg14/www/abq/c17_updated_proposed_fdis.pdf" target="_blank">https://web.archive.org/web/20181230041359/http://www.open-std.org/jtc1/sc22/wg14/www/abq/c17_updated_proposed_fdis.pdf</a> <a href="#fnref:31" class="footnote-back-ref">↩</a></p></li> <li id="fn:32"><p>"Bit-casting object representations". www.open-std.org. Archived from the original on 2018-08-18. Retrieved 2018-11-10. <a href="https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2018/p0476r2.html" target="_blank">https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2018/p0476r2.html</a> <a href="#fnref:32" class="footnote-back-ref">↩</a></p></li> <li id="fn:33"><p>Goldberg 1991, p. 7. - Goldberg, David (1991). "What every computer scientist should know about floating-point arithmetic". ACM Computing Surveys. 23 (1): 5–48. doi:10.1145/103162.103163. S2CID 222008826. <a href="https://doi.org/10.1145%2F103162.103163" target="_blank">https://doi.org/10.1145%2F103162.103163</a> <a href="#fnref:33" class="footnote-back-ref">↩</a></p></li> <li id="fn:34"><p>Goldberg 1991, pp. 15–20. - Goldberg, David (1991). "What every computer scientist should know about floating-point arithmetic". ACM Computing Surveys. 23 (1): 5–48. doi:10.1145/103162.103163. S2CID 222008826. <a href="https://doi.org/10.1145%2F103162.103163" target="_blank">https://doi.org/10.1145%2F103162.103163</a> <a href="#fnref:34" class="footnote-back-ref">↩</a></p></li> <li id="fn:35"><p>E x {\displaystyle E_{x}} should be in the range [ 1 , 254 ] {\displaystyle [1,254]} for x {\displaystyle x} to be representable as a normal number. <a href="/wiki/Normal_number_(computing)" target="_blank">/wiki/Normal_number_(computing)</a> <a href="#fnref:35" class="footnote-back-ref">↩</a></p></li> <li id="fn:36"><p>The only real numbers that can be represented exactly as floating point are those for which M x {\textstyle M_{x}} is an integer. Other numbers can only be represented approximately by rounding them to the nearest exactly representable number. <a href="#fnref:36" class="footnote-back-ref">↩</a></p></li> <li id="fn:37"><p>Goldberg 1991, p. 16. - Goldberg, David (1991). "What every computer scientist should know about floating-point arithmetic". ACM Computing Surveys. 23 (1): 5–48. doi:10.1145/103162.103163. S2CID 222008826. <a href="https://doi.org/10.1145%2F103162.103163" target="_blank">https://doi.org/10.1145%2F103162.103163</a> <a href="#fnref:37" class="footnote-back-ref">↩</a></p></li> <li id="fn:38"><p>McEniry 2007, p. 3. - McEniry, Charles (August 2007). "The Mathematics Behind the Fast Inverse Square Root Function Code" (PDF). Archived from the original (PDF) on 2015-05-11. <a href="https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf" target="_blank">https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf</a> <a href="#fnref:38" class="footnote-back-ref">↩</a></p></li> <li id="fn:39"><p>Since x {\displaystyle x} is positive, S x = 0 {\displaystyle S_{x}=0} . <a href="#fnref:39" class="footnote-back-ref">↩</a></p></li> <li id="fn:40"><p>Hennessey & Patterson 1998, p. 305. - Hennessey, John; Patterson, David A. (1998). Computer Organization and Design (2nd ed.). San Francisco, CA: Morgan Kaufmann Publishers. ISBN 978-1-55860-491-9. <a href="https://archive.org/details/computerorganiz000henn" target="_blank">https://archive.org/details/computerorganiz000henn</a> <a href="#fnref:40" class="footnote-back-ref">↩</a></p></li> <li id="fn:41"><p>Hardy 1908, p. 323. - Hardy, Godfrey (1908). A Course of Pure Mathematics. Cambridge, UK: Cambridge University Press. <a href="http://www.gutenberg.org/files/38769" target="_blank">http://www.gutenberg.org/files/38769</a> <a href="#fnref:41" class="footnote-back-ref">↩</a></p></li> <li id="fn:42"><p>McEniry 2007, p. 6. - McEniry, Charles (August 2007). "The Mathematics Behind the Fast Inverse Square Root Function Code" (PDF). Archived from the original (PDF) on 2015-05-11. <a href="https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf" target="_blank">https://web.archive.org/web/20150511044204/http://www.daxia.com/bibis/upload/406Fast_Inverse_Square_Root.pdf</a> <a href="#fnref:42" class="footnote-back-ref">↩</a></p></li> <li id="fn:43"><p>Eberly 2001, p. 2. - Eberly, David (2001). 3D Game Engine Design. Morgan Kaufmann. ISBN 978-1-55860-593-0. <a href="https://archive.org/details/isbn_9781558605930" target="_blank">https://archive.org/details/isbn_9781558605930</a> <a href="#fnref:43" class="footnote-back-ref">↩</a></p></li> <li id="fn:44"><p>Lomont 2003, p. 10. - Lomont, Chris (February 2003). "Fast Inverse Square Root" (PDF). Retrieved 2009-02-13. <a href="http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf" target="_blank">http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf</a> <a href="#fnref:44" class="footnote-back-ref">↩</a></p></li> <li id="fn:45"><p>Lomont 2003, p. 10. - Lomont, Chris (February 2003). "Fast Inverse Square Root" (PDF). Retrieved 2009-02-13. <a href="http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf" target="_blank">http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf</a> <a href="#fnref:45" class="footnote-back-ref">↩</a></p></li> <li id="fn:46"><p>Lomont 2003, pp. 10–11. - Lomont, Chris (February 2003). "Fast Inverse Square Root" (PDF). Retrieved 2009-02-13. <a href="http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf" target="_blank">http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf</a> <a href="#fnref:46" class="footnote-back-ref">↩</a></p></li> <li id="fn:47"><p>Matthew Robertson (2012-04-24). "A Brief History of InvSqrt" (PDF). UNBSJ. <a href="https://mrober.io/papers/rsqrt.pdf" target="_blank">https://mrober.io/papers/rsqrt.pdf</a> <a href="#fnref:47" class="footnote-back-ref">↩</a></p></li> <li id="fn:48"><p>Kadlec, Jan (2010). "Řrřlog::Improving the fast inverse square root" (personal blog). Archived from the original on 2018-07-09. Retrieved 2020-12-14. <a href="https://web.archive.org/web/20180709021629/http://rrrola.wz.cz/inv_sqrt.html" target="_blank">https://web.archive.org/web/20180709021629/http://rrrola.wz.cz/inv_sqrt.html</a> <a href="#fnref:48" class="footnote-back-ref">↩</a></p></li> <li id="fn:49"><p>Moroz et al. 2018. - Moroz, Leonid V.; Walczyk, Cezary J.; Hrynchyshyn, Andriy; Holimath, Vijay; Cieslinski, Jan L. (January 2018). "Fast calculation of inverse square root with the use of magic constant analytical approach". Applied Mathematics and Computation. 316 (C). Elsevier Science Inc.: 245–255. arXiv:1603.04483. doi:10.1016/j.amc.2017.08.025. S2CID 7494112. <a href="https://arxiv.org/abs/1603.04483" target="_blank">https://arxiv.org/abs/1603.04483</a> <a href="#fnref:49" class="footnote-back-ref">↩</a></p></li> <li id="fn:50"><p>Muller, Jean-Michel (December 2020). "Elementary Functions and Approximate Computing". Proceedings of the IEEE. 108 (12): 2146. doi:10.1109/JPROC.2020.2991885. ISSN 0018-9219. S2CID 219047769. <a href="https://ieeexplore.ieee.org/document/9106347" target="_blank">https://ieeexplore.ieee.org/document/9106347</a> <a href="#fnref:50" class="footnote-back-ref">↩</a></p></li> <li id="fn:51"><p>Ruskin, Elan (2009-10-16). "Timing square root". Some Assembly Required. Archived from the original on 2021-02-08. Retrieved 2015-05-07. <a href="https://web.archive.org/web/20210208132927/http://assemblyrequired.crashworks.org/timing-square-root/" target="_blank">https://web.archive.org/web/20210208132927/http://assemblyrequired.crashworks.org/timing-square-root/</a> <a href="#fnref:51" class="footnote-back-ref">↩</a></p></li> </ol>