Everything I know about the fast inverse square root algorithm

If your computer was built after 1999, it probably supports the SSE instruction set. It contains the _mm_rsqrt_ps instruction, which is faster and will give you four reciprocal square roots at once: <a href="https:&#x2F;&#x2F;www.intel.com&#x2F;content&#x2F;www&#x2F;us&#x2F;en&#x2F;docs&#x2F;intrinsics-guide&#x2F;index.html#text=_mm_rsqrt_ps&amp;ig_expand=5642&amp;ssetechs=SSE" rel="nofollow">https:&#x2F;&#x2F;www.intel.com&#x2F;content&#x2F;www&#x2F;us&#x2F;en&#x2F;docs&#x2F;intrinsics-guid...</a><p>That being said, the techniques discussed here are not totally irrelevant (yet). There still exists some hardware with fast instructions for float&#x2F;int conversion, but lacking rsqrt, sqrt, pow, log instructions, which can all be approximated with this nice trick.

To nitpick the article a bit:<p>&gt; It&#x27;s important to note that this algorithm is very much of its time. Back when Quake 3 was released in 1999, computing an inverse square root was a slow, expensive process. The game had to compute hundreds or thousands of them per second in order to solve lighting equations, and other 3D vector calculations that rely on normalization. These days, on modern hardware, not only would a calculation like this not take place on the CPU, even if it did, it would be fast due to much more advanced dedicated floating point hardware.<p>Calculations like this definitely take place on the CPU all the time. It&#x27;s a common misconception that games and other FLOP-heavy apps want to offload all floating-point operations to the GPU. In fact, it really only makes sense to offload large <i>uniform</i> workloads to the GPU. If you&#x27;re doing one-off vector normalization--say, as part of the rotation matrix construction needed to make one object face another--then you&#x27;re going to want to stay on the CPU, because the CPU is faster at that. In fact, the CPU would remain faster at <i>single</i> floating point operations even if you didn&#x27;t take the GPU transfer time into account--the GPU typically runs at a slower clock rate and relies on parallelism to achieve its high FLOP count.

I wrote a implementation in MMIX:<p><pre><code> % Constants FISRCON GREG #5FE6EB50C7B537A9 THREHAF GREG #3FF8000000000000 % Save half of the original number OR $2,$0,0 INCH $2,#FFF0 % Bit level hacking SRU $1,$0,1 SUBU $0,FISRCON,$1 % First iteration FMUL $1,$2,$0 FMUL $1,$1,$0 FSUB $1,THREHAF,$1 FMUL $0,$0,$1 % Second iteration FMUL $1,$2,$0 FMUL $1,$1,$0 FSUB $1,THREHAF,$1 FMUL $0,$0,$1 </code></pre> This implementation makes an assumption that the original number is greater than 2^-1021.

If you&#x27;re interested, Wikipedia also has a decent discussion of this function and it&#x27;s history. <a href="https:&#x2F;&#x2F;en.wikipedia.org&#x2F;wiki&#x2F;Fast_inverse_square_root" rel="nofollow">https:&#x2F;&#x2F;en.wikipedia.org&#x2F;wiki&#x2F;Fast_inverse_square_root</a>

I collected a few of these here:<p><a href="https:&#x2F;&#x2F;github.com&#x2F;ncruces&#x2F;fastmath&#x2F;blob&#x2F;main&#x2F;fast.go">https:&#x2F;&#x2F;github.com&#x2F;ncruces&#x2F;fastmath&#x2F;blob&#x2F;main&#x2F;fast.go</a><p>Also see this StackOverflow:<p><a href="https:&#x2F;&#x2F;stackoverflow.com&#x2F;questions&#x2F;32042673&#x2F;optimized-low-accuracy-approximation-to-rootnx-n" rel="nofollow">https:&#x2F;&#x2F;stackoverflow.com&#x2F;questions&#x2F;32042673&#x2F;optimized-low-a...</a>

Time for nitpicks, sorry.<p>The formulas for the floats have typos. They should read (-1)^S, not -1^S (which always equals -1).<p>Interpreting the raw bit patterns isn&#x27;t a piecewise linear approximation of the logarithm. The lines between the data points on the blue graph don&#x27;t actually exist, it&#x27;s not possible for a bit to be half set to 1. It&#x27;s rather a discrete version of the logarithm: the only data points that exist - where the red and blue lines meet - are literally equal to the (scaled, shifted) logarithm.<p>Other than that, nice post!

