В 2010 году отечественная индустрия игр была на взлёте, появлялись всё больше и больше новых и интересных игр! Яна и Алексей тоже решили сделать свою крутую 3д игру “Дрожь Земли”, они достаточно быстро придумали, что там нужно будет бегать и стрелять в монстров, собирать оружие и классно прыгать. Яна нарисовала графику, а Алексей занялся написанием движка. Через какое-то время всё было закончено, но вот незадача среди ребят не оказалось ни одного любителя дискры и матанализа, а пары АКиОС тогда ещё нарушали правила ВШЭ(на них никто не ходил) и игра оказалась очень медленной и зависающей. Спустя 13 лет вам в руки попал тот самый код и вы решили, что сможете исправить и разобраться в коде. Покапавшись, вы узнали, что Алексей забыл математику, но пытался как-то реализовать функцию, которая будет уметь быстро считать обратный кубический корень от числа с некой апроксимацией. Ваша задача переписать эту функцию float fast_reverse_cube(float number) и применить её к массиву чисел, который подаётся на стандартный вход. Стандартную библиотеку math.h использовать запрещено.

Вам так же может быть полезно почитать про то, как это сделали ребята в другой игре, а так же более подробную статью про ту же игру.

На вход подается последовательность вещественных чисел. Для каждого выведите ответ

Для подсчёта обратного кубического корня циклы использовать запрещено

Examples

Input

1 2 3 4
Output

0.999990 0.793701 0.693354 0.62996
 

Submit a solution

Language:	gcc - GNU C 11.3.0

File	
Send!	

Sanity Phailed.me
Just another WordPress site. Learn Shit The Insane Way! - Math Problem of the Day.
0x5f400000: Understanding Fast Inverse Sqrt the Easy(ish) Way!
In this post, we attempt to unravel the mysteries of the “magic” constant found in the fast inverse sqrt method in an intuitive fashion.

28 October, 2014 - Algorithm, Article, C, Numerical Analysis - Tags : floating point, numerical analysis 2 Comments
  
Table Of Content
TL;DR – Start Here!

1. Introduction

WTF?

Obsession

2. Switching Strategy

Converting between a Float and its Bits in memory

Looking at b2f(f2b(x)*b)

3. Homing In on 0x5F400000

Adding Alpha

0x3F800000 – Generalized fast power method

4. Epilogue – The Hard Parts

Closing In On 0x5f3759df

Refinement – 0x5f376000

Newton Iteration – 0x5f375a1a

Closing Thoughts

Requested by averageisnothalfway: Differences between a1a and df

TL;DR – Start Here!
1. Let \g{f2b} be the float-to-bit function (that takes in a floating point number and outputs a 32-bit long representing its IEEE 754 encoding used on your everyday computer) and \g{b2f} be the bit-to-float function, then the fast inverse square root method can be rewritten as

    \[\sqrt{\tilde x^{-1}} \approx \g{b2f}(\g{0x5f3759df} - \g{f2b}(x)/2)\]

2. It turns out that \g{b2f}(\alpha + \g{f2b}(x)\times \beta) = C(x,\alpha, \beta) \times x^{\beta} for some non-linear function C(x,\alpha,\beta) \approx C(\alpha,\beta) that doesn’t vary a lot with x, so it becomes an optimization game to find \alpha such that C(\alpha, \beta) = 1.
3. We know that x^{\beta} = 1 when x = 1 no matter what \beta is, so if we want to find \alpha such that \g{b2f}(\alpha + \beta\g{f2b}(x)) \approx x^{\beta}, all we need to do is set x = 1 because then

    \[\g{b2f}(\alpha + \beta\g{f2b}(1)) = 1 = \g{b2f}(\g{f2b}(1))\]

so \alpha = \underbrace{\g{f2b}(1)}_{\g{0x3f800000}}(1-\beta), which let’s us make arbitrary fast \beta^{th} power methods on demand.
4. ???
5. Profit!

1. Introduction
The first time I saw the magical fast inverse square root method, I was a freshman in college and the code pretty much killed my brain. I’d just learned about algorithms and data structures and I thought that binary search trees were like the coolest things ever. Turns out I didn’t know them well enough to tackle this problem back then.

WTF?
If you’re not familiar with the fast inverse square root method, here’s what it looks like in C (source: wikipedia)

view sourceprint?
01.
float Q_rsqrt( float number )
02.
{
03.
long i;
04.
float x2, y;
05.
const float threehalfs = 1.5F;
06.
 
07.
x2 = number * 0.5F;
08.
y  = number;
09.
i  = * ( long * ) &y;                        // evil floating point bit level hacking
10.
i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
11.
y  = * ( float * ) &i;
12.
y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
13.
//      y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed
14.
 
15.
return y;
16.
}
[Note: The use of pointer aliasing here is illegal, as Maristic pointed out in the followup discussion on reddit. Use unions instead.]
Okay, so first you declare i, x2, y, and a constant threehalfs. Next, you set x2 to be \frac{x}{2} and set y = x. So far so good!

Okay, next, i gets the dereferenced value of the float* pointer to y casted into a long* pointer instead?? Okay…

Next, i = 0x5f3759df – (i/2)??? Huh?

y = *(float*)&i????

y = y * ( threehalfs – ( x2 * y * y ) )??????

