快速开平方根倒数算法(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^{23}-1]$
参数相互转化	$e=E-B,m=\frac{M}{L}$	$E=e+B,M=m\times L$
有了x如何求参数e,E	$e=J_e(x)=科学计数法后半部分2的指数e$	$E=I_E(x)=e+B$
有了x如何求参数m, M	$m=J_m(x)=科学计数法前半部分的小数m$	$M=I_M(x)=\frac{m}{L}$
有了参数e,m(E, M)如何反求x	$x=argJ(e,m)=(1+m)\times 2^e$	$x=argI(E,M)=(1+\frac{M}{L})\times 2^{E-B}$

3. 同一储存单元32bits的两种不同意义

以上的表格中，我们详细的对比了对于同一个浮点数x，1.人类或是2.计算机是如何理解、储存及还原它的。
然而，不仅x的理解方式有两种，同样有两种理解方式的还有一个32 bits的数字。当我们在计算机的存储单元中发现了一个32位的二进制数，在没有被告知它的类型是1.整形或是2.浮点型时，我们是并不能判断次32 bits的数字到底代表什么意义的。

下面的表格就对比了同一个32 bits的数字的两种不同释义。注意，我们可以把这个32 bits数字划分为三部分：S是这个32bits的从左起第0位，E是第1-第8位，M是第9-第31位。E, M的意义完全跟上面相同，同样的还有B,L，我们都延续第2部分的概念。
$$S_0{E}_1{E}_2...E_7{E}_8{M}_9{M}_{10}...M_{30}{M}_{31}$$

任意一个32 bits的二进制数	整数	浮点数
t是如何生成的	非科学计数法的一个整数化为二进制数，如果不够32位就在前面添加0	根据IEEE 754将浮点数存成SEM表示方式
知道t的类型如何还原	$x=Int(E,M)=(E+\frac{M}{L})\times L$	$x=Float(E,M)=(1+\frac{M}{L})\times 2^{E-B}$

注意到， 1. 这里的$Float=argI$，即二者的函数(function，也可理解为功能)是一样的。2. 注意整数的$Int(E,M)$和浮点数$Float(E,M)$的形式都很相似，都是$(a+\frac{M}{L})\times b$的形式。

有了上面注意点里面的1、2两点内容，我们就会思考，这是巧合么？我们又能利用这两点巧合做到什么么？ 那么接下来的重头戏即将呈现！！！不过在此之前，我们停下来回顾一下上面的准备工作都准备了哪些：

1. 一个任意的浮点数x都可以被Encoding成为两种形式。一种是我们人脑作为encoder记住其e, m；另一种是计算机作为encoder把此浮点数按SEM方法存储成一个32bits的二进制数存入存储单元。 注意，这两种Encoding的形式都是一一对应(双射)的！！！
2. 一个任意的位于存储单元中的32bits 二进制数字t，我们都可以将其Decoding为两个不同的数字。一种方法是告诉我们x是整数，那我们就要按照其为整数进decoding；另一种是告诉我们t是浮点数，那我们就按照其为浮点数进行decoding。同样，这两种Decoding的方法也都是一一对应(双射)的！！！

4. 公式推导

掌握以上知识，我们试图着手将输入的自变量x从乘方级运算(开平方根就是乘$-\frac{1}{2}$次方)通过映射变为$I(x)$后变为加减乘除级别运算。此处的x在计算机是以SEM方法存储，因此通过回顾第2部分的表格 $$x=argI(E,M)=(1+\frac{M}{L})\times 2^{E-B}\text{\hspace{1em}(5)}$$
以2为底取对数$lo{g}_2$
$$lo{g}_2(x)=lo{g}_2[(1+\frac{M}{L})\times 2^{E-B}]$$ $$=(E-B)+lo{g}_2(1+\frac{M}{L})$$ $$=[E+lo{g}_2(1+\frac{M}{L})]-B\text{\hspace{1em}(6)}$$
OK，稍稍稍微又卡住了，又推不下去了。那么我们需要一点近似，如果把$\frac{M}{L}$当做一个整体$m$, m就是刚才第二部分的m，其意义是二进制科学计数法的小数部分, $m\in [0,1)$。那么$F_1(m)=lo{g}_2(1+m)$与$F_2(m)=m$两条直线在[0, 1]上及其相似：

