Rethinking 0x5F3759DF for GPUs

The old spell was 0x5F3759DF. We recast it for GPUs, FP16/BF16, and RMSNorm, and discovered that modern magic may need a lookup table.

Kirato Yoshihara

Introduction

The most famous magic number in game programming is probably 0x5F3759DF. It appeared in the fast inverse square root routine popularized by Quake III Arena [1], where a floating-point number was reinterpreted as an integer, shifted, subtracted from a mysterious constant, and refined with one Newton step. For its time, it was a beautiful scalar CPU trick: approximate 1/\sqrt{x} without paying the full cost of a square root.

But GPUs changed the question. Modern NVIDIA hardware already has fast reciprocal-square-root instructions, and ML workloads often run in FP16, BF16, or even lower precision. So instead of asking whether the original Quake trick is still useful as-is, we wanted to ask a slightly different question: what should a Q_rsqrt-style approximation look like in the GPU and low-precision ML era?

The short answer is that FP32 is no longer where the old trick shines. On A100, the classic bit hack does not clearly beat native GPU math. But in FP16 and BF16, the story becomes more interesting. Because the input spaces are small enough to search directly, we can find format-specific magic constants, then improve them with tiny mantissa-conditioned lookup tables. The result is not a replacement for native rsqrt in accuracy, but a compact accuracy-throughput trade-off for low-precision GPU kernels.

What Q_rsqrt Was Doing

The original routine is short enough to fit on a postcard.

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

The strange line is not the Newton step. It is this one:

i = 0x5f3759df - (i >> 1);

A positive float stores a number using an exponent and a mantissa. Very roughly, the exponent tells us the scale of the number, which means it behaves a little like \log_2 x. If we reinterpret the float bits as an integer, the integer is not equal to \log_2 x, but it carries a surprisingly useful affine approximation to it.

float32 bit layout

    sign | exponent | mantissa
     1   |    8     |    23

    x ≈ (1 + mantissa) × 2^(exponent - 127)

Here, 127 is the exponent bias for FP32. As we will see later, FP16 and BF16 use different biases, but the exact same core mechanism applies. To compute 1/\sqrt{x}, we want something like - \frac{1}{2}\log_2 x in log space. The shift i >> 1 approximately divides the encoded exponent information by two, and subtracting from the magic constant flips the sign and places the result back into the range of a floating-point number. Reinterpreting those bits as a float gives an initial guess.

That initial guess is only rough, so the code applies one Newton refinement:

y = y * (threehalfs - x2 * y * y);

This step is what turns a crude bit-level estimate into a useful approximation. The magic constant chooses the starting point. Newton's method cleans it up. For a more detailed mathematical derivation of the bit-level approximation, see McEniry's note [2].

FP32 on A100: The Old Hack No Longer Wins

The first question was whether the original trick still matters for FP32 on modern GPUs. We compared five variants on an NVIDIA A100-SXM4-40GB: plain 1.0f / sqrtf(x), CUDA rsqrtf(x), the original Quake constant 0x5f3759df with one Newton step, Lomont's improved constant 0x5f375a86 [3] with one Newton step, and native rsqrt.approx.ftz.f32 followed by one Newton step. Each kernel was warmed up before timing, then measured over 10 runs.

For the improved constant, we used Lomont's constant, which was obtained by numerically searching for a better constant after one Newton step. Lomont also discusses why the best initial approximation is not necessarily the best constant after Newton refinement.

Table 1. FP32 inverse-square-root baselines on an NVIDIA A100. Errors are measured relative to torch.rsqrt. CUDA rsqrtf is the native baseline, so its error is zero by definition. In this setup, 1.0f / sqrtf(x) produced the same outputs as torch.rsqrt. The Quake-style variants match native throughput but lose accuracy.

The result was not subtle. On A100, the classic Quake-style bit hack does not clearly beat native GPU math in FP32. The Quake and Lomont variants run at roughly the same speed as sqrtf and rsqrtf, but their maximum relative error is around 1.75 \times 10^{-3}.

For context, we also measured the native GPU rsqrt accuracy against a float64 reference in the later low-precision experiments. In that setting, native GPU rsqrt has maximum relative error around 1.25 \times 10^{-7}, while Q_rsqrt-style approximations after one Newton step are in the 10^{-3} range. So native rsqrt is orders of magnitude more accurate.

This does not make the old idea useless. It just changes where the interesting question lives. The broader idea of using floating-point bit patterns to form square-root or inverse-square-root initial guesses also appears in older numerical software notes, including the fdlibm square-root notes by Kahan and Ng [4]. For FP32, the scalar CPU-era trick no longer has an obvious advantage. But for FP16, BF16, and ML-shaped kernels, the input space is smaller, the precision target is different, and the design space opens up again.

