搜索
bottom↓
回复: 121

【转】神一样的算法,一个数开平方并取倒,0x5f3759df

  [复制链接]

出0入0汤圆

发表于 2012-1-22 22:47:34 | 显示全部楼层 |阅读模式
http://www.douban.com/note/93460299/

Quake-III Arena (雷神之锤3)是90年代的经典游戏之一。该系列的游戏不但画面和内容不错,而且即使计算机配置低,也能极其流畅地运行。这要归功于它3D引擎的开发者约翰-卡马克(John Carmack)。事实上早在90年代初DOS时代,只要能在PC上搞个小动画都能让人惊叹一番的时候,John Carmack就推出了石破天惊的Castle Wolfstein, 然后再接再励,doom, doomII, Quake...每次都把3-D技术推到极致。他的3D引擎代码资极度高效,几乎是在压榨PC机的每条运算指令。当初MS的Direct3D也得听取他的意见,修改了不少API。

   最近,QUAKE的开发商ID SOFTWARE 遵守GPL协议,公开了QUAKE-III的原代码,让世人有幸目睹Carmack传奇的3D引擎的原码。
   这是QUAKE-III原代码的下载地址:
  http://www.fileshack.com/file.x?fid=7547

(下面是官方的下载网址,搜索 “quake3-1.32b-source.zip” 可以找到一大堆中文网页的
ftp://ftp.idsoftware.com/idstuff/source/quake3-1.32b-source.zip)

   我们知道,越底层的函数,调用越频繁。3D引擎归根到底还是数学_运算。那么找到最底层的数学_运算函数(在game/code/q_math.c), 必然是精心编写的。里面有很多有趣的函数,很多都令人惊奇,估计我们几年时间都学不完。

在game/code/q_math.c里发现了这样一段代码。它的作用是将一个数开平方并取倒,经测试这段代码比(float)(1.0/sqrt(x))快4倍:
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 ); // what the fuck?
   y = * ( float * ) &i;
   y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
   // y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

   #ifndef Q3_VM
   #ifdef __linux__
     assert( !isnan(y) ); // bk010122 - FPE?
   #endif
   #endif
   return y;
}

    函数返回1/sqrt(x),这个函数在图像处理中比sqrt(x)更有用。
    注意到这个函数只用了一次叠代!(其实就是根本没用叠代,直接运算)。编译,实验,这个函数不仅工作的很好,而且比标准的sqrt()函数快4倍!要知道,编译器自带的函数,可是经过严格仔细的汇编优化的啊!
  
   这个简洁的函数,最核心,也是最让人费解的,就是标注了“what the fuck?”的一句
      i = 0x5f3759df - ( i >> 1 );

再加上y = y * ( threehalfs - ( x2 * y * y ) );
两句话就完成了开方运算!而且注意到,核心那句是定点移位运算,速度极快!特别在很多没有乘法指令的RISC结构CPU上,这样做是极其高效的。

算法的原理其实不复杂,就是牛顿迭代法,用x-f(x)/f'(x)来不断的逼近f(x)=a的根。

简单来说比如求平方根,f(x)=x^2=a ,f'(x)= 2*x,f(x)/f'(x)=x/2,把f(x)代入

x-f(x)/f'(x)后有(x+a/x)/2,现在我们选a=5,选一个猜测值比如2,
那么我们可以这么算
5/2 = 2.5; (2.5+2)/2 = 2.25; 5/2.25 = xxx; (2.25+xxx)/2 = xxxx ...
这样反复迭代下去,结果必定收敛于sqrt(5),没错,一般的求平方根都是这么算的
但是卡马克(quake3作者)真正牛B的地方是他选择了一个神秘的常数0x5f3759df 来计算那个猜测值
就是我们加注释的那一行,那一行算出的值非常接近1/sqrt(n),这样我们只需要2次牛 顿迭代就可以达到我们所需要的精度.
好吧 如果这个还不算NB,接着看:


普渡大学的数学家Chris Lomont看了以后觉得有趣,决定要研究一下卡马克弄出来的
这个猜测值有什么奥秘。Lomont也是个牛人,在精心研究之后从理论上也推导出一个
最佳猜测值,和卡马克的数字非常接近, 0x5f37642f。卡马克真牛,他是外星人吗?

