Fast Square Root Approximation

Created: Jun 29, 2025

Last modified: Jul 07, 2025

Back in 1999, Quake III Arena shipped with a function that became a legend: the fast inverse square root, which computes $1/\sqrt{x}$ with a weird hexadecimal constant and some "evil" bit hacking. Here I take the same trick and bend it to compute the plain square root $\sqrt{x}$. I used it to get 0 ms runtimes in both C and Python on this LeetCode problem.

Let's be honest: this is a learning experiment, not something you need today. Modern CPUs have instructions like rsqrtss that do this in hardware (although they don't always give the same results on different CPUs). But understanding why the trick works is a lot of fun.

# Floats are secretly logarithms

Everything starts with how a 32-bit float is stored: 1 sign bit, 8 bits of exponent $E$ and 23 bits of mantissa $M$. Now, the curious part: if you take those 32 bits and read them as a plain integer, which we'll call $I(x)$, that integer turns out to be almost a scaled and shifted logarithm of the number:

The small correction $\sigma$ tunes the approximation. I won't go through every step of why this holds; this article on breq.dev explains it wonderfully. What matters to us: logarithms turn powers into multiplications, and a square root is just a power of ½.

# Where does 0x5F3759DF come from?

The famous constant in Quake's code comes straight out of that formula, with $\sigma = 0.0450466$:

Mystery solved. And if that constant was built for $1/\sqrt{x}$, we can build our own for $\sqrt{x}$ the same way.

# Our own magic number

Let's call the result we want $\Gamma = \sqrt{x}$. Taking logarithms, the square root becomes a simple halving:

Now we replace each logarithm with its "read the bits as an integer" approximation:

And solve for $I(\Gamma)$, the bits of our answer:

Look at that: the first part doesn't depend on $x$ at all, so we can compute it once. That's our magic number:

And halving an integer is just shifting its bits one place to the right. So the whole square root guess fits in one line:

i = 0x1FBD1DF5 + (i >> 1)

Compare it with Quake's line, i = 0x5f3759df - (i >> 1). Same idea, but with a plus instead of a minus, because we want $x^{1/2}$ and not $x^{-1/2}$, and a magic number three times smaller, because the constant gets multiplied by ½ instead of 3/2.

# See it in action

Type any number and follow the trick step by step: first its 32 bits, then the bits after adding the magic number to half of them, and finally how good that first guess $y_0$ is. Try big numbers, small ones, perfect squares... The guess is always within about 5% of the real root.

Notice the table: after the first guess, the error doesn't just shrink, it collapses. That's the second half of the trick.

# Polishing the guess: Newton... or Heron?

A 5% error isn't good enough to get the integer part right, so we refine the guess with Newton's method. To use it, we need an equation whose root is the number we're after. If $y = \sqrt{x}$, then:

Newton's formula says: from a guess $y_n$, follow the tangent line down to zero to get a better one. Plugging in our $f$ and simplifying:

And here's a fun surprise: that last line is much older than Newton. It's the Heron's method (also called the Babylonian method), described by Heron of Alexandria almost 2000 years ago. The idea is beautifully simple: if $y$ is a bit too big, then $x / y$ is a bit too small, and the real root is somewhere in between. So take the average, and repeat.

It converges quadratically: roughly, every step doubles the number of correct digits. For $x = 2$, the error goes from 0.06 to 0.001, then to 0.0000006, then to 0.0000000000001. That's why one or two steps are enough when you start from a good guess, and that's exactly what the bit hack gives us. In code, each step is a single line:

y = 0.5 * (y + x / y)

# From C pointers to Python

The original C code reads a float's bits with a pointer cast, the famous "evil floating point bit level hacking":

i  = * ( long * ) &y;    // float bits -> integer
y  = * ( float * ) &i;   // integer -> float bits

Python has no pointers, but the struct module can pack a number into its 4 raw bytes and unpack it back, which gets us the same result:

# float bits -> integer
i = int.from_bytes(struct.pack("!f", x), byteorder="big", signed=False)
# integer -> float bits
y = struct.unpack("!f", int(i).to_bytes(4, byteorder="big", signed=False))[0]

# The final code

Putting the pieces together: the bit hack for the first guess, a couple of Heron steps, and the integer part. In Python:

import struct

class Solution:

    def mySqrt(self, x: int) -> int:
        bytes_values = struct.pack("!f", x)
        i = int.from_bytes(bytes_values, 'big', signed=False)
        i = 0x1FBD1DF5 + (i >> 1)
        y = int(i).to_bytes(4, byteorder="big", signed=False)
        y = struct.unpack("!f", y)[0]
        y = 0.5 * (y + x / y)
        y = 0.5 * (y + x / y)
        # y = 0.5 * (y + x / y)
        # y = 0.5 * (y + x / y)
        # y = 0.5 * (y + x / y)
        return int(y)

And in C:

#include <stdio.h>
#include <stdint.h>
#include <math.h>

int mySqrt(int x) {
    long i;
    float y, _x;
    _x = 1.0f * x;
    y  = _x;
    i  = * ( long * ) &_x;
    i  = 0x1FBD2B25 + ( i >> 1 );
    y  = * ( float * ) &i;
    y  = 0.5 * (y + _x / y);
    y  = 0.5 * (y + _x / y);
    //y  = 0.5 * (y + _x / y);
    //y  = 0.5 * (y + _x / y);
    //y  = 0.5 * (y + _x / y);
    return (int)y; 
}

Both run in O(1) time and O(1) space: a fixed number of operations, no matter how big $x$ is.

# What happened on LeetCode

# Why is C less accurate?

These graphs show the absolute error of each version over the same inputs and the same number of steps. First C:

And Python:

Python wins clearly, and the reason is precision. Both versions start from the same 32-bit guess, but then:

With numbers close to $2^{31}$, the root is around 46340, and 7 digits leave very little room for decimals. If the real root is something like 46340.9999, a float can round it up to 46341, and the integer part comes out one too high. That's exactly what happened in the three failing cases.

# Wrapping up

This little experiment mixes two ideas from very different eras: a 1990s hack that reads a float's bits as a logarithm to get a decent guess for free, and a 2000-year-old averaging method that turns that guess into an exact answer in a couple of steps. Neither is very useful alone; together, they're surprisingly fast.