Fast inverse square root in Python
Approximates 1/√x by reinterpreting the float bits as an integer and applying a Newton step.
- File
fast_inverse_square_root.py, 60 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/fast-inverse-square-root/. The authors are listed in prov/10WFYW.jsonl.
"""Fast inverse square root (the Quake III routine), in emulated binary32."""
import math
import struct
MAGIC = 0x5F3759DF
_F32 = struct.Struct("<f")
_U32 = struct.Struct("<I")
SMALLEST = 2.0**-125 # x and x/2 are normal binary32 numbers
LARGEST = (2.0 - 2.0**-23) * 2.0**127 # the largest finite binary32 number
def f32(v: float) -> float:
"""Round a binary64 value to the nearest binary32 value (ties to even)."""
return _F32.unpack(_F32.pack(v))[0]
def bits32(x: float) -> int:
"""The 32-bit pattern of the binary32 value x, as an unsigned integer."""
return _U32.unpack(_F32.pack(x))[0]
def float32(i: int) -> float:
"""The binary32 value whose pattern is the unsigned 32-bit integer i."""
return _F32.unpack(_U32.pack(i & 0xFFFFFFFF))[0]
def fast_inverse_square_root(x: float, magic: int = MAGIC, newton_steps: int = 1) -> float:
"""Approximate 1/sqrt(x) with the bit trick and Newton steps in binary32.
x is first rounded to binary32. Precondition: the rounded x lies in
[2^(-125), largest binary32]; otherwise ValueError is raised. Each
operation is a correctly rounded binary32 operation (a binary64 operation
on binary32 inputs, rounded to binary32, gives the same result). With the
default magic constant and one step, the relative error y*sqrt(x) - 1 lies
in [-0.0017525, 2e-7]. Cost: O(newton_steps) operations, O(1) space.
"""
x = f32(x)
if not (SMALLEST <= x <= LARGEST):
raise ValueError("x must be a binary32 value in [2^(-125), max]")
x2 = f32(0.5 * x) # @line L1
i = bits32(x) # @line L2
i = magic - (i >> 1) # @line L3
y = float32(i) # @line L4
for _ in range(newton_steps):
y = f32(y * f32(1.5 - f32(f32(x2 * y) * y))) # @line L5
return y # @line L6
def initial_guess(x: float, magic: int = MAGIC) -> float:
"""The guess of line L4 alone (no Newton step), for a binary32 x in the domain."""
return fast_inverse_square_root(x, magic, newton_steps=0)
def relative_error(x: float, y: float) -> float:
"""y*sqrt(x) - 1, computed in binary64 (its own error is about 1e-16)."""
return y * math.sqrt(x) - 1.0
test_fast_inverse_square_root.py 105 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/fast-inverse-square-root/. The authors are listed in prov/10WFYW.jsonl.
import math
import random
from fractions import Fraction
import pytest
from fast_inverse_square_root import (
LARGEST,
SMALLEST,
bits32,
f32,
fast_inverse_square_root,
float32,
initial_guess,
relative_error,
)
LOW, HIGH = -0.0017525, 2e-7 # the proved bound of proofs/correctness.md
LOW0, HIGH0 = -0.0343758, 0.0339603 # the proved bound on the initial guess
def test_example_table() -> None:
rows = [
(0.15625, 0x3E200000, 0x402759DF, 2.6148602962493896, 2.5254862308502197),
(1.0, 0x3F800000, 0x3F7759DF, 0.9662150740623474, 0.9983071684837341),
(2.0, 0x40000000, 0x3F3759DF, 0.7162150740623474, 0.7069300413131714),
(3.7297205924987793, 0x406EB3BE, 0x3F000000, 0.5, 0.5168924331665039),
(4.0, 0x40800000, 0x3EF759DF, 0.4831075370311737, 0.49915358424186707),
(100.0, 0x42C80000, 0x3DD359DF, 0.10319875925779343, 0.09984488040208817),
]
for x, bx, g, y0, y1 in rows:
assert bits32(x) == bx
assert 0x5F3759DF - (bx >> 1) == g
assert initial_guess(x) == y0 == float32(g)
assert fast_inverse_square_root(x) == y1
def test_worst_case_input() -> None:
x = 2.0 * (1 + Fraction(7254974, 2**23))
xf = float(x)
assert f32(xf) == xf
assert initial_guess(xf) == 0.5
e0 = relative_error(xf, 0.5)
assert abs(e0 - (-0.0343758)) < 1e-7
e1 = relative_error(xf, fast_inverse_square_root(xf))
assert LOW <= e1 <= -0.0017522
def test_scaling_by_four() -> None:
rng = random.Random(7)
for _ in range(300):
x = f32(rng.uniform(1.0, 4.0))
g = initial_guess(x)
for k in range(-60, 60):
xs = x * 4.0**k
assert initial_guess(xs) == g * 2.0**-k
def test_domain_edges() -> None:
for x in (SMALLEST, LARGEST, 2.0**-124, 2.0**127, 1.0, f32(math.pi)):
e0 = relative_error(x, initial_guess(x))
e1 = relative_error(x, fast_inverse_square_root(x))
assert LOW0 <= e0 <= HIGH0
assert LOW <= e1 <= HIGH
for bad in (0.0, -0.0, -1.0, math.inf, -math.inf, math.nan, 2.0**-126, 1e-40, 2.0**-149):
with pytest.raises(ValueError):
fast_inverse_square_root(bad)
def test_random_against_oracle() -> None:
rng = random.Random(20240501)
worst = 0.0
for _ in range(20000):
x = f32(2.0 ** rng.uniform(-124.9, 127.9))
if not (SMALLEST <= x <= LARGEST):
continue
y = fast_inverse_square_root(x)
e1 = relative_error(x, y)
assert LOW <= e1 <= HIGH
e0 = relative_error(x, initial_guess(x))
assert LOW0 <= e0 <= HIGH0
oracle = 1.0 / math.sqrt(x)
assert abs(y - oracle) <= 0.0017525 * oracle
worst = min(worst, e1)
assert worst < -0.0015 # the bound is nearly reached on random inputs
def test_dense_slice_near_the_worst_case() -> None:
start = bits32(3.7297205924987793) - 3000
for i in range(start, start + 6000):
x = float32(i)
assert LOW0 <= relative_error(x, initial_guess(x)) <= HIGH0
assert LOW <= relative_error(x, fast_inverse_square_root(x)) <= HIGH
def test_two_steps_and_other_constants() -> None:
rng = random.Random(99)
for _ in range(2000):
x = f32(rng.uniform(1.0, 4.0))
assert abs(relative_error(x, fast_inverse_square_root(x, newton_steps=2))) < 5e-6
assert abs(relative_error(x, fast_inverse_square_root(x, 0x5F375A86))) < 0.0017513
assert abs(relative_error(x, initial_guess(x, 0x5F37642F))) < 0.034213
What the next level needs
- the tests did not run in this build