This is a good post explaining a lot of interesting concepts, with a section that has surprisingly bad algebra.<p>&gt; The exact steps to go from the first form to this one are numerous to say the least, but I&#x27;ve included it in full for completeness.<p>The algebra included after that line has many unnecessary steps, as well as multiple cancelling sign errors. Most notably the second line to the third line does not distribute the negative sign correctly. Picking up after the second line, I would simply have:<p>y_n+1 = y_n + (1 - x * y_n^2) &#x2F; y_n^2 * (y_n^3 &#x2F; 2)<p>y_n+1 = y_n + (1 - x * y_n^2) * (y_n &#x2F; 2)<p>y_n+1 = y_n + y_n &#x2F; 2 - x * y_n^3&#x2F;2<p>y_n+1 = 3&#x2F;2 * y_n - x * y_n^3&#x2F;2<p>y_n+1 = y_n (3&#x2F;2 * y_n - x * y_n^2 &#x2F; 2)<p>y_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n)<p>This is fewer than 1&#x2F;3 of the included steps (between the second line and the end), and has the advantage of being correct during the intermediate steps. I don&#x27;t think any of my steps are anything other than self-evident to someone who understands algebra, but I&#x27;d be happy to entertain competing perspectives.

Here’s something new that isn’t mentioned and might be worth adding to the repo. (Edit I’m wrong, it is there, I just didn’t recognize it.) The magic number in the famous code snippet is not the optimal constant. You can do maybe 0.5% better relative error using a different constant. Maybe at the time it was infeasible to search for the absolutely optimal number, but now it’s relatively easy. I also went down this rabbit hole at some point, so I have a Jupyter notebook for finding the optimal magic number for this (1&#x2F;x^2) and also for (1&#x2F;x). Anyone wanna know what the optimal magic number is?

I thought the most interesting bit in this article was the link to &quot;How Java&#x27;s Floating-Point Hurts Everyone Everywhere&quot;:<p><a href="https:&#x2F;&#x2F;people.eecs.berkeley.edu&#x2F;~wkahan&#x2F;JAVAhurt.pdf" rel="nofollow">https:&#x2F;&#x2F;people.eecs.berkeley.edu&#x2F;~wkahan&#x2F;JAVAhurt.pdf</a><p>Written by William Kahan, also known as the Old Man of Floating-Point:<p><a href="https:&#x2F;&#x2F;news.ycombinator.com&#x2F;item?id=29042853">https:&#x2F;&#x2F;news.ycombinator.com&#x2F;item?id=29042853</a> - An Interview with the Old Man of Floating-Point (1998)

I really liked this video about it when I saw it a while ago: <a href="https:&#x2F;&#x2F;www.youtube.com&#x2F;watch?v=p8u_k2LIZyo" rel="nofollow">https:&#x2F;&#x2F;www.youtube.com&#x2F;watch?v=p8u_k2LIZyo</a>

I was building 3D engines a few years before the days of Quake, and having recently watched videos about optimizing the trig code in Super Mario 64, and these other faster algorithms, it makes me wonder how much more can be squeezed out of old hardware.<p>I was optimizing my assembler to the nth degree, but optimizing the algorithm is always going to be the real winner.

Something interesting is I&#x27;ve seen this mythologized as evidenced of the &quot;cracked engineer&quot; theory, wherein random engineers will accidentally stumble upon incredibly complex discoveries in the course of their day to day work and expect no fanfare for this, besting the scientists and researchers who take the traditional route but aren&#x27;t spurred by necessity. People, like xkcd, purported this was either Carmack or some random employee who figured it out in a few hours when stuck on a problem then waltzed onto the next one. In truth, as noted in the document, it has a history that goes back decades in academia, this was simply the most notable time it was implemented.<p>I think it speaks to an issue I see in the software engineering world where people assume that collaboration is for low-IQ people and all great innovation comes from some super-genius working on their own for long enough, and isn&#x27;t required to share how they arrived at their findings. I&#x27;m sure the mythology of this algorithm is propelled somewhat by the enigmatic character of writing &quot;what the fuck&quot; next to adding the constant, implying a mystical element to its utility that was arrived at without needing to clarify it in some long boring research paper.

Note that you’re actually using 3 distinct magic numbers - not only the single 0x5f3759df, but also 0.5 and 1.5;<p>So for further accuracy we can instead have:<p>conv.i = 0x5F1FFFF9 - ( conv.i &gt;&gt; 1 ); conv.f <i>= 0.703952253f </i> ( 2.38924456f - x * conv.f * conv.f ); return conv.f;<p>[from Wikipedia]

Thanks for stealing hours of my productive time going down this rabbit hole.

Noob here. Please ELI 5 what this algo used in quake 3 ?.<p>&gt; Solve lighting equations<p>Is it the light or laser effect when the gun is fired ?

I tried this in a raytracter today while learning rust and it made the scene render all white. Lol

Did John not write the code? Or did he copy paste it from somewhere?

What’s up with the diagonal lines in the graph backgrounds?

He posts c code and says that it computes the inverse square root, but the code he posted is undefined c, so while it could compute the inverse square root of it&#x27;s input it could also do a long list of other stuff.<p>One of the many reasons you should stay away from real world languages when talking about algorithms unless you are an expert in the language.