# 平方根倒数速算法

## 算法的切入点

${\displaystyle \|{\boldsymbol {v}}\|={\sqrt {v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}}}$ 可求得向量v的欧几里得范数，此算法正类如对欧几里得空间的两点求取其欧几里得距离
${\displaystyle {\boldsymbol {\hat {v}}}={\boldsymbol {v}}/\|{\boldsymbol {v}}\|}$ 求得的就是标准化的向量，若以${\displaystyle {\boldsymbol {x}}}$ 代表${\displaystyle v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}$ ，则有${\displaystyle {\boldsymbol {\hat {v}}}={\boldsymbol {v}}/{\sqrt {x}}}$

## 代码概览

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（对浮点数的邪恶位元hack）
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;
}


### 将浮点数转化为整数

 符号位 0 1 1 1 1 1 1 1 = 127 0 0 0 0 0 0 1 0 = 2 0 0 0 0 0 0 0 1 = 1 0 0 0 0 0 0 0 0 = 0 1 1 1 1 1 1 1 1 = −1 1 1 1 1 1 1 1 0 = −2 1 0 0 0 0 0 0 1 = −127 1 0 0 0 0 0 0 0 = −128
8位二进制补码示例

### “魔术数字”

S(ign,符号) E(xponent,指数) M(antissa,尾数)
1 位 b位 (n-1-b)
n位[文 13]

${\displaystyle m_{x}={\frac {M_{x}}{L}}}$
${\displaystyle e_{x}=E_{x}-B}$ ，其中${\displaystyle B=2^{b-1}-1}$

${\displaystyle y={\frac {1}{\sqrt {x}}}}$

${\displaystyle \log _{2}{(y)}=-{\frac {1}{2}}\log _{2}{(x)}}$
${\displaystyle \log _{2}(1+m_{y})+e_{y}=-{\frac {1}{2}}\log _{2}{(1+m_{x})}-{\frac {1}{2}}e_{x}}$

${\displaystyle m_{y}+\sigma +e_{y}=-{\frac {1}{2}}m_{x}-{\frac {1}{2}}\sigma -{\frac {1}{2}}e_{x}}$

${\displaystyle M_{y}+(E_{y}-B)L=-{\frac {3}{2}}\sigma {L}-{\frac {1}{2}}M_{x}-{\frac {1}{2}}(E_{x}-B)L}$

${\displaystyle E_{y}L+M_{y}={\frac {3}{2}}(B-\sigma )L-{\frac {1}{2}}(E_{x}L+M_{x})}$

${\displaystyle I_{y}=E_{y}L+M_{y}=R-{\frac {1}{2}}(E_{x}L+M_{x})=R-{\frac {1}{2}}I_{x}}$ ，其中${\displaystyle R={\frac {3}{2}}(B-\sigma )L}$ .

### 牛顿法提高精度

${\displaystyle y={\frac {1}{\sqrt {x}}}}$ 所求的是y的平方根倒数，以之构造以y为自变量的函数，有${\displaystyle f(y)={\frac {1}{y^{2}}}-x=0}$

${\displaystyle y_{n+1}={\frac {y_{n}(3-xy_{n}^{2})}{2}}}$ ，其中${\displaystyle f(y)={\frac {1}{y^{2}}}-x}$ ${\displaystyle f'(y)={\frac {-2}{y^{3}}}}$

## 历史与考究

id Software的创始人约翰·卡马克。这段代码虽非他所作，但常被认为与他相关。

《雷神之锤III》的代码直到QuakeCon 2005才正式放出，但早在2002年（或2003年）时，平方根倒数速算法的代码就已经出现在Usenet与其他论坛上了[1]。最初人们猜测是卡马克写下了这段代码，但他在回复询问他的邮件时否定了这个观点，并猜测可能是先前曾帮id Software优化雷神之锤的资深汇编程序员Terje Mathisen写下了这段代码；而在Mathisen的邮件裡，他表示，在1990年代初，他只曾作过类似的實作，确切来说这段代码亦非他所作。现在所知的最早實作是由Gary Tarolli在SGI Indigo中實作的，但他亦坦承他仅对常数R的取值做了一定的改进，实际上他也不是作者。在向以发明MATLAB而闻名的Cleve Moler查证后，Rys Sommefeldt则认为原始的算法是公司的Greg Walsh所发明，但他也没有任何决定性的证据能证明这一点[5]

## 注释

1. ^ 由于现代计算机系统对长整型的定义有所差异，使用长整型会降低此段代码的可移植性。具体来说，由此段浮点转换为长整型的定义可知，如若这段代码正常运行，所在系统的长整型长度应为4字节（32位），否则重新转为浮点数时可能会变成负数；而由于C99标准的广泛应用，在现今多数64位计算机系统（除使用LLP64数据模型的Windows外）中，长整型的长度都是8字节[文 6][3]
2. ^ 此处“浮点数”所指为标准化浮点数，也即有效数字部分必须满足${\displaystyle \textstyle 0\leq m_{x}<1)}$ ，可参见David Goldberg. What Every Computer Scientist Should Know About Floating-Point Arithmetic. ACM Computing Surveys. March 1991, 23 (1): 5–48. doi:10.1145/103162.103163.
3. ^ Lomont 2003确定R的方式则有所不同，他先将R分解为${\displaystyle \textstyle R_{1}}$ ${\displaystyle \textstyle R_{2}}$ ，其中${\displaystyle \textstyle R_{1}}$ ${\displaystyle \textstyle R_{2}}$ 分别代表R的有效数字域和指数域[文 7]

## 参考

### 脚注

1. Sommefeldt, Rys. Origin of Quake3's Fast InvSqrt(). Beyond3D. 2006-11-29 [2009-02-12]. （原始内容存档于2009-02-09） （英语）.
2. ^ quake3-1.32b/code/game/q_math.c. Quake III Arena. id Software. [2010-11-15]. （原始内容存档于2011-02-17）.
3. ^ Meyers, Randy. The New C: Integers in C99, Part 1. drdobbs.com. 2000-12-01 [2010-09-04]. （原始内容存档于2012-05-27）.
4. ^ Ruskin, Elan. Timing square root. Some Assembly Required. 2009-10-16 [2010-09-13]. （原始内容存档于2010-08-25）.
5. ^ Sommefeldt, Rys. Origin of Quake3's Fast InvSqrt() - Part Two. Beyond3D. 2006-12-19 [2008-04-19]. （原始内容存档于2008-05-13）.