What. The. Fucking. Fuck.

Obsession
Before we go on, I should probably warn you that this code was a constant obsession of mine for the better part of my adult life. Because of this, a lot of what I’m about to say may sound preachy, cheesy, or even outright annoying. I’ll try to tone it down, but truth be told, I am proud of this work; I’ve come to see the fast inverse square root method as my baby and because of that, I’ve become quite possessive of it. As much as this is a post on the possible derivation of 0x5f3759df, this is also a memoir of my college life (minus the fun parts, like parties with buzz and services as Zoom Escorts).

From the very beginning, I tried to wrap my head around this strange game of algebraic reasoning, and for the next 4 years, my academic life largely revolved around this mystery: I spent my winter break freshman year learning C and figured out that the “bit level” hack gave you the array of bits that your computer encoded floating point numbers as; sophomore year, I took a first course on numerical methods and understood how iterating the sequence

    \[y_{k+1} = y_k \pa{\frac{3}{2} - \frac{x}{2}y_k^2}\]

ad infinitum will eventually converge to y_{\infinity} = \frac{1}{\sqrt{x}}; then, during the summer before my junior year, I worked through a lot of tedious algebra dealing with the quirks of IEEE754 and found another bit-level hack to quickly compute the square root a la the spirit of the fast inverse square root method, but this time exploiting the distribution of floats and a convenient analogy to second order taylor expansions.

Indeed, my obsession with this method has given me an appetite for numerical analysis, programming languages, software verification, compilers, and even combinatorics and dynamical systems.

Without it, I would have never explored most of the areas of Noise Control science that I currently consider my home.

But through it all, there was always a constant question at the back of my mind that I’ve failed to answer time and time again: where did the 0x5f3759df come from?

I never got to a point where I can claim that I have a satisfactory answer to that question, but I think I did come to discover something along the road: a different magical constant 0x5f400000 that I believe lay at the very beginning of the history of this method, and it is the joy of discovering this constant that I would like to gift to my readers.

2. Switching Strategy
Towards the end of my college years, I realized that there was a big problem in the way I approached the problem of looking at 0x5f3759df. Most of the traditional research into the subject has been through the perspective that whoever wrote the algorithm went after the most optimal method possible. In the process of hunting for the jackpot, a lot of great intuition was left behind because those ideas weren’t efficient or effective enough. This is what led me down the path of looking for inefficient but elegant methods that could inspire the fast inverse square root method.

But first, let’s ramp up on some preliminaries. A lot of work and paper on this matter jump immediately to the representation of floating point numbers; I will completely forgo that Uber Jobs Near Me instead of relying on other things to do the job.

Converting between a Float and its Bits in memory
Let’s first consider two functions: \g{f2b}(n), which stands for float-to-bits will return the IEEE 754 bit array representation of the floating number n, effectively performing *(long*)&y; and \g{b2f}(b), which stands for bits-to-float and takes a bit array and turns it back into a float. Note that \g{f2b} = \g{b2f}^{-1} and vice versa so that \g{b2f}(\g{f2b}(n)) = n.

It turns out that the magical fast inverse square root method is equivalent to

    \[\g{qsqrt}(n) = \g{b2f}\pa{\g{0x5f3759df} - \frac{\g{f2b}(n)}{2}}\]

Looking at b2f(f2b(x)*b)
With the conversion functions out of the way, I next turned my attention to everything but the magic constant. Specifically, what are the effects of multiplying inside the \g{b2f}. I fired up my favorite plotting library and looked at

    \[\g{b2f}(\g{f2b}(x) \times 0.5).\]

After a little bit of wrestling around with the plot, I realized that when plotting \g{b2f}(\g{f2b}(x) \times 0.5) in log-log style, it has the same shape as x^{0.5} = \sqrt{x}. This means that

    \[\g{b2f}(\g{f2b}(x) \times 0.5) \approx C \times \sqrt{x}\]

for some constant C.

I fitted a “regression” line (this isn’t a real regression however, don’t be mislead) through the data onto \sqrt{x} and found that the relative error is quite small as can be seen in the figure below.

Values
Constant(8.131516293641283e-20) * x^0.5
b2f(f2b(x)*0.5)
0.1
1
10
1e-20
2e-20
4e-20
1e-19
2e-19
4e-19
1e-18
Highcharts.com
Log-Log plot of C*x0.5 versus b2f(f2b(x)*0.5));
Values
Constant(8.131516293641283e-20) * x^0.5
b2f(f2b(x)*0.5)
x^0.5
0.1
1
10
1e-74
1e-37
1
1e+37
1e+74
Highcharts.com
Log-Log plot of x0.5
Values
Plot of the relative error
relative error
0.1
1
10
0
0.1
0.2
0.3
Highcharts.com
    \[\g{change }\beta\]

0.5
In fact, if you tick through the slider at the bottom of the figure, you’ll come to find that for any arbitrary \beta, it seems to be the case that

    \[\g{b2f}(\g{f2b}(x) \times \beta) \approx C_{\beta} \times x^{\beta}\]

for some constant C_{\beta}.

3. Homing In on 0x5F400000
Adding Alpha
Let’s now consider the following function:

    \[\gamma(x, \alpha,\beta) = \g{b2f}\pa{\g{f2b}(x)\times \beta + \alpha}\]

