平方根倒数快速计算

2021/8/14 23:09:00

本文主要是介绍平方根倒数快速计算,对大家解决编程问题具有一定的参考价值,需要的程序猿们随着小编来一起学习吧!

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×10²,数字的首位一定不为0,即不会出现0.2×10²的情况。于是,在二进制里相应的最高位一定也不会为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×2^3+0×2^2+0×2^1+1×2^0+1×2^-1+0×2^{-2}+1×2^{-3}\),其规范浮点数表达为 \(1.001101×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语言献花,向牛顿献花,向约翰·卡马克(代码编写者)献花



这篇关于平方根倒数快速计算的文章就介绍到这儿,希望我们推荐的文章对大家有所帮助,也希望大家多多支持为之网!


扫一扫关注最新编程教程