传奇并没有在这里结束。Lomont计算出结果以后非常满意,于是拿自己计算出的起始
值和卡马克的神秘数字做比赛,看看谁的数字能够更快更精确的求得平方根。结果是
卡马克赢了... 谁也不知道卡马克是怎么找到这个数字的。

最后Lomont怒了,采用暴力方法一个数字一个数字试过来,终于找到一个比卡马克数
字要好上那么一丁点的数字,虽然实际上这两个数字所产生的结果非常近似,这个暴
力得出的数字是0x5f375a86。

Lomont为此写下一篇论文,"Fast Inverse Square Root"。



论文下载地址:
http://www.math.purdue.edu/~clomont/Math/Papers/2003/InvSqrt.pdf
http://www.matrix67.com/data/InvSqrt.pdf
点击此处下载 ourdev_714253IQ3KPV.pdf(文件大小:148K) (原文件名:InvSqrt.pdf)



参考:<IEEE Standard 754 for Binary Floating-Point Arithmetic><FAST INVERSE SQUARE ROOT>


最后,给出最精简的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;
}
大家可以尝试在PC机、51、AVR、430、ARM、上面编译并实验,惊讶一下它的工作效率。


==============增加================
看matrix67大侠的解释
http://www.matrix67.com/blog/archives/362

这样的代码速度肯定飞快,我就不用多说了;但算法的原理是什么呢?其实说穿了也不是很神,程序首先猜测了一个接近1/sqrt(number)的值,然后两次使用牛顿迭代法进行迭代。根号a的倒数实际上就是方程1/x^2 - a = 0的一个正实根,它的导数是-2/x^3。运用牛顿迭代公式x' = x - f(x)/f'(x),式子化简为x' = x * (1.5 - 0.5a * x^2)。迭代几次后,x的值将趋于1/sqrt(a)。
    但这段代码真正牛B的是那个神秘的0x5f3759df,因为0x5f3759df - (i >> 1)出人意料地接近根号y的倒数。人们都不知道这个神秘的常数是怎么来的,只能把它当作神来膜拜。这个富有传奇色彩的常数到底咋回事,很少有人说得清楚。这篇论文(就是上面的附件)比较科学地解释了这个常数。

阿莫论坛20周年了!感谢大家的支持与爱护!!

月入3000的是反美的。收入3万是亲美的。收入30万是移民美国的。收入300万是取得绿卡后回国,教唆那些3000来反美的!

出0入0汤圆

发表于 2012-1-22 22:56:06 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-1-22 23:00:24 | 显示全部楼层
真不错。学习了

出0入0汤圆

发表于 2012-1-22 23:06:16 | 显示全部楼层
大过年的发帖不容易,给你个大红包。

出0入0汤圆

发表于 2012-1-22 23:08:01 | 显示全部楼层
mark先

出0入25汤圆

发表于 2012-1-22 23:08:30 | 显示全部楼层
这才是程序猿!国内有多少程序猿能踏实下心来反复优化自己的代码?这就是国内和国外技术的差距,基础差距

出0入0汤圆

发表于 2012-1-22 23:36:11 | 显示全部楼层
非常好!!顶!

出15入190汤圆

发表于 2012-1-22 23:38:30 | 显示全部楼层
多谢楼主带着大家领略奇人奇事

出0入42汤圆

发表于 2012-1-22 23:44:39 | 显示全部楼层
全部代码

quake3 源代码ourdev_714257X7RZZA.zip(文件大小:5.46M) (原文件名:quake3-1.32b-source.zip)

出0入0汤圆

发表于 2012-1-22 23:45:25 | 显示全部楼层
火星人,欢迎来到地球。

出0入0汤圆

发表于 2012-1-22 23:58:20 | 显示全部楼层
我靠,我还以为就我一人在此溜达,没想到还有这么多仁兄在!祝大家新年快乐哈!

出0入0汤圆

发表于 2012-1-23 00:03:54 | 显示全部楼层
精彩!!

出0入0汤圆