Notice here that we now have added an additive term \alpha to the equation. Again, we can rewrite the fast inverse square root method as

    \[\sqrt{\tilde x^{-1}}= \gamma\pa{x, \g{0x5f3759df}, -\frac{1}{2}}\]

A natural question one might consider is the effect of varying \alpha on the behavior of \gamma. It turns out that on the log-log graph of \gamma(x, \alpha, \beta), varying \alpha does not change the slope of the line, but only its y-intercept. In other words, varying \alpha will only affect magnitude of the constant C_{\beta}, as we can see below.

Values
b2f(0 - f2b(x)/2)
1/sqrt(x)
0.1
1
10
1e-38
1e-19
1
10000000000000000000
1e+38
Highcharts.com
Log-Log plot of b2f(0-f2b(x)/2)
    \[\g{change }\alpha\]

0
This then suggests that \beta determines the degree of \gamma while \alpha determines its magnitude. Symbolically:

    \[\gamma(x, \alpha, \beta) \approx C_\beta(\alpha) x^{\beta}\]

where C_{\beta}(\alpha) is some constant.

0x3F800000 – Generalized fast power method
We are now left at an interesting point in the puzzle. We have on one hand the knowledge that

    \[\gamma(x, \alpha, \beta) \approx C_\beta(\alpha) x^{\beta}\]

for some function C_\beta(\alpha) and we want to determine \alpha such that

    \[\gamma(x,\alpha,\beta) \approx C_\beta(\alpha) x^{\beta} = x^{\beta}\]

Canceling out the x^{\beta} on either side, we are left with the functional equation

    \[C_\beta(\alpha) = 1\]

So right now we need to find \alpha such that C_\beta(\alpha) = 1. But for x = 1, x^{\beta} = 1 for any \beta, so the naive thing to do here is to reduce our approximation to

    \begin{align*} \gamma(1, \alpha, \beta) &\approx C_\beta(\alpha) = 1 \\ \g{b2f}(\g{f2b}(1)\times \beta + \alpha) &\approx 1 \\ & = \g{b2f}(\g{f2b}(1)) \\ \g{f2b}(1)\times \beta + \alpha &\approx \g{f2b}(1) \\ \alpha &\approx \g{f2b}(1) \times (1-\beta) \\ &= \g{0x3f800000} \times (1-\beta) \end{align*}

Let’s take a moment and appreciate this statement:

    \[\alpha \approx \g{0x3f800000} \times (1-\beta)\]

This says that for any exponent \beta (positive or negative, integer or fractional), we can approximate x^{\beta} with

    \[\g{b2f}\pa{\g{f2b}(x)\times \beta + \g{0x3f800000} \times (1-\beta)}\]

That’s mind blowing because we can immediately generate “fast \beta^{th} power” methods using nothing but \beta!

view sourceprint?
1.
float fast_pow(float x, float n) {
2.
long bit = (*((long*)&x))*n + 0x3f800000*(1-n);
3.
return *((float*)&bit);
4.
}
Going back to inverse square roots, we have \alpha = \g{0x3f800000} \times (1- (-0.5)) = \g{0x5f400000}, which is the namesake of this article.

4. Epilogue – The Hard Parts
Closing In On 0x5f3759df
After discovering the existence of 0x5f400000, things became much clearer. However, it was still quite off from SGI’s 0x5f3759df or Lomont’s 0x5f375a86, and so the question of how this constant came to be still nagged at the back of my head. Nevertheless, I am certain that 0x5f400000 appeared at least in some shape or form in the earlier versions of this numerical recipe.

Many others have previously discussed the possible origins of 0x5f3759df, and in the rest of this post, I will offer my own views on the matter.

For the most part, I believe that the discovery of 0x5f3759df is an iterative process: it starts with 0x5f40, then it gets refined to 0x5f375X, and finally someone must have decided on 0x5f3759df. As engineers are practitioners by nature, I don’t believe that anyone that worked on 0x5f3759df cared for its optimality; therefore, I believe that intuition and visualization play the most important part in its discovery.

In this section, I will first outline one path to refine 0x5f40 to 0x5f375X through intuition about the maximum error bounds of the recipe, followed by a discussion on the final discovery as a product of analyzing one step of newton’s iteration.

Refinement – 0x5f376000
If we look at the graph of the relative error of \gamma\pa{x, \g{0x5f400000}, -\frac{1}{2}} on a logarithmic x-axis below, you will notice a couple of peculiar things.

Values
Plot of the relative error of 0x5f400000 and friends
relative error of 0x5f400000
relative error of 0x5f390000
relative error of 0x5f3759DF
0.1
1
10
-0.1
-0.05
0
0.05
0.1
Highcharts.com
relative error of 0x5f400000x = 0.5509525579382814, y = -0.0755725444596338
First and foremost, changing the value of the magic constant by small amounts seems to lift the relative error up. Next, the error is periodic with respect to \log_2(x) with a period of 2 (so that \g{re}(x) = \g{re}(2^{\log_2(x)+2}) = \g{re}(4x)). Finally, the position of the maximum and the minimum of the relative error is largely unchanged when we vary the magic constant.

We can use these observations to our advantage. Notice that the error for \g{magic} = \g{0x5f400000} is completely negative. However, if we were to use a smaller magic constant, we can shift part of the error up into the positive zone, which reduces the absolute error as can be seen below.