1. $F_1(0)=F_2(0)=0$且$F_1(1)=F_2(1)=1$，意味着两条直线在[0, 1]两个端点重合
2. 两条直线的最大差 $$Max(F_1(m)-F_2(m))=0.086071,m\in [0,1)$$
3. 两条直线的平均差 $${\int }_0^1{F}_1(m)-F_2(m)dx/(1-0)=0.057304959,m\in [0,1)$$

注意， 上述第2条、第3条可以看我另一篇博客：数学公式之求 log2(1+x)-x的积分
因此，由于$F_1(m)=lo{g}_2(1+m)$与$F_2(m)=m$两条直线在[0, 1]上的相似性，我们可以用$m$来做$lo{g}_2(1+m)$的近似，只不过需要加一个误差变量$\alpha$, 即$$lo{g}_2(1+m)=m+\alpha ,m\in [0,1)\text{\hspace{1em}(7)}$$
既然有了(7), 将其带入(6)，我们又可以继续前进了！！！
$$log(x)=[E+lo{g}_2(1+\frac{M}{L})]-B$$ $$=[E+\frac{M}{L}+\alpha ]-B$$ $$=\frac{(E+\frac{M}{L})\times L}{L}+\alpha -B\text{\hspace{1em}(8)}$$
我的天！我们看到了什么啊！一个似曾见过的数字$(E+\frac{M}{L})\times L$。 这个东西好像在哪里见过，赶紧去回顾一下上面第3部分的表格，发现它的意义是32bits 二进制数decoding为整数。

4. 本文核心

有了(8)，我们便得到了本文的核心。考虑到其重要性，我单开了一个文章标题来介绍它
$$log(小数x)=(小数先按SEM方法存储为32bits，再按照整数方法还原)\times \frac{1}{L}+\alpha -B$$
即
$$log(x)=Float(I_E(x),I_M(x))\times \frac{1}{L}+\alpha -B\text{\hspace{1em}(9)}$$
那么(9)式就是我们本文的核心，如果能看懂(9)，就说明你明白了本文在想干什么了。注意：$L,B$都是常量，而$\alpha$是一个变量；$I_E,I_M$在第2部分表格有介绍，$Float(E,M)$在第3部分表格有介绍。

5. 开平方根倒数

关键性的技术突破了，那接下来就很容易了。对$x$进行inverse square root $$lo{g}_2(x^{-1/2})=-\frac{1}{2}\times lo{g}_2[(1+\frac{M}{L})\times 2^{E-B}]$$ $$=-\frac{1}{2}(\alpha -B)-\frac{1}{2}\times \frac{Float(I_E(x),I_M(x))}{L}\text{\hspace{1em}(10)}$$

什么？ 对小数开平方根的倒数 取对数后，只需要把其化为IEEE-754标准encoding存入存储单元，再按照整数标准decoding把此32bits 二进制数还原为整数，那么对这个整数只需要进行加减乘除法的运算就相当于对原数进行乘方开方的运算。

