C语言之美——平方根倒数快速计算
前言
由于特殊原因，陆陆续续接触陀螺仪很长一段时间，对于各种解析算法的运算速率有了切身体会，不断追求更快、更准。最近，发现了一份比较特殊的平方根倒数速算法，一下子来了兴趣，要知道，陀螺仪解析算法里，倒数可是很常见的啊。下面来看一看这一份优美的代码，足以体现C语言的独特美感。
源码
你看不懂就对了<(￣ˇ￣)/，不然谁听我下面的哔哔呢?
float rsqrt(float number)
{long i; //32-bit numberfloat x2, y; //32-bit decimal numberconst float threehalfs = 1.5F; //1.5(also 32-bit)x2 = number * 0.5F;y = number;i = *(long *)&y; // evil floating point bit hacki = 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 removereturn 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×10²，数字的首位一定不为0，即不会出现0.2×10²的情况。于是，在二进制里相应的最高位一定也不会为0，只能为1。
研究754标准那帮人特牛，嗯~ o(￣▽￣)o，直接把确定的1省去默认，23个小数位全为真·小数位！（哈哈哈），简单来说，又省了一位确定的，将他分配给小数表示。
到现在，我们可以又下面的公式 （敲黑板(☄⊙ω⊙)☄）：
( 1 + m 2 23 ) × 2 e − 127 \left ( {1+\frac {m} {{2}^{23}}} \right )\times {2}^{e-127} (1+223m)×2e−127
其中m表示23位的小数部分，e表示8位的指数部分 由于指数有正负，所以减去127
下面我们跟着将浮点数公式取对数变换，最终得到：
l o g 2 ( 1 + m 2 23 ) + e − 127 {log}_{2}\left ( {1+\frac {m} {{2}^{23}}} \right )+e-127 log2(1+223m)+e−127
然后咋办呢~ (￣▽￣) ~*？我们这样：
l o g 2 ( 1 + x ) ≈ x {log}_{2}\left ( {1+x} \right )\approx x log2(1+x)≈x
这一波神之操作，很多人会说底数应该为e才对，不过呢，你把2带进去会发现在两端图像重合，加一个修正系数就好了
即：
l o g 2 ( 1 + x ) ≈ x + u {log}_{2}\left ( {1+x} \right )\approx x+u log2(1+x)≈x+u
（u=0.430时，在0~1的区间内，误差最小）
于是公式又能化简为：
1 2 23 ( m + 2 23 × e ) + u − 127 \frac {1} {{2}^{23}}\left ( {m+{2}^{23}\times e} \right )+u-127 2231(m+223×e)+u−127
了解这些后我们就可以做接下来的神奇操作了。
奇葩的位操作
哼哼（￣ー￣）ノ~~マ☆’.・.・:★，很显然，在754标准里，float类型是不适合位操作的，不过long类型绝对适合，为嘛呢？
如二进制数 1001.101，我们可以根据上面的表达式表达为： 1 × 2 3 + 0 × 2 2 + 0 × 2 1 + 1 × 2 0 + 1 × 2 − 1 + 0 × 2 − 2 + 1 × 2 − 3 1×2^3+0×2^2+0×2^1+1×2^0+1×2^-1+0×2^{-2}+1×2^{-3} 1×23+0×22+0×21+1×20+1×2−1+0×2−2+1×2−3，其规范浮点数表达为 1.001101 × 2 3 1.001101×2^3 1.001101×23。
由上面的等式，我们可以得出：向左移动二进制小数点一位相当于这个数除以 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 = 1 y l=\frac {1} {\sqrt {y}} l=y1
化简得 l o g ( l ) = − 1 2 l o g ( y ) log\left ( {l} \right )=-\frac {1} {2}log\left ( {y} \right ) log(l)=−21log(y)
由于我们已推导 1 2 23 ( m + 2 23 × e ) + u − 127 \frac {1} {{2}^{23}}\left ( {m+{2}^{23}\times e} \right )+u-127 2231(m+223×e)+u−127
代入化简得
( m l + 2 23 × e l ) = 3 2 × 2 23 ( 127 − u ) − 1 2 ( m y + 2 23 × e y ) \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 ) (ml+223×el)=23×223(127−u)−21(my+223×ey)
其中 3 2 × 2 23 ( 127 − u ) = 0 x 5 f 3759 d f \frac {3} {2}\times {2}^{23}\left ( {127-u} \right )=0x5f3759df 23×223(127−u)=0x5f3759df
y = *(float *)&i;
逆向操作，类型转化即可
y = y * (threehalfs - (x2 * y * y)); //1st iteration
牛顿迭代， f ( y ) = 1 y 2 − x f\left ( {y} \right )=\frac {1} {{y}^{2}}-x f(y)=y21−x，用它来推导即可
//y = y * (threehalfs - (x2 * y * y)); // 2st iteration,can be remove
迭代两次，精度更高，取决于你自己啰！！！
总结
喵的，强啊！！！
七夕节到了，向C语言献花，向牛顿献花，向约翰·卡马克（代码编写者）献花
