【转载】世界上最快的浮点开方算法
【转载】http://blog.csdn.net/wangsanhuai2010/article/details/5570730任何一个3D引擎都是通过其内部的数学模型和实现工具来展现它的力量与速度的，Quake III中使用了一个非常有意思的技巧来计算平方根倒数（inverse square root）
Carmack's 不寻常平方根倒数
第一个跳出来的便是对函数Q_rsqrt中对0x5f3759df的使用，这个数计算了一个浮点数的inverse square root，但是为什么这个函数有这样的功能呢？
观察q_math.c原本的函数：
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 );
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;
}
它不仅有效，甚至在某些CPU上，Carmack的Q_rsqrt 比(float)(1.0/sqrt(x)的计算快4倍，尽管sqrt()通常使用的是FSQRT的汇编指令！
在另一个文件code/common/cm_trace.c 中，我们发现了更简洁的对同样HACK的实现。这一次，它被用来计算一个float - sqrt(x)的平方根。注意，其中的唯一不同是在返回值上－－用返回*y取代了返回y。
float SquareRootFloat(float number) {
long i;
float x, y;
const float f = 1.5F;
x = number * 0.5F;
y= number;
i= * ( long * ) &y;
i= 0x5f3759df - ( i >> 1 );
y= * ( float * ) &i;
y= y * ( f - ( x * y * y ) );
y= y * ( f - ( x * y * y ) );
return number * y;
}
牛顿对根的近似值
上面的代码执行了众所周知的牛顿对根的近似值，像绝大多数其它迭代求近似值的计算一样，牛顿近似值假定是迭代的；每一次迭代都增强了它的准确度直至达到需要的准确度。
在牛顿近似值中的一般想法是我们我们猜测一个数x的平方根值y，我们可能通过一个简单的操作用x/y来拉平y来取得更好的猜测，使其更接近实际的平方根，例如，我们像下面这样计算2的平方根，我们假定初始的猜测是1：
2/1 = 2 ;(2 + 1) / 2 = 1.5
2/1.5 = 1.3333; ( 1.5 + 1.3333 ) / 2 = 1.4167
2/1.4167 = 1.4117;( 1.4167 + 1.4117 ) / 2 = 1.4142
And so on...
如前面所提到的，牛顿的近似值是一个大家所熟知的用以快速计算平方根的方法。但是，Carmack在初始的猜测中就选取的不寻常的值，它彻底加强了准确度并且将Quake III中计算所要的值的迭代次数降到了1次！
魔数
这个函数中真正有意思的方面是神奇的常量0x5f3759df，用来计算初始猜测的，在
i= 0x5f3759df - ( i >> 1 );
因此，把输入除以2并从神奇常量中减去。这个常数工作起来几乎是完美的－－对于一个 low relative error of 10^-3来说只要一次牛顿近似值迭代就够了。如评论中第二次迭代中展示的，这个近似值对Quake III引擎来说已经足够了。
结果，这个神奇的常数0x5f3759df成了一个迷了，在文章"Fast Inverse Square Root"，普度大学的数学家Chris Lomont研究了这个常数，用了几种精细的技术，Lomont想自己用数学方法求出这个常数来，结果令人惊奇－－Lomont用数学方法计算出来的最佳常数(0x5f37642f)有一点点不同，并且除了理论上强一些之外，它产生的结果并没有源代码中使用的原始常数好！确实，John Carmack 一定用了天才般的黑盒来找到这个常数。
只是仅仅从数字上来找的方法中，Lomont找到了一个更好的常数，这个数比原始的那个强了那么一点点。然而，实践中两个常数产生了大概相同的结果，Lomont提出这个使用了更好的常数的函数：
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f375a86- (i>>1);
x = *(float*)&i;
x = x*(1.5f-xhalf*x*x);
return x;
}
-----------------------------------------------------------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------------------------------------------------------
关于本文，东拼西凑来的，今天也已经很晚了。只是最近我在研究图形学算法的过程中，一些数字的数学的方法常常让我感到震惊，我惊叹于数学竟是如此的神奇。黑盒的应用是一门大学问，那些伟大的但是并不一定很出名的数学家常常利用黑盒另外开辟了区别于经典算法的道路，但却能带来非常人能想像的效率。在研究这些强人的作品的时候，能得到一种很美妙很利落的体会，好的算法就就像好的音乐，当欣赏他们的时候，只想拿起吉它，轻轻的弹上一首……
以其说计算机创造了21世纪，还不如说算法与数据结构创造了21世纪。前人的智慧也常让我感到无比的惊叹，而已往我却没有过这种深刻的体会，应视为较近距离接触了吧。最近在研究卡尔曼滤波器、小波变换理论、数据融合算法、以及复杂数据结构的过程中，真正体会到了什么才叫旷世绝学，虽然现在自己也只是略懂丁点儿，但这些来自前人的伟大思想，却让我振奋，让我更加明智。
最快的浮点开方算法，先看看 数学就是那么的神奇。 学习了，神奇的0x5f3759df。 很有意思。。 神奇的常数 这个神奇的数以前在本坛出现过，忘记帖子叫啥名了。。 作者也说了，这是求近似值，精度实在是有限
我试过用来做姿态解算里的归一化，相对F103没有FPU来说确实是快了一些，但是对于有DSP的F4来说不仅慢而且不精确。 maxwelllls 发表于 2015-3-1 21:03
作者也说了，这是求近似值，精度实在是有限
我试过用来做姿态解算里的归一化，相对F103没有FPU来说确实是快 ...
先记下来，等有空用我手头的f407试一下你说的对不对。 这个以前在8051上用过，确实快了很多，神奇的数学 按10进制算，这个神奇常数比最佳常数小了2640，即：
1597463007 - 1597465647 = -2640
对应相对值约为：
2640/1597463007≈1.653e-6
神奇常数小了约百万分之1.653 。
先收藏了，谢谢分享！ 请教楼主，预计效率提高多少倍？ ................ 收藏了，谢谢！ 这么神奇？赶紧记录下来。 感谢分享，看看 mark, 0x5f3759df 浮点开方运算 记下试试看 在一个项目上用过，速度快了很多 为何冠名最快，有人实测过吗？不过还是可以记下来留待以后应用。 神奇的数学，看不懂{:lol:} 数学是神奇的，这个算法收藏了，呵呵 这么神奇，做个记号。 数值分析相关的内容。 没有看懂,用的时候再细看吧 学习。。。做好记号 看不懂。。。看来得好好学习了{:titter:} 以前遊戲雷神之錘的解法 bangbangji 发表于 2015-3-1 20:30
这个神奇的数以前在本坛出现过，忘记帖子叫啥名了。。
doom启示录里面的Quake-III Arena里面使用的3D引擎使用了这个数，
我也是在本论坛上第一次知道0x5f3759df这个数的； 神奇的0x5f3759df …… 搜了一下发现一句注释：
i = 0x5f3759df - ( i >> 1 ); // what the fuck? 在AVR单片机MEGA16上试过，这个方法比最普通的SQRT函数来说，一般可以记节约75%的时间。 学习了，受教了，，， 很早论坛就有这个算法的帖子了，虽然我表示看不懂{:lol:} 谢谢分享 用过，确实很快 顶起~~~~ 抽空好好研究下算法 累啊！！ 好东西，收藏 收藏一下！！！ 好东西，看看试下 我数学不是很懂啊，唉 有字节序问题吗？ 还是牛人多啊，我先学习了！ 还是挺有意思的，有机会了在自己的程序中用一下 感谢分享！ 很流弊的样子，顶起来{:smile:} 最快的浮点开方算法，先看看 世界上最快的浮点开方算法 收藏一下 最快的浮点开方算法 很早就知道这个东西了, 雷神之锤源码开源的时候就广为流传, 用这个方法要看你的精度要求 嗯，一直使用的，精度还可以接受 谢谢，慢慢学习下 最精简的1/sqrt()函数：
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating VALUE
i = 0x5f375a86- (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits BACK to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
return x;
}
===============
mark 好东西，验证一下。 厉害！有空试试。 maxwelllls 发表于 2015-3-1 21:03
作者也说了，这是求近似值，精度实在是有限
我试过用来做姿态解算里的归一化，相对F103没有FPU来说确实是快 ...
f4的dsp库也是用这个数字来开方的 不明觉厉 xwy19920211 发表于 2015-3-4 20:03
f4的dsp库也是用这个数字来开方的
不是用这个数，用sqrt算的是很精确的数，跟卡西欧的计算器算出来的一样 数学永远都是这么神奇。。。。 谢谢分享 数学成绩差的表示看不懂啊 有木有0x5f3759df的数学原理啊，表示这样看只知道有用但是终究不是自己的东西 mack,不知道0x5f3759df这个数是怎么发现的。 不错！谢谢分享！！
卡马克大神的故事,恩 这里有关于这个常数的说明
http://zh.wikipedia.org/zh/平方根倒数速算法 平方根倒数速算法,好东西！ 数学，在程序算法中涉及很多。 数学学得好
最精简的1/sqrt()函数：
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating VALUE
i = 0x5f375a86- (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits BACK to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
return x;
}
这个不错，以用在交流电压采集真有效值计算程序中 maxwelllls 发表于 2015-3-4 20:36
不是用这个数，用sqrt算的是很精确的数，跟卡西欧的计算器算出来的一样
很有可能会用到，数学库也是用迭代算法算的，这个方法可以快速得到一个近似值，收敛更快 感謝分享! 真有在MCU 上用的,这个速度快是肯定的,不过 个精度不能通用,原来的用途是3d和图形处理,要求低的多 {:curse:}第一次听说这个，好神奇的说 用在单精度的，好像用于双精度的也有类似的一个数值 好神奇的数字！ 不懂原理，但觉得很神奇，又学到了东西了，有时间也试一下，看看效果，多谢！ 参考一下 谢谢分享 很神奇，想不到的快。
for(i=1;i<200000;i++)
{
for(j=1;j<10000;j++)
f=Q_rsqrt( j ) ;
if(f>1000000)printf("x=%d=>%f\n",j,f);
}
printf ("AAA \n");
for(i=1;i<200000;i++)
{
for(j=1;j<10000;j++)
f=1.0/sqrt((float)j);
if(f>1000000)printf("x=%d=>%f\n",j,f);
}
printf ("BBB \n");
getchar();
以上程序，在PC上试过，我的是I3（4130），时间差不多，第一个14秒，第二个16秒左右的样子。 有的算法境界高的另众生无法膜拜 不懂，未用过这个神数。 mark，最快的浮点开方算法，这种好东西收起来备用 感谢分享…… 以前只知道那个牛顿迭代方法开跟号，基本上看RP的那种，不知道这个优化过的是不是也要看RP。。 mark，好东西收起来，日后有用{:lol:} 居然还可以这样开方的,卡大英明 以后也许能用到，谢谢 学习了，标记下也许以后会用到。 神奇的数字，数硬合体更神奇！ 太神奇了，谢谢啊！
太神奇了，谢谢啊！ Thank you !!!
页:
[1]