Values
Plot of the absolute relative error of 0x5f400000 and friends
abs relative error of 0x5f400000
abs relative error of 0x5f390000
abs relative error of 0x5f3759DF
0.1
1
10
0
0.025
0.05
0.075
0.1
Highcharts.com
We are trying to reduce the error bounds as much as possible (under the infinity norm optimization as we would say). In order to do so, we need to ensure that the maximum point plus the minimum point in the non-absolute relative error plot is as close to 0 as possible.

Here’s what we’re going to do. First, we’re going to zoom in onto the section of the graph within the interval x \in [1,4] because it is periodic. Next, we will look at the maximum point and minimum point and look at their sum.

Values
Plot of the relative error of 0x5f400000
max 1.00000max 1.00000
min 2.67586min 2.67586
sum = -0.08865725622311937sum = -0.08865725622311937
relative error of 0x5f400000
Series 2
1
-0.2
-0.1
0
0.1
0.2
Highcharts.com
    \[\g{change magic constant }m\]

5f400000
You will see that the optimal constant here lies somewhere close to 0x5f376000, which is about 10^3 away from the SGI constant.

Newton Iteration – 0x5f375a1a
In our final trick tonight, we will refine this constant down further to within 60 away from SGI’s 0x5f3759df.

One of the important things to notice from the quake code is that the second newton iteration was commented away. This suggests that whoever worked on the latest version of the fast inverse square root method must have optimized the constant with the effects of newton iteration applied. Let’s look at how well the constants do with the help of an iteration.

Values
Plot of the relative error of 0x5f400000 and friends after 1 newton iteration
relative error of 0x5f400000
relative error of 0x5f3759DF
0.1
1
10
0
0.005
0.01
0.015
Highcharts.com
Wow, using 0x5f3759df as the magic constant improved the accuracy nearly tenfolds! This likely suggests that whoever worked on quake’s version of fast inverse square root method sought after 4 digits of accuracy. In order to achieve this with the vanilla constant 0x5f400000, they would need to run two newton iterations.

Let’s dig in a little bit and try to find the magic constant that minimizes the maximum error after one newton iteration.

Values
Plot of the relative error of 0x5f376000 after 1 newton step.
min y = 1.403799362309499e-7min y = 1.403799362309499e-7
max y = 0.0017647334370672896max y = 0.0017647334370672896
relative error of 0x5f376000
1
0
0.0005
0.001
0.0015
0.002
Highcharts.com
    \[\g{change magic constant }m\]

5f376000
Incidentally, this occurs when the right most two peeks in the above graph collide. We find that the magic constant 0x5f375a1a works best here, which is only 59 away from 0x5f3759df.

We can look at the mean error as well (the one-norm):

Values
Mean error of 0x5f376000 after 1 newton step: 0.0009416359957725362
relative error of 0x5f376000
1
0
0.0005
0.001
0.0015
0.002
Highcharts.com
    \[\g{change magic constant }m\]

5f376000
but as the bulk of the mass of the error is tucked under the two humps, the 1-norm/mean-metric (and in fact, the first few reasonable L-norms) will tend to disregard the significance of the outer edges and their maximum-bounds. Therefore, these will attempt to decrease the constant by as much as possible.

Finally, let’s look at the maximum relative error after two iterations:

Values
Relative error of 0x5f376000 after 2 newton steps
min y = 2.9531932455029164e-14min y = 2.9531932455029164e-14
max y = 0.000004668678215136168max y = 0.000004668678215136168
relative error of 0x5f376000
1
0
0.000001
0.000002
0.000003
0.000004
0.000005
0.000006
Highcharts.com
    \[\g{change magic constant }m\]

5f376000
This ended up giving us the same magic constant of 0x5f375a1a, which I guess makes sense.

Closing Thoughts
In this article, we’ve come across a couple of different constants: 0x5f400000, 0x54376000, 0x54375a1a. None of these are exactly the 0x543759df that we’ve been searching for, but the last one is really close, which leads me to believe that the discovery of 0x543759df went through a similar process: 1. intuition, and 2. refinement.

The dark art of fast reciprocal square roots have been largely done away with nowadays with faster intrinsics. If anything, this article contributes little technical value besides proposing a half-assed solution to a great puzzle. Nevertheless, I sincerely hope that, having made it this far, you’ve also found some value in this research. Good Luck!

Requested by averageisnothalfway: Differences between a1a and df
Values
Plot of the relative error of a1a and df 1 newton iteration
relative error of 0x5f375a1a
relative error of 0x5f3759DF
Difference in Relative Error (positive segments is where df wins)
1
-0.0005
0
0.0005
0.001
0.0015
0.002
Highcharts.com

 
Posted By Lee
Woosh
Next Post: >>
<< Previous Post: Counting Binary Trees (The Hard Way)

RECENT POSTS

0x5f400000: Understanding Fast Inverse Sqrt the Easy(ish) Way!
Counting Binary Trees (The Hard Way)
Calculating Modified Fibonacci Numbers, Quickly and Exactly
Counting Words
Annoying Mclaurin Series
RECENT COMMENTS

trevor jones on Learn Lua the Hard Way
zersa ogaha on Somewhat Fast
hersa dgaha on Somewhat Fast
Discuz on Learn Lua the Hard Way – Tables
Max on Polygonal Collision Detection
ARCHIVES