发表于 2012-1-23 00:18:19 | 显示全部楼层
i&#160;=&#160;0x5f3759df&#160;-&#160;(&#160;i&#160;>>&#160;1&#160;);&#160;//&#160;what&#160;the&#160;fuck?

出0入0汤圆

发表于 2012-1-23 02:09:13 | 显示全部楼层
有料!!


------------------------------------------------------------------------------------------
回复【13楼】techonly  
well, it's the big difference between engineers from china and us.
i think some engineers in us can be extrmely excel at their fields that we can not only regard them as engineers, but also as scientists/geniuses.
-----------------------------------------------------------------------
不必妄自菲薄,只要把心放在自己所做的东西上,精亦求精,並对自己所做的东西像是自己小孩一羕感到自豪就是了.

出0入0汤圆

发表于 2012-1-23 07:54:48 | 显示全部楼层
学习一下

出0入0汤圆

发表于 2012-1-23 08:58:40 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-1-23 09:00:48 | 显示全部楼层
quake3 里面的暴力枪还是很还念的

出0入0汤圆

发表于 2012-1-23 09:16:54 | 显示全部楼层
mark,牛X!新年快乐!

出0入0汤圆

发表于 2012-1-23 09:36:06 | 显示全部楼层
niuX  ~~~

出0入0汤圆

发表于 2012-1-23 10:18:52 | 显示全部楼层
欢迎lz从火星来地球过新年

出0入0汤圆

发表于 2012-1-23 10:35:38 | 显示全部楼层
mark

出675入8汤圆

发表于 2012-1-23 10:38:03 | 显示全部楼层
Mark

出0入0汤圆

发表于 2012-1-23 10:47:54 | 显示全部楼层
红包红包

出0入0汤圆

发表于 2012-1-23 13:03:06 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-1-23 15:20:22 | 显示全部楼层
不错,领教了。

出0入0汤圆

发表于 2012-1-23 17:22:10 | 显示全部楼层
真牛!!

出0入0汤圆

发表于 2012-1-23 18:47:51 | 显示全部楼层
厉害!!!

出330入1862汤圆

发表于 2012-1-23 19:28:39 | 显示全部楼层
mark!

出0入0汤圆

发表于 2012-1-23 19:33:29 | 显示全部楼层
论坛有个里06年的帖子。。。。楼主火星啦。

出0入0汤圆

发表于 2012-1-23 20:04:06 | 显示全部楼层
强啊,mark

出0入0汤圆

发表于 2012-1-23 20:26:36 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-1-23 23:25:10 | 显示全部楼层
w

出5入42汤圆

发表于 2012-1-23 23:41:32 | 显示全部楼层
强帖留名

出0入0汤圆

发表于 2012-1-23 23:49:06 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-1-24 09:35:50 | 显示全部楼层
回复【14楼】mylogin
-----------------------------------------------------------------------

这跟膜拜天才是两回事的。

出0入0汤圆

发表于 2012-1-24 10:31:58 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-1-24 10:37:26 | 显示全部楼层
mark。见到过。

出0入0汤圆

发表于 2012-1-24 10:49:26 | 显示全部楼层
快速 平方根分之1 算法源码 MARK

出65入0汤圆

发表于 2012-1-24 10:53:53 | 显示全部楼层
what the fuck?

出0入0汤圆

发表于 2012-1-24 11:17:50 | 显示全部楼层
看看

出0入0汤圆

发表于 2012-1-24 18:08:43 | 显示全部楼层
MARK,膜拜下

出0入0汤圆

发表于 2012-1-24 20:21:28 | 显示全部楼层
bu dong

出0入0汤圆

发表于 2012-1-24 22:10:37 | 显示全部楼层
mark

出0入17汤圆

发表于 2012-1-24 22:40:06 | 显示全部楼层
回复【13楼】techonly
well, it's the big difference between engineers from china and us.
i think some engineers in us can be extrmely excel at their fields that we can not only regard them as engineers, but also as scientists/geniuses.
-----------------------------------------------------------------------

They are artists indeed.

出0入0汤圆

