[C++知识库]C语言之美——平方根倒数快速计算

C语言之美——平方根倒数快速计算

前言

由于特殊原因，陆陆续续接触陀螺仪很长一段时间，对于各种解析算法的运算速率有了切身体会，不断追求更快、更准。最近，发现了一份比较特殊的平方根倒数速算法，一下子来了兴趣，要知道，陀螺仪解析算法里，倒数可是很常见的啊。下面来看一看这一份优美的代码，足以体现C语言的独特美感。

源码

你看不懂就对了<(￣ˇ￣)/，不然谁听我下面的哔哔呢?

float rsqrt(float number)
{
    long i;                        //32-bit number
    float x2, y;                   //32-bit decimal number
    const float threehalfs = 1.5F; //1.5(also 32-bit)

    x2 = number * 0.5F;
    y = number;
    i = *(long *)&y;           // evil floating point bit 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)); // 2st iteration,can be remove

    return y;
}

IEEE754标准

定点数

定点数，自然是相对于浮点数而言的。
何为定点？何为浮点？

定点，通俗来说就是小数点位置是固定的。
官方一丢丢的解释(>▽<)：

在定点数表达法中，其小数点固定地位于实数所有数字中间的某个位置。例如，货币的表达就可以采用这种表达方式，如 55.00 或者 00.55 可以用于表达具有 4 位精度，小数点后有两位的货币值。由于小数点位置固定，所以可以直接用 4 位数值来表达相应的数值。

但我们不难发现，定点数表达法的缺点就在于其形式过于僵硬，固定的小数点位置决定了固定位数的整数部分和小数部分，不利于同时表达特别大的数或者特别小的数。因此，最终绝大多数现代的计算机系统都采纳了所谓的浮点数表达法。

浮点数

浮点，通俗来说就是小数点位置是浮动的。
官方一点的解释(o゜▽゜)o☆：

浮点数表达法采用了科学计数法来表达实数，即用一个有效数字。一个基数（Base）、一个指数（Exponent）以及一个表示正负的符号来表达实数。浮点数利用指数达到了浮动小数点的效果，从而可以灵活地表达更大范围的实数。

表示方法 （仅限于单精度浮点数）

在IEEE 754标准中浮点数由三部分组成：符号位（sign bit），有偏指数（biased exponent），小数（fraction）。浮点数分为两种，单精度浮点数（single precision）和双精度浮点数（double precision），它们两个所占的位数不同。

单精度浮点数（共32位）：
1个符号位
8个指数位
23个小数位

先来看看下面的例子，简单阐述了浮点数的表示方法：

在简单知道了浮点数表示方法之后，一个有意思的事情又开始了：
我们知道了，浮点数的表示方法其实灵感来源于科学计数法，科学计数法表示的数如1.2×102，数字的首位一定不为0，即不会出现0.2×102的情况。于是，在二进制里相应的最高位一定也不会为0，只能为1。
研究754标准那帮人特牛，嗯~ o(￣▽￣)o，直接把确定的1省去默认，23个小数位全为真·小数位！（哈哈哈），简单来说，又省了一位确定的，将他分配给小数表示。

到现在，我们可以又下面的公式 （敲黑板(?⊙ω⊙)?）：
$$\left(1+\frac{m}{2^{23}}\right)\times 2^{e?127}$$
其中m表示23位的小数部分，e表示8位的指数部分 由于指数有正负，所以减去127

下面我们跟着将浮点数公式取对数变换，最终得到：
$${log}_2\left(1+\frac{m}{2^{23}}\right)+e?127$$

然后咋办呢~ (￣▽￣) ~*？我们这样：
$${log}_2\left(1+x\right)\approx x$$
这一波神之操作，很多人会说底数应该为e才对，不过呢，你把2带进去会发现在两端图像重合，加一个修正系数就好了
即：

$${log}_2\left(1+x\right)\approx x+u$$
（u=0.430时，在0~1的区间内，误差最小）

于是公式又能化简为：

$$\frac{1}{2^{23}}\left(m+2^{23}\times e\right)+u?127$$

了解这些后我们就可以做接下来的神奇操作了。

奇葩的位操作

哼哼（￣ー￣）ノ~~マ☆’.?.?:★，很显然，在754标准里，float类型是不适合位操作的，不过long类型绝对适合，为嘛呢？

如二进制数 1001.101，我们可以根据上面的表达式表达为：$1\times 2^3+0\times 2^2+0\times 2^1+1\times 2^0+1\times 2^{?}1+0\times 2^{?2}+1\times 2^{?3}$，其规范浮点数表达为 $1.001101\times 2^3$。
由上面的等式，我们可以得出：向左移动二进制小数点一位相当于这个数除以 2，而向右移动二进制小数点一位相当于这个数乘以 2。

所以知道了为毛long类型适合位操作了吧。

C语言的技巧

强制转换

不懂强制转换的同学叉出去

 i = (long)y;

将浮点数(单双精度)转换为整数时，将舍弃浮点数的小数部分， 只保留整数部分。

嗯，不知叫啥的方法

 i = *(long *)&y;

就只是改变了类型，内存里的数值没有改变
（简单来说，C语言：“你又骗我！这到底是不是float？”）

牛顿迭代法

详细操作自己去看看百度，爷懒了不想写<(ˉ^ˉ)>
致敬牛顿

解析源码

 i = *(long *)&y;           // evil floating point bit hack

很显然，这是进行类型转化

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

what the fuck? 为啥有一个0x5f3759df，嘛呢？
我们看下面的操作：
先设 $$l=\frac{1}{\sqrt{y}}$$
化简得 $$log\left(l\right)=?\frac{1}{2}log\left(y\right)$$
由于我们已推导 $$\frac{1}{2^{23}}\left(m+2^{23}\times e\right)+u?127$$
代入化简得
$$\left(m_l+2^{23}\times e_l\right)=\frac{3}{2}\times 2^{23}\left(127?u\right)?\frac{1}{2}\left(m_y+2^{23}\times e_y\right)$$
其中 $$\frac{3}{2}\times 2^{23}\left(127?u\right)=0x5f3759df$$

y = *(float *)&i;

逆向操作，类型转化即可

 y = y * (threehalfs - (x2 * y * y)); //1st iteration

牛顿迭代，$$f\left(y\right)=\frac{1}{y^2}?x$$，用它来推导即可

//y = y * (threehalfs - (x2 * y * y)); // 2st iteration,can be remove

迭代两次，精度更高，取决于你自己啰！！！

总结

喵的，强啊！！！
七夕节到了，向C语言献花，向牛顿献花，向约翰·卡马克（代码编写者）献花