October 2014
December 2013
November 2013
October 2013
July 2013
May 2013
March 2013
January 2013
November 2012
August 2012
May 2012
February 2011
CATEGORIES

Algorithm
Article
C
Combinatorics
Lua
Math Problems
Numerical Analysis
OCaml
Programming Language Theory
Python
Reading Diary
Subscribe! Wordpress - Template: Renova - Design by: llowit
2018 Sanity Phailed.me A hearty shoutout to the good folks over at QuickLatex; you guys rock!
 Subscribe! Subscribe  Follow Me! Follow on Twitter  Add Me! Add on Facebook
:)


Main menu
 
WikipediaThe Free Encyclopedia

Search
Donate
Create account
Log in

Personal tools
Contents  hide
(Top)
History
Motivation
Overview of the code
Toggle Overview of the code subsection
Worked example
Avoiding undefined behavior
Algorithm
Toggle Algorithm subsection
Floating-point representation
Aliasing to an integer as an approximate logarithm
First approximation of the result
Newton's method
Accuracy
Subsequent improvements
Toggle Subsequent improvements subsection
Magic number
Zero finding
Obsolescence
See also
Notes
References
Toggle References subsection
Bibliography
Further reading
External links
Fast inverse square root

12 languages
Article
Talk
Read
Edit
View history

Tools
Appearance  hide
Text

Small

Standard

Large
Width

Standard

Wide
Color (beta)

Automatic

Light

Dark
This is a good article. Click here for more information.
From Wikipedia, the free encyclopedia

Lighting and reflection calculations, as in the video game OpenArena, use the fast inverse square root code to compute angles of incidence and reflection.
Fast inverse square root, sometimes referred to as Fast InvSqrt() or by the hexadecimal constant 0x5F3759DF, is an algorithm that estimates 
1
x
{\textstyle {\frac {1}{\sqrt {x}}}}, the reciprocal (or multiplicative inverse) of the square root of a 32-bit floating-point number 
x
{\displaystyle x} in IEEE 754 floating-point format. The algorithm is best known for its implementation in 1999 in Quake III Arena, a first-person shooter video game heavily based on 3D graphics. With subsequent hardware advancements, especially the x86 SSE instruction rsqrtss, this algorithm is not generally the best choice for modern computers,[1] though it remains an interesting historical example.[2]

The algorithm accepts a 32-bit floating-point number as the input and stores a halved value for later use. Then, treating the bits representing the floating-point number as a 32-bit integer, a logical shift right by one bit is performed and the result subtracted from the number 0x5F3759DF, which is a floating-point representation of an approximation of
2
127
{\textstyle {\sqrt {2^{127}}}}.[3] This results in the first approximation of the inverse square root of the input. Treating the bits again as a floating-point number, it runs one iteration of Newton's method, yielding a more precise approximation.

History
William Kahan 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.[4] In the late 1980s, Cleve Moler at Ardent Computer learned about this technique[5] 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 3dfx Interactive circa 1994.[6][7]

Jim Blinn demonstrated a simple approximation of the inverse square root in a 1997 column for IEEE Computer Graphics and Applications.[8] Reverse engineering of other contemporary 3D video games uncovered a variation of the algorithm in Activision's 1997 Interstate '76.[9]

Quake III Arena, a first-person shooter video game, was released in 1999 by id Software and used the algorithm. Brian Hook may have brought the algorithm from 3dfx to id Software.[6] A discussion of the code appeared on the Chinese developer forum CSDN in 2000,[10] and Usenet and the gamedev.net forum spread the code widely in 2002 and 2003.[11] Speculation arose as to who wrote the algorithm and how the constant was derived; some guessed John Carmack.[7] Quake III's full source code was released at QuakeCon 2005, 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 Slashdot.[6]

In 2007 the algorithm was implemented in some dedicated hardware vertex shaders using field-programmable gate arrays (FPGA).[12][13]

Motivation

Surface normals are used extensively in lighting and shading calculations, requiring the calculation of norms for vectors. A field of vectors normal to a surface is shown here.

A two-dimensional example of using the normal 
C
{\displaystyle C} to find the angle of reflection from the angle of incidence; in this case, on light reflecting from a curved mirror. The fast inverse square root is used to generalize this calculation to three-dimensional space.
The inverse square root of a floating point number is used in digital signal processing to normalize a vector, scaling it to length 1 to produce a unit vector.[14] For example, computer graphics programs use inverse square roots to compute angles of incidence and reflection for lighting and shading. 3D graphics 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.[7] This was troublesome for 3D graphics programs before the advent of specialized hardware to handle transform and lighting. Computation of square roots usually depends upon many division operations, which for floating point numbers are computationally expensive. The fast inverse square generates a good approximation with only one division step.

The length of the vector is determined by calculating its Euclidean norm: the square root of the sum of squares of the vector components. 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-dimensional space, so
v
{\displaystyle {\boldsymbol {v}}} would be a vector 
(
v
1
,
v
2
,
v
3
)
{\displaystyle (v_{1},v_{2},v_{3})}.

‖
v
‖
=
v
1
2
+
v
2
2
+
v
3
2
{\displaystyle \|{\boldsymbol {v}}\|={\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}}
is the Euclidean norm of the vector.