换句话说，我们彻底完成了需求分析中的，乘方到加减乘除的降维！！！
现在重新来看我们的问题(1)式
$$y=x^{-1/2}$$ 两边取对数 $$logy=-\frac{1}{2}logx$$ 把(8)带入上式
$$\frac{(E_y+\frac{M_y}{L})\times L}{L}+{\alpha }_y-B=-\frac{1}{2}\times (\frac{(E_x+\frac{M_x}{L})\times L}{L}+{\alpha }_x-B)$$ 整理移项得 $$M_y+E_{y}L=[\frac{3}{2}BL-({\alpha }_a+\frac{{\alpha }_b}{2})L]-\frac{1}{2}(M_x+E_{x}L)\text{\hspace{1em}(11)}$$
如果把$M+EL$当做一个整体，它是一个function $I$；且把(11)式[]方括号中心的全部内容看作R的话，那么(11)式可以变成什么？
$$I(y)=R-\frac{1}{2}I(x)$$
回头一看，此式正是(3)式！！！对于此$I$，我们可再熟悉不过了，它就是第2部分表中的encoder $I$。再帮大家回忆一下：$$E_x=I_E(x)=e+B$$ $$M_x=I_M(x)=\frac{m}{L}$$这样，有了x，我们就可以通过$I$求出为三十二位bits二进制数的$I(x)$。在对$I(x)$进行乘以$-\frac{1}{2}$(或者是-$I(x)$>>1)的基础操作后，我们给其加一个常数R就能得到$I(y)$。之后，我们需要利用同样是第2部分表中的decoder
$$y=argI(E_y,M_y)=(1+\frac{M_y}{L})\times 2^{E_y-B}$$ 就可以将三十二位bits的二进制数还原为一个浮点数。而这个浮点数，恰恰就是我们要求的$y$，也就是$x^{-1/2}$。

由此，我们算是讲完了Fast inverse square root算法中最核心的思想，这是一个从开平方到加减乘除质的飞跃！！！

剩下的一共还有两项工作。

1. 是求
   $$R=[\frac{3}{2}BL-({\alpha }_a+\frac{{\alpha }_b}{2})L]\text{\hspace{1em}(12)}$$
   注意到上式$B,L$都是确定的常数，而${\alpha }_a,{\alpha }_b$是变量，因此我们要做的是试图用一个确定的常量 来近似所有的$\alpha$(也包括${\alpha }_a,{\alpha }_b$)，使得(12)大致成立！！！
2. 由于1.的近似R，因此求出的$y$只是$x^{-1/2}$的近似解，我们需要用牛顿迭代加深这个初始解$y$的精确度。（根据源代码，1次牛顿迭代就够了。换句话说，用了FISR算法以后只需要一次牛顿迭代，这个数值解就已经近似于解析解了。当然你可以做n次都可以，可以但没必要好吧）

6. 工作扫尾

针对遗留下来的2个小问题，第1个问题在我的另一篇文章两条曲线相似度的探究(MSE推广) 有所提及。(博主为写这一篇文章容易嘛，不仅花了一个礼拜绞尽脑汁理解算法，还写了两篇辅助工作的博客。那能看到这里的人都是人才了，楼主深表欣慰，也希望你们点个赞评个论呀~)
由于 $$\alpha =Error(x)=lo{g}_2(1+x)-x,x\in [0,1)$$
因此这里有两种计算$\alpha$的方法。
1.是取Error最大值得一半，$\alpha =\frac{Max(Error(x))}{2}=0.0430357$
2.是取Error函数在[0,1]上面积分的平均，$\alpha =\frac{{\int }_0^{1}Error(x)dx}{1-0}=0.0573049$

然而，这些都只是猜测与理论。经过实验，最终源代码的作者最终采用了$$\alpha =0.0450465\text{\hspace{1em}(13)}$$这个数值（我也不知道怎么来的，应该是按照性能和准确率做实验做出来的）。
那么，将(13)结果带入(12):
$$R=[\frac{3}{2}BL-(\alpha +\frac{\alpha }{2})L]$$ $$=1.5\times L\times (B-1.5\alpha )$$ $$=1.5\times 2^{23}\times (127-1.5\times 0.0450465)$$ $$\approx 1,597,463,007$$ $$=0x5f3759df\text{\hspace{1em}(14)}$$
至此，我们就解决了遗留工作1. $R=0x5f3759df$，这也就是代码

	i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 

中的那个magic number。
懂了么？明白这个magic number怎么来的了么？！

至于第2个遗留工作就不是本文的讨论内容了，那是牛顿老人家的研究成果，详情自行google或百度即可。

三、总结

额，太累了，还没想好，下次再说吧。就五一快乐！(8点终于可以回家了)