【转载】世界上最快的浮点开方算法

JJKwong

【转载】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原本的函数：

1. float Q_rsqrt( float number )   
   
2. {  
   
3.   long i;   
   
4.   float x2, y;   
   
5.   const float threehalfs = 1.5F;   
   
6.   x2 = number * 0.5F;   
   
7.   y  = number;  
   
8.   i  = * ( long * ) &y;  // evil floating point bit level hacking   
   
9.   i  = 0x5f3759df - ( i >> 1 );  
   
10.   y  = * ( float * ) &i;  
    
11.   y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration   
    
12.   y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed   
    
13.   return y;   
    
14. }

复制代码

   

它不仅有效，甚至在某些CPU上，Carmack的Q_rsqrt 比(float)(1.0/sqrt(x)的计算快4倍，尽管sqrt()通常使用的是FSQRT的汇编指令！
在另一个文件code/common/cm_trace.c 中，我们发现了更简洁的对同样HACK的实现。这一次，它被用来计算一个float - sqrt(x)的平方根。注意，其中的唯一不同是在返回值上－－用返回*y取代了返回y。  
　　

1. float SquareRootFloat(float number) {  
   
2.     long i;   
   
3.     float x, y;   
   
4.     const float f = 1.5F;   
   
5.     x = number * 0.5F;  
   
6.     y  = number;   
   
7.     i  = * ( long * ) &y;   
   
8.     i  = 0x5f3759df - ( i >> 1 );   
   
9.     y  = * ( float * ) &i;   
   
10.     y  = y * ( f - ( x * y * y ) );   
    
11.     y  = y * ( f - ( x * y * y ) );   
    
12.     return number * y;   
    
13. }

复制代码

牛顿对根的近似值
上面的代码执行了众所周知的牛顿对根的近似值，像绝大多数其它迭代求近似值的计算一样，牛顿近似值假定是迭代的；每一次迭代都增强了它的准确度直至达到需要的准确度。
在牛顿近似值中的一般想法是我们我们猜测一个数x的平方根值y，我们可能通过一个简单的操作用x/y来拉平y来取得更好的猜测，使其更接近实际的平方根，例如，我们像下面这样计算2的平方根，我们假定初始的猜测是1：
　　

1. 2/1 = 2 ;  (2 + 1) / 2 = 1.5   
   
2. 2/1.5 = 1.3333; ( 1.5 + 1.3333 ) / 2 = 1.4167   
   
3. 2/1.4167 = 1.4117;  ( 1.4167 + 1.4117 ) / 2 = 1.4142   
   
4. And so on...   

复制代码

如前面所提到的，牛顿的近似值是一个大家所熟知的用以快速计算平方根的方法。但是，Carmack在初始的猜测中就选取的不寻常的值，它彻底加强了准确度并且将Quake III中计算所要的值的迭代次数降到了1次！

魔数
　　这个函数中真正有意思的方面是神奇的常量0x5f3759df，用来计算初始猜测的，在
　　

1. 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提出这个使用了更好的常数的函数：

1. float InvSqrt(float x)   
   
2. {   
   
3.   float xhalf = 0.5f*x;   
   
4.   int i = *(int*)&x;   
   
5.   i = 0x5f375a86- (i>>1);   
   
6.   x = *(float*)&i;   
   
7.   x = x*(1.5f-xhalf*x*x);   
   
8.   return x;   
   
9. }

复制代码

   

-----------------------------------------------------------------------------------------------------------------------------------------------

-----------------------------------------------------------------------------------------------------------------------------------------------

-----------------------------------------------------------------------------------------------------------------------------------------------

　　关于本文，东拼西凑来的，今天也已经很晚了。只是最近我在研究图形学算法的过程中，一些数字的数学的方法常常让我感到震惊，我惊叹于数学竟是如此的神奇。黑盒的应用是一门大学问，那些伟大的但是并不一定很出名的数学家常常利用黑盒另外开辟了区别于经典算法的道路，但却能带来非常人能想像的效率。在研究这些强人的作品的时候，能得到一种很美妙很利落的体会，好的算法就就像好的音乐，当欣赏他们的时候，只想拿起吉它，轻轻的弹上一首……

　　以其说计算机创造了21世纪，还不如说算法与数据结构创造了21世纪。前人的智慧也常让我感到无比的惊叹，而已往我却没有过这种深刻的体会，应视为较近距离接触了吧。最近在研究卡尔曼滤波器、小波变换理论、数据融合算法、以及复杂数据结构的过程中，真正体会到了什么才叫旷世绝学，虽然现在自己也只是略懂丁点儿，但这些来自前人的伟大思想，却让我振奋，让我更加明智。

dingwei

最快的浮点开方算法，先看看

lihq97

数学就是那么的神奇。

登云钓月

学习了，神奇的0x5f3759df。

jie5313

很有意思。。

huangqi412

神奇的常数

bangbangji

这个神奇的数以前在本坛出现过，忘记帖子叫啥名了。。

maxwelllls

作者也说了，这是求近似值，精度实在是有限
我试过用来做姿态解算里的归一化，相对F103没有FPU来说确实是快了一些，但是对于有DSP的F4来说不仅慢而且不精确。

dawanpi

maxwelllls 发表于 2015-3-1 21:03
作者也说了，这是求近似值，精度实在是有限
我试过用来做姿态解算里的归一化，相对F103没有FPU来说确实是快 ...

先记下来，等有空用我手头的f407试一下你说的对不对。

eedesign

这个以前在8051上用过，确实快了很多，神奇的数学

wiser803

按10进制算，这个神奇常数比最佳常数小了2640，即：
    1597463007 - 1597465647 = -2640
对应相对值约为：
    2640/1597463007≈1.653e-6
神奇常数小了约百万分之1.653 。

yzh

先收藏了，谢谢分享！

styleno1

请教楼主，预计效率提高多少倍？

hanshiruo

................

zhangxun0712

收藏了，谢谢！

Excellence

这么神奇？赶紧记录下来。

cjt5132

感谢分享，看看

Xplain

mark, 0x5f3759df 浮点开方运算

cwei

记下试试看

aasf3

在一个项目上用过，速度快了很多

freefei

为何冠名最快，有人实测过吗？不过还是可以记下来留待以后应用。

jiangkehong

神奇的数学，看不懂

syauxwm

数学是神奇的，这个算法收藏了，呵呵

就在今天

这么神奇，做个记号。

dhbighead

数值分析相关的内容。

vibrate

没有看懂,用的时候再细看吧

1148729990

学习。。。做好记号

永恒520

看不懂。。。看来得好好学习了

fw190d9

以前遊戲雷神之錘的解法

jm2011

bangbangji 发表于 2015-3-1 20:30
这个神奇的数以前在本坛出现过，忘记帖子叫啥名了。。

doom启示录里面的Quake-III Arena里面使用的3D引擎使用了这个数，
我也是在本论坛上第一次知道0x5f3759df这个数的；

kakarotto

神奇的0x5f3759df ……

kakarotto

搜了一下发现一句注释：

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

huangxiaolpbany

在AVR单片机MEGA16上试过，这个方法比最普通的SQRT函数来说，一般可以记节约75%的时间。

captain_steve

学习了，受教了，，，

kickdown

很早论坛就有这个算法的帖子了，虽然我表示看不懂

wsm80828

谢谢分享

fm0826

用过，确实很快

ZL_electric

顶起~~~~

cdlxzlp

抽空  好好研究下  算法 累啊！！

hackdj

好东西，收藏

kneken

收藏一下！！！

xiaomao1490

好东西，看看试下

xiatao1800

我数学不是很懂啊，唉

abutter

有字节序问题吗？

nsj21n

还是牛人多啊，我先学习了！

hawgcq

还是挺有意思的，有机会了在自己的程序中用一下

whm_WUHAN

感谢分享！

DM_2014

很流弊的样子，顶起来

LK9286

最快的浮点开方算法，先看看

xzp21st

世界上最快的浮点开方算法

pop_wu

收藏一下

孤独飞行

最快的浮点开方算法

windrarara

很早就知道这个东西了, 雷神之锤源码开源的时候就广为流传, 用这个方法要看你的精度要求

cyh

嗯，一直使用的，精度还可以接受

pp99sky

谢谢，慢慢学习下

licheng0620

最精简的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

alexcnsz

好东西，验证一下。

wolflsh

厉害！有空试试。

xwy19920211

maxwelllls 发表于 2015-3-1 21:03
作者也说了，这是求近似值，精度实在是有限
我试过用来做姿态解算里的归一化，相对F103没有FPU来说确实是快 ...

f4的dsp库也是用这个数字来开方的

sensitiveplant

不明觉厉

maxwelllls

xwy19920211 发表于 2015-3-4 20:03
f4的dsp库也是用这个数字来开方的

不是用这个数，用sqrt算的是很精确的数，跟卡西欧的计算器算出来的一样

leijiayou

数学永远都是这么神奇。。。。   谢谢分享

kevin_me

数学成绩差的表示看不懂啊

迗眞de瘋子

有木有0x5f3759df的数学原理啊，表示这样看只知道有用但是终究不是自己的东西

yoyotansa

mack,不知道0x5f3759df这个数是怎么发现的。

ddcour

不错！谢谢分享！！

windingway

卡马克大神的故事,恩

max

这里有关于这个常数的说明
http://zh.wikipedia.org/zh/平方根倒数速算法

xiaoao360

平方根倒数速算法,好东西！

crazydtone

数学，在程序算法中涉及很多。

飞剑

数学学得好

Serial

最精简的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;
}  