v
^
=
v
‖
v
‖
{\displaystyle {\boldsymbol {\hat {v}}}={\frac {\boldsymbol {v}}{\left\|{\boldsymbol {v}}\right\|}}}
is the normalized (unit) vector. Substituting 
‖
v
‖
{\textstyle \|{\boldsymbol {v}}\|}, the equation can also be written as:

v
^
=
v
v
1
2
+
v
2
2
+
v
3
2
{\displaystyle {\boldsymbol {\hat {v}}}={\frac {\boldsymbol {v}}{\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}}}
which relates the unit vector to the inverse square root of the distance components. The inverse square root can be used to compute 
v
^
{\displaystyle {\boldsymbol {\hat {v}}}} because this equation is equivalent to

v
^
=
v
1
v
1
2
+
v
2
2
+
v
3
2
{\displaystyle {\boldsymbol {\hat {v}}}={\boldsymbol {v}}\,{\frac {1}{\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}}}
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}}.

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.

Overview of the code
The following code is the fast inverse square root implementation from Quake III Arena, stripped of C preprocessor directives, but including the exact original comment text:[15]

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;
}
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. Common software methods in the early 1990s drew approximations from a lookup table.[16] 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.[17] The algorithm was designed with the IEEE 754-1985 32-bit floating-point specification in mind, but investigation from Chris Lomont showed that it could be implemented in other floating-point specifications.[18]

The advantages in speed offered by the fast inverse square root trick came from treating the 32-bit floating-point word[note 1] as an integer, then subtracting it from a "magic" constant, 0x5F3759DF.[7][19][20][21] This integer subtraction and bit shift results in a bit pattern which, when re-defined 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 Newton's method; however, it is much slower and less accurate than using the SSE instruction rsqrtss on x86 processors also released in 1999.[1][22]

Worked example
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:

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)
Interpreting as IEEE 32-bit representation:

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
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 single iteration of Newton's method, the final result is 
y
=
2.52549
{\displaystyle y=2.52549}, an error of only 0.17%.

Avoiding undefined behavior
According to the C standard, reinterpreting a floating point value as an integer by casting then dereferencing the pointer to it is not valid (undefined behavior). Another way would be to place the floating point value in an anonymous union containing an additional 32-bit unsigned integer member, and accesses to that integer provides a bit level view of the contents of the floating point value. However, type punning through a union is also undefined behavior in C++.

# include <stdint.h> // uint32_t

float Q_rsqrt(float number)
{
  union {
    float    f;
    uint32_t i;
  } conv = { .f = number };
  conv.i  = 0x5f3759df - (conv.i >> 1);
  conv.f *= 1.5F - (number * 0.5F * conv.f * conv.f);
  return conv.f;
}
In modern C++, the recommended method for implementing this function's casts is through C++20's std::bit_cast. This also allows the function to work in constexpr context:

import std;

constexpr std::float32_t Q_rsqrt(std::float32_t number) noexcept
{
  const auto y = std::bit_cast<std::float32_t>(
    0x5f3759df - (std::bit_cast<std::uint32_t>(number) >> 1));
  return y * (1.5f32 - (number * 0.5f32 * y * y));
}
Algorithm
The algorithm computes 
1
x
{\textstyle {\frac {1}{\sqrt {x}}}} by performing the following steps:

Alias the argument 
x
{\displaystyle x} to an integer as a way to compute an approximation of the binary logarithm 
log
2
⁡
(
x
)
{\textstyle \log _{2}(x)}
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)}
Alias back to a float, as a way to compute an approximation of the base-2 exponential
Refine the approximation using a single iteration of Newton's method.
Floating-point representation
Main article: Single-precision floating-point format
Since this algorithm 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 real number 
x
{\displaystyle x} as a single precision float, the first step is to write 
x
{\displaystyle x} as a normalized binary number:[23]

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}}}
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 significand. 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:

x	=±
2
e
x
(1+
m
x
) 
{\displaystyle {\begin{aligned}x&=\pm 2^{e_{x}}(1+m_{x})\end{aligned}}}
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:[24]

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)
E
x
=
e
x
+
B
{\textstyle E_{x}=e_{x}+B} is the "biased exponent", where 
B
=
127
{\displaystyle B=127} is the "exponent bias"[note 2] (8 bits)
M
x
=
m
x
×
L
{\textstyle M_{x}=m_{x}\times L}, where 
L
=
2
23
{\textstyle L=2^{23}}[note 3] (23 bits)
These fields are then packed, left to right, into a 32-bit container.[25]

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:

x
=
+
2
−
3
(
1
+
0.25
)
{\textstyle x=+2^{-3}(1+0.25)}
and thus, the three unsigned integer fields are:

S
=
0
{\displaystyle S=0}
E
=
−
3
+
127
=
124
=
0111
 
1100
2
{\displaystyle E=-3+127=124=0111\ 1100_{2}}
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}}
these fields are packed as shown in the figure below:


Aliasing to an integer as an approximate logarithm
If 
1
x
{\textstyle {\frac {1}{\sqrt {x}}}} were to be calculated without a computer or a calculator, a table of logarithms 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:

If 
x
{\displaystyle x} is a positive normal number:

x
=
2
e
x
(
1
+
m
x
)
{\displaystyle x=2^{e_{x}}(1+m_{x})}
then

