Carmack的快速开平方根倒数算法

基本原理

需求 $y = \frac{1}{\sqrt{x}}$

$l o g (a^{b} \times a^{c}) = b l o g a + c l o g a = (b + c) l o g a$

32位浮点表示法：二进制的科学计数法
符号位1+阶码8（有符号的反码表示幂指数）+小数位23（二进制小数首位必为1，默认，只需表示小数位即可）

字符串形式： $S_{0} E_{1} E_{2} . . . E_{7} E_{8} M_{9} M_{10} . . . M_{30} M_{31} $

正数符号位 $S_{0}$ 忽略，并引入两个常量 $B = 127, L = 2^{23}$ ，关联到指数E和小数M

故：
$x = a r g I (E, M) = (1 + L / M ) \times 2^{E - B}$

对x取对数：
$l o g_{2} (x) = l o g_{2} [(1 + \frac{M}{L}) \times 2^{E - B}]$

泰勒近似：

$l o g_{2} (1 + m) = m + α, m \in [0, 1)$ , $α$ 为误差

$l o g_{2} (x) = l o g_{2} [(1 + \frac{M}{L}) \times 2^{E - B}] = \frac{(E + \frac{M}{L}) \times L}{L} + α - B$

而 $(E + \frac{M}{L}) \times L$ 正是32bits 二进制数decoding为整数

$I (y) = R - \frac{1}{2} I (x)$

$y = a r g I (E_{y}, M_{y}) = (1 + \frac{M_{y}}{L}) * 2^{E_{y} - B}$

$α = E r r o r (x) = l o g_{2} (1 + x) - x, x \in [0, 1)$

源代码作者采取 $α = 0.0450465$

C代码实现

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

python实现

import struct
import math
  
def fast_inverse_square_root(number):
    threehalfs = 1.5
    x2 = number * 0.5
    y = number
    i = struct.unpack('i', struct.pack('f', y))[0]
    i = 0x5f3759df - (i >> 1)
    y = struct.unpack('f', struct.pack('i', i))[0]
    y = y * (threehalfs - (x2 * y * y))
    y  = y * ( threehalfs - ( x2 * y * y ) );   # 2nd iteration, this can be removed
    return y

# 测试
number = 101.0
inverse_sqrt = fast_inverse_square_root(number)
error = inverse_sqrt / (1/math.sqrt(number))
print("Fast Inverse Square Root of", number, "is:", inverse_sqrt, "error", error)

在这段代码中，struct.unpack('i', struct.pack('f', y))[0]这行代码涉及到了struct模块的使用，主要是对浮点数和整数进行打包和解析操作。让我解释一下各个参数的含义：

struct.pack('f', y)：
- struct.pack(format, value)函数用于将值打包为指定格式的字节对象。
- 'f'是格式化字符串，代表将值打包为浮点数。
- y是要打包的值，即输入的浮点数。
struct.unpack('i', ...)：
- struct.unpack(format, data)函数用于从字节对象中解析出值。
- 'i'是格式化字符串，代表解析的值为整数。
- ...是前面打包得到的字节对象，即浮点数y打包后的字节对象。
struct.unpack('i', struct.pack('f', y))[0]：
- 这个表达式的含义是先将浮点数y打包为字节对象，然后再从字节对象中解析出整数值。
- [0]表示取解析结果中的第一个值，因为struct.unpack返回的是一个元组。