Fast Square Root Approximation
Published on 17.01.2026
Almost a year ago I came across an old article named Understanding Quake's Fast Inverse Square Root. It blew my mind. A very fascinating approach to cut down the execution time by a good margin.
# $ ./a.out
number: 1337
+ result of math.h sqrt 36.565
* math.h took 15 ticks(1.5e-05s)
+ result of fast sqrt 36.5668
* fast sqrt took 1 ticks(1e-06s)
Note that this is an approximation, which was good enough for quake. In this post, we'll re-create this algorithm from scratch.
The original code looks like the following.
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
#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}
Not that obvious, huh?
IEEE 754 Standard
"The IEEE Standard for Floating-Point Arithmetic (IEEE 754) is a technical standard for floating-point arithmetic established in 1985 by the Institute of Electrical and Electronics Engineers (IEEE). The standard addressed many problems found in the diverse floating-point implementations that made them difficult to use reliably and portably" -Wikipedia
Floating point numbers follow IEEE 754 standard for efficient memory utilization. Consider a number 43.4, how do you think a computer would store this number in 32 bits of memory?
How about splitting the number into 43 and 0.3, 16 bits each? Well that would allow us to have numbers only upto 2**16 = 65536, which is a waste of memory. But there's an efficient way of storing larger range of numbers with IEEE 754 representation.
Consider the same number 43.3.
[cling]\$ float n = 43.3
(float) 43.3000f
[cling]\$ *(long *)&n
(long) 1110258483
[cling]\$ #include <bitset>
[cling]\$ std::cout << std::bitset<32> (*(long *)&n);
01000010001011010011001100110011
43.3 is represented as 1110258483 in memory, but why? Short answer is to be more efficient with space.
The idea is to represent a floating point number as a floating point arithmetic expression.
\[ n = {s \times b ^ {e}} \]
Here s is the Significand, b is the Base, and e is the Exponent.
An example,
\[ 12.345 = {12345 * 10 ^ {-3}} \]
So let's try this for the number 43.3 with the base as 2.
>>> bin(43)
'0b101011'
>>> n = 0.3
>>> ''.join(map(str, [int((n := (n - (1 * (n > 1))) * 2) > 1) for i in range(32)]))
'01001100110011001100110011001100'
So the number in binary looks like the following.
# 43.3
101011.01001100110011001100110011001100
Now let's convert this to a floating point arithmetic expression, by moving the decimal point to the left.
\[ 1.0101101001100110011001100110011001100 \times 2 ^ {5} \]
Where s = 1.0101101001100110011001100110011001100, b = 2, and e = 5.
If we wanna store this representation in 32 bits, we can ignore the constants, like the base b. But we can also ignore the first digit of s, because we are dealing with Base 2 and the first digit is always 1, hence we can omit it.
So now s - 1 becomes the Mantissa M, and E is the Exponent. In memory the 32 bit number looks like the following.
Sign(1 bit) Exponent(8 bits) Mantissa(23 bits)
# ######## #######################
- Sign S is 1 if the number is negative
- Exponent E is 127 + e. We add 127 so we can have negative exponents. So in our case for the number 43.3, we had the e = 5, but we add 127 to it, so it becomes E = 132
- Mantissa M will be s - 1, where s is the significand. Note that Mantissa is about 23 bits for a 32 bit floating point number, so numbers beyond 23 bits are truncated.
For the number 43.3, the IEEE 754 representation looks like the following.
Sign(1 bit) Exponent(8 bits) Mantissa(23 bits)
0 10000100 01011010011001100110011 // Binary
0 132 2962227 // Decimal
The decimal representation of this number is,
>>> 0b01000010001011010011001100110011
1110258483
We can confirm that in C++.
[cling]\$ float n = 43.3
(float) 43.3000f
[cling]\$ *(long *)&n
(long) 1110258483
With S, M and E, we can get the original floating point number using the following expression.
\[ n = {(-1) ^ S} \times ({1 + {M \over {2 ^ {23}}}}) \times {2 ^ {E-127}} \]
Solving for n,
\[ n = {(-1) ^ 0} \times ({1 + {2962227 \over {2 ^ {23}}}}) \times {2 ^ {132-127}} \]
\[ n = 43.29999923706055 \]
Note I won't be covering rounding or other nuances in IEEE 754 standard, but this should give us a pretty good idea of what's going on under the hood.
We can also get the IEEE 754 representation of the number with S, M and E using,
\[ B = {2 ^ {31}} \times S + {2 ^ {23}} \times E + M \]
Solving for B,
\[ B = {2 ^ {31}} \times 0 + {2 ^ {23}} \times 132 + 2962227 \]
\[ B = 1110258483 \]
Now that we have a general understand of how floating point numbers are stored, we can start with solving for square roots.
Simplifying with Logarithms
Let's simplify the following expression.
\[ y = ({1 + {M \over {2 ^ {23}}}}) \times {2 ^ {E-127}} \]
Note: From here, I'm only considering positive normal finite IEEE 754 binary32 values.
S is omitted because we're only considering positive numbers.
Before we find \( \sqrt{y} \), we can simplify the equation a bit. Since we are dealing with exponents, we can try applying logarithm and see what happens.
\[ \log_2 (y) = \log_2 (({1 + {M \over {2 ^ {23}}}}) \times {2 ^ {E-127}}) \]
\[ \log_2 (y) = \log_2 (({1 + {M \over {2 ^ {23}}}})) + \log_2({2 ^ {E-127}}) \]
\[ \log_2 (y) = \log_2 (({1 + {M \over {2 ^ {23}}}})) + {E-127} \]
We simplified it to a degree, but there's still \( \log_2 \) in the equation. As it turns out that for \( 0 < {M \over {2 ^ {23}}} < 1 \), there's an approximation.
\[ {\log_2 ({1 + x}) \approx x}, {0 < {x} < 1} \]
Here's a function plot for \( f(x) = \log_2 ({1 + x}) \), for \( x \in [0,1) \).
Kinda looks linear doesn't it? Let's plot one for \( f(x) = x \).
It's slightly off, but we can add a bias \( \sigma \) to offset the \( f(x) = x \) such that it would best fit \( \log_2 ({1 + x}) \), meaning with least error.
\[ \log_2 (y) = {M \over {2 ^ {23}}} + \sigma + {E-127} \]
\[ \log_2 (y) = {1 \over 2^{23}} ({E \times {2 ^ {23}}} + M) + \sigma - 127 \]
\( {E \times {2 ^ {23}} + M} \) should look familiar because that's how we get the IEEE 754 representation. And also it turns out, \( \sigma = 0.0430 \) minimizes the worst-case error for points \( [0,1) \).
Solving for Square Roots
Now let's solve for the square root, we'll call it \( r \).
\[ r = \sqrt{y} \]
\[ \log(r) = \log(\sqrt{y}) \]
\[ \log(r) = {1 \over 2}\log(y) \]
\[ {1 \over {2^{23}}}(i_r) + \sigma - 127 = {1 \over 2} ({1 \over {2^{23}}}(i_y) + \sigma - 127), i_x = ({E_x \times {2 ^ {23}}} + M_x), x = r | y \]
\[ {i_r \over 2^{23}} = {i_y \over 2^{24}} + {\sigma - 127 \over 2} - \sigma + 127 \]
\[ {i_r \over 2^{23}} = {i_y \over 2^{24}} + {\sigma - 127 - 2\sigma + 254 \over 2} \]
\[ {i_r} = 2^{23}({i_y \over 2^{24}} + {127 - \sigma \over 2}) \]
\[ {i_r} = ({127 - \sigma}){2^{22}} + {i_y \over 2} \]
\[ {i_r} = ({127 - 0.0430}){2^{22}} + {i_y \over 2}, \sigma = 0.0430 \]
\[ {i_r} = 532496253 + {i_y \over 2} \]
>>> round((127 - 0.0430) * (2**22))
532496253
>>> hex(_)
'0x1fbd3f7d'
The solution can be translated to the following code.
i = 0x1fbd3f7d + (i >> 1);
Here instead of division i/2, we shift right by one i >> 1 to half the number. We do this because division can be slow. But these days compilers are intelligent enough to use bit shifts and multiplication to divide a number if the divisor is known at compile time.
int optimized(int n) {
int x = n / 2, y = n / 3;
return x + y;
}
optimized(int):
mov eax, edi
shr eax, 31
add eax, edi
sar eax
mov edx, eax
movsx rax, edi
sar edi, 31
imul rax, rax, 1431655766
shr rax, 32
sub eax, edi
add eax, edx
ret
Let's test the square root solution for the number 43.3.
[cling]\$ float n = 43.3
(float) 43.3000f
[cling]\$ long i = *(long *)&n
(long) 1110258483
[cling]\$ i = 0x1fbd3f7d + (i >> 1)
(long) 1087625494
[cling]\$ *(float *)&i
(float) 6.62025f
>>> 6.62025 ** 2
43.8277100625
That's a pretty good approximation, but we can make it better with Newton–Raphson method.
Error minimization with Newton–Raphson method
"In numerical analysis, Newton's method, also known as the Newton–Raphson method, named after Isaac Newton and Joseph Raphson, is a root-finding algorithm which produces successively better approximations to the roots (or zeroes) of a real-valued function." -Wikipedia
The idea is to get closer to the approximation by starting off with calculating \( f(x) \) for some initial value of \( x \), then finding its derivative to get closer to the root of \( x \).
Here's an example, consider \( f(x) = (2x-1)^{3} + 0.3 \), and initial \( x_0 = 0.76 \).
\[ x_{n+1} = x_n - {f(x_n) \over f'(x_n)} \]
\[ x_{1} = x_0 - {{(2x-1)^{3} + 0.3} \over {24x^{2} - 24x + 6}} \]
\[ x_{1} = 0.76 - {{0.440608} \over {1.622399999999999}} \]
\[ x_{1} = 0.4884220907297829 \]
We can now plug \( x_{1} = 0.4884220907297829 \) back into the newton's method equation to get \( x_{2} \), and so on, it's iterative.
In our case, \( f(x) = y_a^2 - n \), where \( y_r = \sqrt{n} \) and \( y_r \approx y_a \). Here \( y_r \) is the result, and \( y_a \) is our approximate value.
\[ x_{n+1} = x_n - {f(x_n) \over f'(x_n)} \]
\[ y_{n+1} = y_n - {y_n^2 - n \over 2y_n} \]
\[ y_{n+1} = {2y_n^2 - y_n^2 + n \over 2y_n} \]
\[ y_{n+1} = 0.5 \times ({y_n^2 + n \over y_n}) \]
We can confirm the minimization of the error like the following.
[cling]\$ #include <cmath>
[cling]\$ pow(6.62025f, 2) // 6.62025f is our initial approximation
(double) 43.827710
[cling]\$ 0.5 * ((pow(6.62025f, 2) + 43.3)/6.62025f)
(double) 6.5803943
[cling]\$ pow(6.5803943f, 2)
(double) 43.301589
There we go, we have a better approximation.
Is our magic number the best one?
Our constant 0x1fbd3f7d looks nothing like Quake's 0x5f3759df. Why? Because we solved a different equation. We solved \( \sqrt{y} \), which halves the exponent, so C + (i >> 1). Quake solved \( 1/\sqrt{y} \), which negates and halves it, so C - (i >> 1). Different equation, different magic number.
But is our number even the best one for our equation? Turns out, no. Our \( \sigma = 0.0430 \) minimizes the error of the \( \log_2 \) approximation, but not the error of the whole bit hack. Quake has the same problem, Chris Lomont later found 0x5f37642f which beats Quake's constant before the Newton step.
We can search for a better constant by brute force. Try a bunch of constants, and for each one measure the largest relative error across 200,000 samples over \( x \in [1, 4) \). We only need this range because the error pattern repeats every factor of 4, the trick only cares about exponent parity.
float sqrt_guess(float x, uint32_t C) {
uint32_t i;
memcpy(&i, &x, 4);
i = C + (i >> 1);
float y;
memcpy(&y, &i, 4);
return y;
}
double max_rel_err(uint32_t C) {
double maxe = 0;
for (int k = 0; k < 200000; k++) {
float x = 1.0f + 3.0f * k / 200000;
double e = fabs(sqrt_guess(x, C) / sqrt((double)x) - 1.0);
if (e > maxe) maxe = e;
}
return maxe;
}
Sweep around our derived constant, first in coarse steps of 512, then step of 1 around the winner.
== sqrt: i = C + (i >> 1) ==
derived (sigma=0.0430) 0x1fbd3f7d max rel err = 4.55%
brute forced 0x1fbb4f2e max rel err = 3.47%
brute forced + 1 newton 0x1fbb67b6 max rel err = 0.060%
Our derived constant was decent, but the brute forced one is ~24% more accurate. Same code, just a better number.
To make sure the search isn't broken, I ran it on Quake's equation too. It finds 0x5f37642e, just one off from Lomont's 0x5f37642f. Close enough.
So there's no single magic number, there's a whole family of them. The math gets you in the neighborhood, brute force finds the best house.
Now, putting it all together with the constant that performed the best in our search, the final code looks like the following. Note that since we do a Newton iteration at the end, we use 0x1fbb67b6 and not the raw-guess winner, the Newton step reshapes the error so the best starting constant shifts.
#include <cstdint>
#include <iostream>
#include <bits/stdc++.h>
using namespace std;
int main() {
float n = 43.3;
// use u32 because long can be 64 bits but float is 32 bits
uint32_t i;
std::memcpy(&i, &n, sizeof(float));
// Solved equation for square roots (brute forced constant)
i = 0x1fbb67b6 + (i >> 1);
float y;
std::memcpy(&y, &i, sizeof(float));
// Newton-Rapson iteration
y = (((y * y) + n)/(y)) * 0.5;
std::cout << y << std::endl;
}
Note: Instead of the pointer magic (which produced undefined behaviours), I've used memcpy to get the raw IEEE 754 representation of a floating point number.
Hope this was insightful.
- Quake 3 Fast Inverse Square Root Source Code
- Methods of Computing Square Roots
- Fast Inverse Square Root — A Quake III Algorithm
- Fast Inverse Square Root — Chris Lomont
Feel free to reach out on twitter.