快速开平方根倒数算法(Fast inverse square root)的一点探究

一、写在前面

1. 提示

请朋友们先行阅读Fast inverse square root algorithm(下称FISR)算法。本文内容适用于已经了解此算法基本原理的朋友们，还不了解的读者请移步以下参考资料（推荐程度按照以下顺序）:

(1) 中文且强推: 魔力数字与快速平方根倒数算法
(2) Wikipedia官方资料(需要翻墙) : Wikipedia Fast inverse square root
(3) 中文且做补充资料: 数学之美：平方根倒数速算法中的神奇数字 0x5f3759df
(4) 英文paper原文: InvSqrt.pdf(CHRIS LOMONT)

2. 背景与前情

大学非常厉害的舍友面试阿里的后端时，面试官让他不用函数库，实现一个快速开平方根的算法。当舍友把这个问题抛给我的时候，我满脑子想的都是：二分法！
显然这是大家都能想得到的算法，其时间复杂度为$O(logN)$，更准确点，其实需要做的操作是$lo{g}_2(\frac{X}{error})$次，X为被开平方的数，error是预设的可允许误差。

But！ 据说面试官听到舍友的二分法后，略带失望的摇了摇头(有演绎的成分)，说：再想想。经过5分钟的煎熬，最后面试官妥协了：“那用二分法做吧，其实这个题有线性时间解法”。

这个线性时间的解法一下子就让我充满了兴趣。开平方根，怎么可能有线性时间的算法？？？第一个反应是Impossible，第二个反应是竭尽脑汁思考除了二分法以外的其他方法。在经过长达10分钟的思考后，我终于…
明白了一句话：终日而思矣，不如须臾之所学也！好吧，上网搜得了！

最终的最终，经过一番学习以后，我现在脑子里对此题有了三种解法。

• 二分法(正常学过算法的都应该秒想到的)
• 牛顿迭代法(数值计算学过，应该属于正常人能想得到的非常好用而且有效的方法)
• Fast inverse sqaure root后再取倒数，也就是本文的核心

接下来，我试图从自己理解的角度，阐述一下对于FISR的理解。

二、正文

先把FISR代码贴在这儿，注意注释为what the fuck的那一行：

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

1. 需求分析

问题: 求$y=\frac{1}{\sqrt{x}}$
对此问题，我们变变形，方便以下计算 $$y=\frac{1}{\sqrt{x}}=x^{-1/2}=x^0\times x^{-1/2}\text{\hspace{1em}(1)}$$

已有工具: 因为不让用函数库，因此我们有的就是$+,-,\times ,÷$四种基础工具。
联想: 如果能把平方这种级别的运算“降维”变成$+,-$就好了，怎么做呢？
$⟶$: 数学上最常用的手段就是取对数$log$，比如$log(a^b\times a^c)=bloga+cloga=(b+c)loga$。
那么，如果我们对(1)式取对数，就降维成了
$$logy=log{x}^0-\frac{1}{2}logx\text{\hspace{1em}(2)}$$ 那么，开平方根倒数的计算就变成了$-\frac{1}{2}$的计算，那不就大功告成了！！！
事实上，观察(2)，它跟以下的方程很像 $$I(y)=R-\frac{1}{2}I(x)\text{\hspace{1em}(3)}$$ 或者是那一行代码非常之像
$$i=0x5f3759df-(i>>1)\text{\hspace{1em}(4)}$$ 注意(4)中等是左边的$i$其实就是$I(y)$, 神奇数字$0x5f3759df$其实就是$R$，最后$-(i>>1)$其实就是$-\frac{1}{2}i$。

是不是感觉有点思路了？有了上面非常相似的(2)(3)式，我们现在关心的有两点：1.是求$I$和$argI$，因为我们需要把$x$变为$I(x)$以及把$I(y)$还原为$y$。2.是求R所对应的神奇数字。

2. 必备工具之IEEE-754浮点数表示方法

相应的Wikipedia的介绍和例子放在下面(英文不好的同学们去别的地方搜一下吧，这部分的基础原理不在本文介绍之中)。

注意到，我们被开平方根的数字是正数，所以符号位S=0在本文中永远成立，接下来不再讨论。此外，下文中简称IEEE-754的32 bits表述一个浮点数的方法为SEM法(S符号位：一位二进制数，E指数位：八位二进制数，M尾数位：二十三位二进制数)。

下面的表格对比了：对于同一个浮点数正数x的不同表示方法，方便接下来讨论。注意，没有理解此表格的朋友请理解好了再往下走，切勿囫囵吞枣。

为了以下方便表示，我们引入两个常量
$$B=127,L=2^{23}$$ 至于原因，是因为IEEE 32位浮点数表示中，E有8位，所以$B=2^8-1=127$， 而M有23位，所以$L=2^{23}$。

输入x	二进制科学计数法	SEM 32bits法
应用场景	实际人类将x化为二进制时	根据IEEE754 计算机储存x
x具体表示为	$1.m\times 2^e$ ,($m=0.b_1{b}_2...b_n$)	SEM ,($S=0$, $E:8bits$ $M:23bits$)
参数取值范围	$e是整数,m\in [0,1)$	$E\in [0,2^8-1],M\in [0,2^{}$