log
2
⁡
(
x
)
=
e
x
+
log
2
⁡
(
1
+
m
x
)
{\displaystyle \log _{2}(x)=e_{x}+\log _{2}(1+m_{x})}
and since 
m
x
∈
[
0
,
1
)
{\displaystyle m_{x}\in [0,1)}, the logarithm on the right-hand side can be approximated by[26]

log
2
⁡
(
1
+
m
x
)
≈
m
x
+
σ
{\displaystyle \log _{2}(1+m_{x})\approx m_{x}+\sigma }
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 optimal approximation (the best in the sense of the uniform norm of the error). However, this value is not used by the algorithm as it does not take subsequent steps into account.


The integer aliased to a floating point number (in blue), compared to a scaled and shifted logarithm (in gray).
Thus there is the approximation

log
2
⁡
(
x
)
≈
e
x
+
m
x
+
σ
.
{\displaystyle \log _{2}(x)\approx e_{x}+m_{x}+\sigma .}
Interpreting the floating-point bit-pattern of 
x
{\displaystyle x} as an integer 
I
x
{\displaystyle I_{x}} yields[note 4]

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}}}
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

log
2
⁡
(
x
)
≈
I
x
L
−
(
B
−
σ
)
.
{\displaystyle \log _{2}(x)\approx {\frac {I_{x}}{L}}-(B-\sigma ).}
First approximation of the result
The calculation of 
y
=
1
x
{\textstyle y={\frac {1}{\sqrt {x}}}} is based on the identity

log
2
⁡
(
y
)
=
−
1
2
log
2
⁡
(
x
)
{\displaystyle \log _{2}(y)=-{\tfrac {1}{2}}\log _{2}(x)}
Using the approximation of the logarithm above, applied to both 
x
{\displaystyle x} and 
y
{\displaystyle y}, the above equation gives:

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)}
Thus, an approximation of 
I
y
{\displaystyle I_{y}} is:

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}}
which is written in the code as

i  = 0x5f3759df - ( i >> 1 );
The first term above is the magic number

3
2
L
(
B
−
σ
)
=
0
x
5
F
3759
D
F
{\displaystyle {\tfrac {3}{2}}L(B-\sigma )={\mathtt {0x5F3759DF}}}
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.[27]

Newton's method
Main article: Newton's method





Relative error between direct calculation and fast inverse square root carrying out 0, 1, 2, 3, and 4 iterations of Newton's root-finding method. Note that double precision is adopted and the smallest representable difference between two double precision numbers is reached after carrying out 4 iterations.
With 
y
{\displaystyle y} as the inverse square root, 
f
(
y
)
=
1
y
2
−
x
=
0
{\displaystyle f(y)={\frac {1}{y^{2}}}-x=0}. The approximation yielded by the earlier steps can be refined by using a root-finding method, a method that finds the zero of a function. The algorithm uses Newton's method: 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}-{\frac {f(y_{n})}{f'(y_{n})}}}, where 
f
′
(
y
n
)
{\displaystyle f'(y_{n})} is the derivative of 
f
(
y
)
{\displaystyle f(y)} at 
y
n
{\displaystyle y_{n}}.[28] For the above 
f
(
y
)
{\displaystyle f(y)},

y
n
+
1
=
y
n
(
3
−
x
y
n
2
)
2
{\displaystyle y_{n+1}={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}}
where 
f
(
y
)
=
1
y
2
−
x
{\displaystyle f(y)={\frac {1}{y^{2}}}-x} and 
f
′
(
y
)
=
−
2
y
3
{\displaystyle f'(y)=-{\frac {2}{y^{3}}}}.

Treating 
y
{\displaystyle y} as a floating-point number, y = y*(threehalfs - x/2*y*y); is equivalent to

y
n
+
1
=
y
n
(
3
2
−
x
2
y
n
2
)
=
y
n
(
3
−
x
y
n
2
)
2
.
{\displaystyle y_{n+1}=y_{n}\left({\frac {3}{2}}-{\frac {x}{2}}y_{n}^{2}\right)={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}.}
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 converge to the inverse square root.[29] For the purposes of the Quake III engine, only one iteration was used. A second iteration remained in the code but was commented out.[21]

Accuracy
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.

Subsequent improvements
Magic number
It is not known precisely how the exact value for the magic number was determined. Chris Lomont developed a function to minimize approximation error 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.[30] 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.[30] He concluded by asking whether the exact value of the original constant was chosen through derivation or trial and error.[31] 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.[32]

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,[33] arriving after an exhaustive search at

conv.i = 0x5F1FFFF9 - ( conv.i >> 1 );
conv.f *= 0.703952253f * ( 2.38924456f - x * conv.f * conv.f );
return conv.f;
A complete mathematical analysis for determining the magic number is now available for single-precision floating-point numbers.[34][35]

Zero finding
Intermediate to the use of one vs. two iterations of Newton's method in terms of speed and accuracy is a single iteration of Halley's method. 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

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),}
where the implementation should calculate 
x
y
n
2
{\displaystyle xy_{n}^{2}} only once, via a temporary variable.

Obsolescence
Subsequent additions by hardware manufacturers have made this algorithm redundant for the most part. For example, on x86, Intel introduced the SSE instruction rsqrtss in 1999. In a 2009 benchmark on the Intel Core 2, this instruction took 0.85ns per float compared to 3.54ns for the fast inverse square root algorithm, and had less error.[1]

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 CORDIC.

