Fast inverse square root


Fast inverse square root, sometimes referred to as ' or by the hexadecimal constant ', is an algorithm that estimates, the reciprocal of the square root of a 32-bit floating-point number 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, though it remains an interesting historical example.
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, which is a floating-point representation of an approximation of. 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

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. In the late 1980s, Cleve Moler at Ardent Computer learned about this technique 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.
Jim Blinn demonstrated a simple approximation of the inverse square root in a 1997 column for IEEE Computer Graphics and Applications. Reverse engineering of other contemporary 3D video games uncovered a variation of the algorithm in Activision's 1997 Interstate '76.
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. A discussion of the code appeared on the Chinese developer forum CSDN in 2000, and Usenet and the gamedev.net forum spread the code widely in 2002 and 2003. Speculation arose as to who wrote the algorithm and how the constant was derived; some guessed John Carmack. Quake IIIs 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.
In 2007 the algorithm was implemented in some dedicated hardware vertex shaders using field-programmable gate arrays.

Motivation

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. 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. 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 would be a vector. Then,
is the Euclidean norm of the vector, and the normalized vector is
where the fraction term is the inverse square root of.
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 C code is the fast inverse square root implementation from Quake III Arena, stripped of C preprocessor directives, but including the exact original comment text:

float Q_rsqrt

At the time, the general method to compute the inverse square root was to calculate an approximation for, 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. 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. 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.
The advantages in speed offered by the fast inverse square root trick came from treating the 32-bit floating-point word as an integer, then subtracting it from a "magic" constant,. 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.

Worked example

As an example, the number can be used to calculate. 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:
0101_1111_0011_0111_0101_1001_1101_1111 The magic number 0x5F3759DF
0100_0000_0010_0111_0101_1001_1101_1111 The result of 0x5F3759DF -
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, which has an error of about 3.4%. After one iteration of Newton's method, the final result is, 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 can cause undefined behavior in case the sizes of integer and float do not match on the given architecture. This can be avoided by using alternative type punning techniques such as C's unions or C++20's std::bit_cast.

Algorithm

The algorithm computes by performing the following steps:
  1. Alias the argument to an integer as a way to compute an approximation of the binary logarithm
  2. Use this approximation to compute an approximation of
  3. Alias back to a float, as a way to compute an approximation of the base-2 exponential
  4. Refine the approximation using a single iteration of Newton's method.

    Floating-point representation

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 as a single precision float, the first step is to write as a normalized binary number:
where the exponent is an integer, and 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:
where means, so. From this form, three unsigned integers are computed:
  • , the "sign bit", is if is positive and negative or zero
  • is the "biased exponent", where is the "exponent bias"
  • , where
Thus: and.
These fields are then packed, left to right, into a 32-bit container.
As an example, consider again the number. Normalizing yields:
and thus, the three unsigned integer fields are:
these fields are packed as shown in the figure below:
The number is represented in binary as:
Also, since this algorithm works on real numbers, is only defined for. The code thus assumes and.
The number, given to calculate the square root, could be rewritten as:

Aliasing to an integer as an approximate logarithm

If were to be calculated without a computer or a calculator, a table of logarithms would be useful, together with the identity, which is valid for every base. 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 is a positive normal number:
then
and since, the logarithm on the right-hand side can be approximated by
where is a free parameter used to tune the approximation. For example, yields exact results at both ends of the interval, while yields the optimal approximation. However, this value is not used by the algorithm as it does not take subsequent steps into account.
Thus there is the approximation
Interpreting the floating-point bit-pattern of as an integer yields
It then appears that is a scaled and shifted piecewise-linear approximation of, as illustrated in the figure on the right. In other words, is approximated by