发表于 2012-1-24 23:04:25 | 显示全部楼层
刚刚在8051上试了一下,速度比用math库快3倍,但是精度不怎么样,只能精确到3位有效数字左右。用库函数则完全一致(对比于windows自带的计算器)

出0入0汤圆

发表于 2012-1-24 23:15:33 | 显示全部楼层
mark 平方根算法

出0入0汤圆

发表于 2012-1-25 00:46:40 | 显示全部楼层
mark…真心觉得作者功底在那里

出0入0汤圆

发表于 2012-1-25 00:51:28 | 显示全部楼层
mark先

出10入12汤圆

发表于 2012-1-25 09:37:11 | 显示全部楼层
不错

出0入0汤圆

发表于 2012-1-25 10:10:26 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-1-25 10:28:32 | 显示全部楼层
MMMMMMMMMMARK

出0入0汤圆

发表于 2012-1-25 10:44:25 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-1-25 11:14:07 | 显示全部楼层
回复【52楼】theloong
mmmmmmmmmmark
-----------------------------------------------------------------------
键盘卡住了?

出0入0汤圆

发表于 2012-1-25 12:02:40 | 显示全部楼层
算法学习。

出0入0汤圆

发表于 2012-1-25 12:08:48 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-1-26 00:05:45 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-1-26 03:07:38 | 显示全部楼层
mark

出0入0汤圆

发表于 2012-8-1 16:21:39 | 显示全部楼层
神奇的算法

出0入0汤圆

发表于 2012-8-1 16:45:33 | 显示全部楼层
i = 0x5f3759df - ( i >> 1 ); // what the fuck?

出0入0汤圆

发表于 2012-8-1 17:03:23 | 显示全部楼层
顶一个     

出0入0汤圆

发表于 2012-8-8 11:58:13 | 显示全部楼层
怎么这么久都没人来看一看啊!顶一个

出0入34汤圆

发表于 2012-8-11 11:17:12 | 显示全部楼层
长见识了,佩服!~  佩服!~

出0入4汤圆

发表于 2012-8-11 11:46:56 | 显示全部楼层
好复杂的算法啊。

出0入0汤圆

发表于 2012-8-15 17:28:34 | 显示全部楼层
what the fuck?
神奇

出35入0汤圆

发表于 2012-8-15 18:07:28 | 显示全部楼层
咱们国内的人都是中国人,很少有火星来的.

出0入0汤圆

发表于 2012-8-15 18:22:13 | 显示全部楼层
真神奇的常数,这个要标记一下。

出0入0汤圆

发表于 2012-8-19 23:20:21 | 显示全部楼层
哈哈,看完,我笑了,真实“what the fuck”,不是一般的人能搞出来的

出0入0汤圆

发表于 2012-8-19 23:38:14 | 显示全部楼层
这么好的东西怎么没几个人回啊。
顶一个。

出0入0汤圆

发表于 2012-8-20 17:29:28 | 显示全部楼层
测了一下,玩STM32的可以无视了。
#include "stm32f10x_lib.h"
#include "math.h"
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;
}
int main(void)
{
  float AA=28;
  float BB;
  BB=InvSqrt(AA);
  BB=1/BB;
  AA=28;
  BB=1/sqrtf(AA);
  BB=1/BB;
}
这就是全部程序。STM32软仿,8mhz主频,用这个函数需要29.12us,用math里的sqrtf并求倒数需要22.12us。
math里sqrtf是float 函数,sqrt是double函数,调用float并求倒数的话需要41us。
不知道这和STM32自带硬件乘、除法器是否有关

出0入0汤圆

发表于 2012-8-20 17:34:25 | 显示全部楼层
这个真是厉害啊,感谢楼主带了这么精彩的东西~

出0入0汤圆

发表于 2012-8-20 21:59:19 | 显示全部楼层
数学家就是不一样啊。如此如此之节约。

出0入0汤圆

发表于 2012-8-20 22:56:00 | 显示全部楼层
咋挖出来的啊 牛贴不沉啊,MARK

出0入0汤圆

发表于 2012-8-28 16:45:01 | 显示全部楼层

出0入0汤圆

发表于 2012-8-28 17:37:46 | 显示全部楼层
不错!收藏

出0入0汤圆