Searching FP16/BF16 Magic Constants

The FP32 space is too large to brute-force casually, but FP16 and BF16 are small enough to search directly. For each format, we enumerated all positive normal inputs, tried every 16-bit candidate constant C, and evaluated the Q_rsqrt-style approximation:

// Conceptual brute-force search for a 16-bit floating-point format
// The input set contains all positive normal FP16 or BF16 values
float best_error = INFINITY;
uint16_t best_C = 0;

for (uint32_t C = 0; C < 65536; ++C) {
    float max_error = 0.0f;

    for (uint16_t x_bits : positive_normal_inputs) {
        float x = reinterpret_as_float16_or_bfloat16(x_bits);

        uint16_t y_bits = (uint16_t)(C - (x_bits >> 1));
        float y = (float)reinterpret_as_float16_or_bfloat16(y_bits);

        // Newton refinement, computed in FP32
        y = y * (1.5f - 0.5f * (x * y) * y);

        float ref = 1.0f / sqrtf(x);
        float rel_error = fabsf(y - ref) / fabsf(ref);

        max_error = max(max_error, rel_error);
    }

    if (max_error < best_error) {
        best_error = max_error;
        best_C = (uint16_t)C;
    }
}

In practice we ran this with vectorized NumPy/PyTorch code, but the search is conceptually just a minimax scan over all 16-bit constants and all positive normal inputs.

We searched two objectives: the quality of the initial bit-level approximation, and the quality after one Newton refinement. The objective was the maximum relative error over all positive normal values in the format.

Table 2. Brute-force magic constant search over all positive normal FP16 and BF16 inputs. Relative errors are measured against a float64 reference, 1.0f / sqrtf(x), evaluated on the quantized FP16/BF16 input values. The best constant for the initial approximation is not always the best constant after Newton refinement.

The BF16 result is the most visually satisfying one. The best constant is 0x5f37, which looks like the prefix of the original Quake constant:

 FP32 Quake constant:  0x5f3759df
 BF16 searched value:  0x5f37

This does not mean BF16 simply inherits the FP32 constant. It means the same bit-level structure survives in a smaller format. For BF16, the exponent layout is close enough to FP32 that the optimal 16-bit constant lands exactly on the familiar-looking prefix.

Another small lesson is that the best initial approximation is not necessarily the best constant after Newton refinement. For FP16, the initial minimax constant is 0x59bb, while the best constant after one Newton step is 0x59ba. This mirrors the older FP32 story: the magic number should be optimized for the computation we actually run, not only for the initial guess.

From One Magic Constant to a Tiny LUT

The original Quake trick uses one global constant. Every input, regardless of its mantissa, goes through the same transformation:

y_bits = C - (x_bits >> 1);

This is elegant, but it is also rigid. The approximation error is not uniform across the mantissa. Some mantissa regions are consistently overestimated, and others are consistently underestimated. So instead of asking for one magic constant to fit the entire format, we tried a tiny mantissa-conditioned correction table:

The exact exponent and mantissa widths differ between FP16 and BF16, but the LUT idea is the same: take the top few mantissa bits, map them to a bucket, and add a small correction.

// Conceptual tiny-LUT correction
uint16_t mantissa = x_bits & mantissa_mask;
uint16_t bucket = mantissa >> (mantissa_bits - bucket_bits);

uint16_t y_bits = C - (x_bits >> 1);
y_bits = y_bits + lut[bucket];

float y = reinterpret_as_float16_or_bfloat16(y_bits);

// Newton refinement in FP32
y = y * (1.5f - 0.5f * (x * y) * y);

For example, a 16-entry LUT uses 4 mantissa bits, and a 32-entry LUT uses 5 mantissa bits. These tables are tiny: for int32 corrections, 16 entries are only 64 bytes, and 32 entries are only 128 bytes. This keeps the spirit of the original hack while giving it a little local flexibility.

Table 3. Tiny-LUT correction improves Q_rsqrt-style FP16/BF16 approximations. Relative errors are measured against a float64 reference over all positive normal inputs.

The effect is small but consistent. For FP16, the 32-entry LUT reduces the maximum error from 1.7939 \times 10^{-3} to 1.3694 \times 10^{-3}. For BF16, it reduces the maximum error from 1.8992 \times 10^{-3} to 1.4965 \times 10^{-3}. The mean error is roughly halved in both formats.