这个不错，以用在交流电压采集真有效值计算程序中

myxiaonia

maxwelllls 发表于 2015-3-4 20:36
不是用这个数，用sqrt算的是很精确的数，跟卡西欧的计算器算出来的一样

很有可能会用到，数学库也是用迭代算法算的，这个方法可以快速得到一个近似值，收敛更快

bossmaster_chan

感謝分享!

mcu5i51

真有在MCU 上用的,这个速度快是肯定的,不过 个精度不能通用,原来的用途是3d和图形处理,要求低的多

windancerhxw

第一次听说这个，好神奇的说

Achin

用在单精度的，好像用于双精度的也有类似的一个数值

2avr

好神奇的数字！

最笨的企鹅

不懂原理，但觉得很神奇，又学到了东西了，有时间也试一下，看看效果，多谢！

timom333

参考一下 谢谢分享

sdsdc

很神奇，想不到的快。

sdsdc

        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秒左右的样子。

lzlym182

有的算法境界高的另众生无法膜拜

zlutian

不懂，未用过这个神数。

youmu

mark，最快的浮点开方算法，这种好东西收起来备用

what007

感谢分享……

bli19

以前只知道那个牛顿迭代方法开跟号，基本上看RP的那种，不知道这个优化过的是不是也要看RP。。

qq854149876

mark，好东西收起来，日后有用

kevin_ares

居然还可以这样开方的,卡大英明

my二月兰

以后也许能用到，谢谢

hricchaya

学习了，标记下也许以后会用到。

chouxier99318

神奇的数字，数硬合体更神奇！

lionboa

太神奇了，谢谢啊！

gloryglory

太神奇了，谢谢啊！

liang16888

Thank you !!!