See also
Methods of computing square roots § Approximations that depend on the floating point representation
Magic number
Notes
 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.
 
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.
 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.
 Since 
x
{\displaystyle x} is positive, 
S
x
=
0
{\displaystyle S_{x}=0}.
References
 Ruskin, Elan (2009-10-16). "Timing square root". Some Assembly Required. Archived from the original on 2021-02-08. Retrieved 2015-05-07.
 feilipu. "z88dk is a collection of software development tools that targets the 8080 and z80 computers". GitHub.
 Munafo, Robert. "Notable Properties of Specific Numbers". mrob.com. Archived from the original on 16 November 2018.
 "sqrt implementation in fdlibm - See W. Kahan and K.C. Ng's discussion in comments in lower half of this code".
 Moler, Cleve (19 June 2012). "Symplectic Spacewar". MATLAB Central - Cleve's Corner. MATLAB. Retrieved 2014-07-21.
 Sommefeldt, Rys (2006-12-19). "Origin of Quake3's Fast InvSqrt() - Part Two". Beyond3D. Retrieved 2008-04-19.
 Sommefeldt, Rys (2006-11-29). "Origin of Quake3's Fast InvSqrt()". Beyond3D. Retrieved 2009-02-12.
 Blinn 1997, pp. 80–84.
 Peelar, Shane (1 June 2021). "Fast reciprocal square root... in 1997?!".
 "Discussion on CSDN". Archived from the original on 2015-07-02.
 Lomont 2003, p. 1-2.
 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.
 Middendorf 2007, pp. 155–164.
 Blinn 2003, p. 130.
 "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.
 Eberly 2001, p. 504.
 Lomont 2003, p. 1.
 Lomont 2003.
 Lomont 2003, p. 3.
 McEniry 2007, p. 2, 16.
 Eberly 2001, p. 2.
 Fog, Agner. "Lists of instruction latencies, throughputs and micro-operation breakdowns for Intel, AMD and VIA CPUs" (PDF). Retrieved 2017-09-08.
 Goldberg 1991, p. 7.
 Goldberg 1991, pp. 15–20.
 Goldberg 1991, p. 16.
 McEniry 2007, p. 3.
 Hennessey & Patterson 1998, p. 305.
 Hardy 1908, p. 323.
 McEniry 2007, p. 6.
 Lomont 2003, p. 10.
 Lomont 2003, pp. 10–11.
 Matthew Robertson (2012-04-24). "A Brief History of InvSqrt" (PDF). UNBSJ.
 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.
 Moroz et al. 2018.
 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.
Bibliography
Blinn, Jim (July 1997). "Floating Point Tricks". IEEE Computer Graphics & Applications. 17 (4): 80. doi:10.1109/38.595279.
Blinn, Jim (2003). Jim Blinn's Corner: Notation, notation notation. Morgan Kaufmann. ISBN 1-55860-860-5.
Eberly, David (2001). 3D Game Engine Design. Morgan Kaufmann. ISBN 978-1-55860-593-0.
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.
Hardy, Godfrey (1908). A Course of Pure Mathematics. Cambridge, UK: Cambridge University Press. ISBN 0-521-72055-9.
Hennessey, John; Patterson, David A. (1998). Computer Organization and Design (2nd ed.). San Francisco, CA: Morgan Kaufmann Publishers. ISBN 978-1-55860-491-9.
Lomont, Chris (February 2003). "Fast Inverse Square Root" (PDF). Retrieved 2009-02-13.
McEniry, Charles (August 2007). "The Mathematics Behind the Fast Inverse Square Root Function Code" (PDF). Archived from the original (PDF) on 2015-05-11.
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. Open access icon
Striegel, Jason (2008-12-04). "Quake's fast inverse square root". Hackszine. O'Reilly Media. Archived from the original on 2009-02-15. Retrieved 2013-01-07.
IEEE Computer Society (1985). "754-1985 - IEEE Standard for Binary Floating-Point Arithmetic". Institute of Electrical and Electronics Engineers.
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.
Further reading
Kushner, David (August 2002). "The wizardry of Id". IEEE Spectrum. 39 (8): 42–47. doi:10.1109/MSPEC.2002.1021943.
External links
0x5f3759df, further investigations into accuracy and generalizability of the algorithm by Christian Plesner Hansen
Origin of Quake3's Fast InvSqrt()
Quake III Arena source code (archived in Software Heritage)
Implementation of InvSqrt in DESMOS
"Fast Inverse Square Root — A Quake III Algorithm" (YouTube)
vte
Quake series
vte
id Software (games)
Categories: Quake (series)Source codeRoot-finding algorithms
This page was last edited on 12 November 2024, at 17:23 (UTC).
Text is available under the Creative Commons Attribution-ShareAlike 4.0 License; additional terms may apply. By using this site, you agree to the Terms of Use and Privacy Policy. Wikipedia® is a registered trademark of the Wikimedia Foundation, Inc., a non-profit organization.
Privacy policyAbout WikipediaDisclaimersContact WikipediaCode of ConductDevelopersStatisticsCookie statementMobile view
Wikimedia FoundationPowered by MediaWiki