In other words, the natural successor to one magic number may be a very small magic table.

Accuracy and Throughput on A100

The tiny LUT improves accuracy, but only matters if it does not cost too much on GPU. So we benchmarked the FP16 and BF16 variants on an NVIDIA A100-SXM4-40GB. Each method was warmed up for 100 runs, then timed over 10 runs.

Table 4. Accuracy and throughput trade-off on A100. Relative errors are measured against a float64 reference over all positive normal FP16 and BF16 inputs.

The main result is that the low-precision Q_rsqrt-style kernels are close to native rsqrt throughput. In FP16, magic only and magic + Newton are slightly faster than native rsqrt in this benchmark. In BF16, they are essentially tied with native rsqrt.

The 32-entry LUT adds a small overhead. For FP16, magic + Newton took 1.0632 ms, while 32-entry LUT + Newton took 1.0820 ms. For BF16, the same comparison was 1.0558 ms versus 1.0796 ms. In exchange, the mean relative error was roughly halved, and the maximum relative error dropped from 1.7939 \times 10^{-3} to 1.3694 \times 10^{-3} for FP16, and from 1.8992 \times 10^{-3} to 1.4965 \times 10^{-3} for BF16.

Native rsqrt is still far more accurate. The point is not to replace it as a high-accuracy primitive. The point is that, in low precision, a tiny table can move the Q_rsqrt-style approximation along a useful accuracy-throughput frontier with only a small extra cost.

RMSNorm-like Inputs

The previous tables evaluate all positive normal values. That is useful for worst-case behavior, but ML kernels do not usually call inverse square root on arbitrary positive numbers. In RMSNorm, the input is more structured:

r = mean(x_i ** 2) + eps
scale = 1 / sqrt(r)

So we also tested a simple RMSNorm-like distribution. we sampled activations from a normal distribution, computed r, quantized r to FP16 or BF16, and then measured the relative error of each inverse-square-root approximation.

x = sigma * torch.randn(batch, hidden_dim)
r = x.square().mean(dim=1) + eps

r_fp16 = r.to(torch.float16)
r_bf16 = r.to(torch.bfloat16)

scale = invsqrt_approx(r_fp16_or_bf16)

As the hidden dimension grows, r concentrates more tightly around its expected value, so this test is closer to the values seen inside real normalization layers than a uniform sweep over all positive floats.

Table 5. RMSNorm-like accuracy test with hidden_dim = 4096 and sigma = 1.0. Relative errors are measured against a float64 reference on quantized FP16/BF16 inputs.

The LUT still helps on this ML-shaped input distribution. For FP16, the 32-entry LUT reduces the mean relative error from 1.5320 \times 10^{-3} to 9.7337 \times 10^{-4}. For BF16, it reduces the mean relative error from 1.5466 \times 10^{-3} to 9.4188 \times 10^{-4}.

The tail also improves. The p99 error drops from 1.7333 \times 10^{-3} to 1.2731 \times 10^{-3} for FP16, and from 1.8327 \times 10^{-3} to 1.4585 \times 10^{-3} for BF16.

This is the most relevant result for ML workloads. The tiny LUT is not only improving a synthetic worst-case sweep. It also improves the error profile on the kind of values that appear inside normalization layers.

Conclusion and Limitations

The original 0x5F3759DF trick was one magic constant for scalar CPUs. On modern GPUs, that exact story no longer holds for FP32: native rsqrt is fast and far more accurate.

But in low precision, the idea still has an interesting shape. For FP16 and BF16, format-specific constants such as 0x59ba and 0x5f37 give reasonable Q_rsqrt-style approximations, and tiny mantissa-conditioned LUTs improve both mean and tail error with only a small throughput cost on A100.

This does not claim to beat native rsqrt in accuracy, nor does it prove an end-to-end Transformer speedup. The experiments are mainly A100-based and should be read as an exploratory design-space study. Still, the result suggests a nice modern twist: the successor to one magic number may be a tiny magic table.

Citation

References

1. id Software. Quake III Arena Source Code. 1999. https://github.com/id-Software/Quake-III-Arena.
2. Matthew Robertson McEniry. The Mathematics Behind the Fast Inverse Square Root Function Code. 2007. https://0x5f37642f.com/documents/McEniryMathematicsBehind.pdf.
3. Chris Lomont. Fast Inverse Square Root. 2003. https://www.lomont.org/papers/2003/InvSqrt.pdf.
4. William Kahan and K. C. Ng. sqrt implementation notes in fdlibm. 1993. https://netlib.org/fdlibm/e_sqrt.c.