Carmack的快速开平方根倒数算法(Fast inverse square root)

基本原理

需求\(y =\frac{1} {\sqrt{x} }\)

\(log(a^b×a^c)=bloga+cloga=(b+c)loga\)

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

-20240511163945890.webp)

字符串形式：\(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=argI(E,M)=(1+L/M​)×2^{E−B}\)

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

泰勒近似：

\(log_2​(1+m)=m+α,m∈[0,1)\), \(\alpha\)为误差

\(log_2(x)=log_2[(1+\frac{M}{L} )\times 2^{E-B}]= \frac{(E+\frac{M}{L}) \times L}{L} +\alpha -B\)

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

-20240511171551171.webp)

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

\(y=argI(E_y,M_y)=(1+\frac{M_y}{L})*2^{E_y-B}\)

\(\alpha=Error(x)=log_2(1+x)-x, x∈[0,1)\)

源代码作者采取\(α=0.0450465\)

-20240511172133421.webp)

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模块的使用，主要是对浮点数和整数进行打包和解析操作。让我解释一下各个参数的含义：

1. struct.pack('f', y)：

   • struct.pack(format, value)函数用于将值打包为指定格式的字节对象。
   • 'f'是格式化字符串，代表将值打包为浮点数。
   • y是要打包的值，即输入的浮点数。

2. struct.unpack('i', ...)：

   • struct.unpack(format, data)函数用于从字节对象中解析出值。
   • 'i'是格式化字符串，代表解析的值为整数。
   • ...是前面打包得到的字节对象，即浮点数y打包后的字节对象。

3. struct.unpack('i', struct.pack('f', y))[0]：

   • 这个表达式的含义是先将浮点数y打包为字节对象，然后再从字节对象中解析出整数值。
   • [0]表示取解析结果中的第一个值，因为struct.unpack返回的是一个元组。

综合起来，这行代码的作用是将浮点数y转换为字节对象（即二进制表示），然后再从字节对象中解析出整数值。这个操作是算法中的一个关键步骤，用于对输入的浮点数进行位操作和处理。