【学习笔记】快速平方根倒数算法

来源：互联网发布：vs for mac 离线安装编辑：程序博客网时间：2024/05/16 18:48

很以前久，看《DOOM启示录》的时候就看到，当年卡马克大神在《雷神之锤》中使用了一个神奇的数字，能够通过位操作快速计算平方根倒数y=1x√。但是当时并没有深究，这两天偶然看到了这篇文章，终于将我这个多年的“未解之谜”解开了，特此做下笔记，图片取自原文。

代码实现

首先是这个神奇算法的代码：

float fast_inverse_sqrt(float x){    float half_x = 0.5 * x;    int i = *((int *)&x); // 以整数方式读取X    i = 0x5f3759df - (i>>1); // 神奇的步骤    x = *((float *)&i); // 再以浮点方式读取i    x = x*(1.5 - (half_x * x * x)); // 牛顿迭代一遍提高精度    return x;}

其中最最重要的就是第二步，通过这个神奇的步骤，能够得到x平方根倒数的近似值。前后穿插着用整型读取浮点数等等，实在是不知所云。但实际上，其中包含着很有意思的数学背景（每次一提到数学，我就觉得好方）。

背后的数学

首先来重新温习一下机组里面的浮点数基础知识（机组60分默默飘过）。

单精度浮点数

单精度浮点数有32bit，包括1位符号位s，8位指数位e，23位尾数位m（赞一下CSDN的markdown编辑器的latex数学公式的功能~)：

s e⋯e8 m⋯m23

其表示的数值为：(1+m′)2e′，这里的m’代表尾数部分的数值，m’只包含小数位，即m′=0.m⋯m23，e’则是指数部分的数值，只有整数位，e′=e⋯e8−B。为了表示负的指数，指数部分要减去一个偏移量B。
再规定M表示以整数方式解释尾数二进制位的值，E表示相同的方式解释指数部分的值。可以得到以下的转换关系：

m′=MLL=223
e′=E−BB=127

那么对于一个浮点数的二进制位s e⋯e8 m⋯m23以整数方式解读的话，其值就为：I=EL+M（因为平方根输入只能为正，所以默认s为0）。

平方根倒数近似计算

下面开始进入正题：
对于输入x，我们要计算的是y=1x√，将该式子两边取对数：log2(y)=−12log2(x)。然后我们将x，y都替换为浮点数表示形式：

log2((1+m′y)2e′y)=−12log2((1+m′x)2e′x)
log2(1+m′y)+e′y=−12(log2(1+m′x)+e′x)

观察发现，上式等号两边均有log2(1+v)形式的式子，v的范围是0到1，而该式在这个定义域内的函数图像非常接近一条直线：

这里写图片描述

因此我们近似的用x+σ替换得到：

m′y+σ+e′y=−12(m′x+σ+e′x)

σ为预先计算的一个常数。接下来将数值转为整形方式读取二进制位的形式：

MyL+σ+Ey−B=−12(MxL+σ+Ex−B)

移项合并后得到：

32L(σ−B)+My+LEy=−12(Mx+LEx)

恰好我们又发现，在等式的两边都含有以整数方式解释浮点数的表示I=M+LE，进一步化简有：

Iy=−12Ix+κκ=32L(B−σ)

由此，我们得到了一个比较有用的近似公式，κ便是神奇的数字0x5f3759df 。该公式表明计算x的平方根倒数时，我们只需要用整型方式读取x，再代入上式计算得到结果的整形解释值，然后再用浮点数方式读取即可！我当时就震惊了。

牛顿迭代提高精度

经过上一步，我们有了初步的近似结果，为了进一步提高精度，我们再用牛顿迭代法计算一次。牛顿迭代法描述的是，给定连续函数f(x)，对于方程f(x)=0，确定任意初始的x0，我们可以通过下面这个迭代式计算得到其解：

xn+1=xn−f(xn)f′(xn)

我们要计算的式子为y=1x0√，x0为已知（这个x0是指输入，不是牛顿迭代的初始值），将其转换为关于y的方程就是f(y)=1y2−x0=0，根据牛顿迭代可得：

yn+1=yn−1y2n−x0−2y3n
yn+1=32yn−12x0y3n

也就是最后一步的代码了。