发表于 2012-9-5 16:27:53 | 显示全部楼层
这个要MARK

出0入0汤圆

发表于 2012-9-12 16:29:37 | 显示全部楼层
太强,膜拜

出0入0汤圆

发表于 2012-9-23 21:13:13 | 显示全部楼层
好,学习了。

出0入0汤圆

发表于 2012-10-13 22:35:15 | 显示全部楼层
再留一次爪印。

出0入0汤圆

发表于 2012-10-21 22:27:35 | 显示全部楼层
太牛B了。

出0入0汤圆

发表于 2012-10-21 23:35:19 | 显示全部楼层
神一样的算法

出0入0汤圆

发表于 2012-10-22 00:17:21 来自手机 | 显示全部楼层
算法让人头疼的东西,不过楼主我学习到了。

出0入0汤圆

发表于 2012-10-25 23:45:12 | 显示全部楼层
mark!!!!!!!!!

出0入0汤圆

发表于 2012-10-25 23:51:33 | 显示全部楼层
正在膜拜中

出0入0汤圆

发表于 2012-10-31 10:44:06 | 显示全部楼层
what the fuck?

出0入0汤圆

发表于 2012-10-31 11:44:01 | 显示全部楼层
强啊,mark

出0入0汤圆

发表于 2012-10-31 12:05:48 | 显示全部楼层
mark
!!!!!!!!

出0入0汤圆

发表于 2012-10-31 12:24:16 | 显示全部楼层
牛啊!!

出0入0汤圆

发表于 2012-10-31 12:44:08 | 显示全部楼层
谁能把里面有用的函数都摘出来

出0入0汤圆

发表于 2012-10-31 13:04:14 | 显示全部楼层
神神神神神神神神神神神神

出0入0汤圆

发表于 2012-10-31 13:42:09 | 显示全部楼层
mark,牛B

出0入0汤圆

发表于 2012-10-31 13:55:38 | 显示全部楼层
收藏+MARK

出0入0汤圆

发表于 2012-10-31 16:26:27 | 显示全部楼层
what the fuck?

出0入0汤圆

发表于 2012-10-31 16:54:00 | 显示全部楼层
看了热血沸腾啊

出100入0汤圆

发表于 2012-10-31 17:11:14 | 显示全部楼层
mark 有意思,慢慢读~

出0入0汤圆

发表于 2012-10-31 17:28:10 | 显示全部楼层
厉害!!!

出350入8汤圆

发表于 2012-10-31 17:55:37 | 显示全部楼层
学习学习!

出20入0汤圆

发表于 2012-10-31 19:28:19 | 显示全部楼层
牛X  膜拜啊!

出0入0汤圆

发表于 2012-10-31 20:02:31 | 显示全部楼层
mark~~~~~~~~~

出0入0汤圆

发表于 2013-10-1 11:28:42 | 显示全部楼层
牛逼的算法以后再算开发再也不用sqrt了

出0入0汤圆

发表于 2013-10-1 16:26:44 | 显示全部楼层
数学是伟大的
回帖提示: 反政府言论将被立即封锁ID 在按“提交”前,请自问一下:我这样表达会给举报吗,会给自己惹麻烦吗? 另外:尽量不要使用Mark、顶等没有意义的回复。不得大量使用大字体和彩色字。【本论坛不允许直接上传手机拍摄图片,浪费大家下载带宽和论坛服务器空间,请压缩后(图片小于1兆)才上传。压缩方法可以在微信里面发给自己(不要勾选“原图),然后下载,就能得到压缩后的图片】。另外,手机版只能上传图片,要上传附件需要切换到电脑版(不需要使用电脑,手机上切换到电脑版就行,页面底部)。
您需要登录后才可以回帖 登录 | 注册

本版积分规则

手机版|Archiver|amobbs.com 阿莫电子技术论坛 ( 粤ICP备2022115958号, 版权所有:东莞阿莫电子贸易商行 创办于2004年 (公安交互式论坛备案:44190002001997 ) )

GMT+8, 2024-4-26 02:56

© Since 2004 www.amobbs.com, 原www.ourdev.cn, 原www.ouravr.com

快速回复 返回顶部 返回列表