转自:http://www.cnblogs.com/pkuoliver/archive/2010/10/06/1844725.html

源码下载地址:http://diducoder.com/sotry-about-sqrt.html

好吧,我承认我标题党了,不过既然你来了,就认真看下去吧,保证你有收获。

我们平时经常会有一些数据运算的操作,需要调用sqrt,exp,abs等函数,那么时候你有没有想过:这个些函数系统是如何实现的?就拿最常用的sqrt函数来说吧,系统怎么来实现这个经常调用的函数呢?

虽然有可能你平时没有想过这个问题,不过正所谓是“临阵磨枪,不快也光”,你“眉头一皱,计上心来”,这个不是太简单了嘛,用二分的方法,在一个区间中,每次拿中间数的平方来试验,如果大了,就再试左区间的中间数;如果小了,就再拿右区间的中间数来试。比如求sqrt(16)的结果,你先试(0+16)/2=8,8*8=64,64比16大,然后就向左移,试(0+8)/2=4,4*4=16刚好,你得到了正确的结果sqrt(16)=4。然后你三下五除二就把程序写出来了:

float SqrtByBisection(float n) //用二分法
{ if(n<0) //小于0的按照你需要的处理 return n; float mid,last; float low,up; low=0,up=n; mid=(low+up)/2; do{if(mid*mid>n)up=mid; else low=mid;last=mid;mid=(up+low)/2; }while(abs(mid-last) > eps);//精度控制return mid;
} 

然后看看和系统函数性能和精度的差别(其中时间单位不是秒也不是毫秒,而是CPU Tick,不管单位是什么,统一了就有可比性) 

从图中可以看出,二分法和系统的方法结果上完全相同,但是性能上整整差了几百倍。为什么会有这么大的区别呢?难道系统有什么更好的办法?难道。。。。哦,对了,回忆下我们曾经的高数课,曾经老师教过我们“牛顿迭代法快速寻找平方根”,或者这种方法可以帮助我们,具体步骤如下:

求出根号a的近似值:首先随便猜一个近似值x,然后不断令x等于x和a/x的平均数,迭代个六七次后x的值就已经相当精确了。 
例如,我想求根号2等于多少。假如我猜测的结果为4,虽然错的离谱,但你可以看到使用牛顿迭代法后这个值很快就趋近于根号2了: 
(       4  + 2/4        ) / 2 = 2.25 
(     2.25 + 2/2.25     ) / 2 = 1.56944.. 
( 1.56944..+ 2/1.56944..) / 2 = 1.42189.. 
( 1.42189..+ 2/1.42189..) / 2 = 1.41423.. 
....
这种算法的原理很简单,我们仅仅是不断用(x,f(x))的切线来逼近方程x^2-a=0的根。根号a实际上就是x^2-a=0的一个正实根,这个函数的导数是2x。也就是说,函数上任一点(x,f(x))处的切线斜率是2x。那么,x-f(x)/(2x)就是一个比x更接近的近似值。代入 f(x)=x^2-a得到x-(x^2-a)/(2x),也就是(x+a/x)/2。

相关的代码如下:

float SqrtByNewton(float x)
{float val = x;//最终float last;//保存上一个计算的值do{last = val;val =(val + x/val) / 2;}while(abs(val-last) > eps);return val;
}

然后我们再来看下性能测试:

哇塞,性能提高了很多,可是和系统函数相比,还是有这么大差距,这是为什么呀?想啊想啊,想了很久仍然百思不得其解。突然有一天,我在网上看到一个神奇的方法,于是就有了今天的这篇文章,废话不多说,看代码先:

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 y0x = *(float*)&i; // convert bits BACK to floatx = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracyx = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracyx = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracyreturn 1/x;
}

然后我们最后一次来看下性能测试:

这次真的是质变了,结果竟然比系统的还要好。。。哥真的是震惊了!!!哥吐血了!!!一个函数引发了血案!!!血案,血案。。。

到现在你是不是还不明白那个“鬼函数”,到底为什么速度那么快吗?不急,先看看下面的故事吧:

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 hackingi   = 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#endifreturn 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的根。

没错,一般的求平方根都是这么循环迭代算的但是卡马克(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

参考:<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 y0x = *(float*)&i; // convert bits BACK to floatx = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracyreturn x;
}  

大家可以尝试在PC机、51、AVR、430、ARM、上面编译并实验,惊讶一下它的工作效率。

前两天有一则新闻,大意是说 Ryszard Sommefeldt 很久以前看到这么样的一段 code (可能出自 Quake III 的 source code):

float InvSqrt (float x)
{float xhalf = 0.5f*x;int i = *(int*)&x;i = 0x5f3759df - (i>>1);x = *(float*)&i;x = x*(1.5f - xhalf*x*x);return x;
}

他一看之下惊为天人,想要拜见这位前辈高人,但是一路追寻下去却一直找不到人;同时间也有其他人在找,虽然也没找到出处,但是 Chris Lomont 写了一篇论文 (in PDF) 解析这段 code 的算法 (用的是 Newton’s Method,牛顿法;比较重要的是后半段讲到怎么找出神奇的 0x5f3759df 的)。 
PS. 这个 function 之所以重要,是因为求 开根号倒数 这个动作在 3D 运算 (向量运算的部份) 里面常常会用到,如果你用最原始的 sqrt() 然后再倒数的话,速度比上面的这个版本大概慢了四倍吧… XD 
PS2. 在他们追寻的过程中,有人提到一份叫做 MIT HACKMEM 的文件,这是 1970 年代的 MIT 强者们做的一些笔记 (hack memo),大部份是 algorithm,有些 code 是 PDP-10 asm 写的,另外有少数是 C code (有人整理了一份列表)

好了,故事就到这里结束了,希望大家能有有收获:),我把源码也提供下载了,有兴趣的朋友们可以自己运行下试试看。

源码下载地址:http://diducoder.com/sotry-about-sqrt.html

转载于:https://www.cnblogs.com/kira2will/p/4014916.html

转:一个Sqrt函数引发的血案相关推荐

  1. 一个Sqrt函数引发的血案

    我们平时经常会有一些数据运算的操作,需要调用sqrt,exp,abs等函数,那么时候你有没有想过:这个些函数系统是如何实现的?就拿最常用的sqrt函数来说吧,系统怎么来实现这个经常调用的函数呢? 虽然 ...

  2. 一个“Spring轮子”引发的“血案”(4)

    事件的升级,国产技术社区中所存现出来的浮躁.世态炎凉,在我的一篇<80前>一文,终于引爆了出来. <80前>一文是一个长篇,是我吃这么20多年饭,读了不少书的一些所思.所误.所 ...

  3. 一个@Component注解引发的血案

    一个注解@Component引发的血案 首先,我们这个是用springboot架构来实现的业务 这是项目包结构和配置文件结构这是定时需要执行的任务 这是我执行PromoCodeCMCJob这个定时器的 ...

  4. 一个由正则表达式引发的血案

    阿里妹导读:周末快到了,今天为大家送上一篇很有意思的小文章,具有提神醒脑之功效.作者是来自阿里巴巴LAZADA产品技术部的申徒童鞋. 血案由来 近期我在为Lazada卖家中心做一个自助注册的项目,其中 ...

  5. 一个数字键盘引发的血案——移动端H5输入框、光标、数字键盘全假套件实现...

    https://juejin.im/post/5a44c5eef265da432d2868f6 为啥要写假键盘? 还是输入框.光标全假的假键盘? 手机自带的不用非得写个假的,吃饱没事干吧? 装逼?炫技 ...

  6. 答网友问:一个abs函数引发的问题

    某日,网友吃泡面不加开水加我好友,问了一个关于abs函数的问题.在keil中,使用abs计算浮点数的绝对值是没问题的,同样的代码,放到gcc交叉编译器中,却得不到预期结果.趁夜深人静,看了些资料,帮网 ...

  7. 记一次prometheus监控pod内存使用率错误使用sum函数引发的血案

    prometheus监控pod内存使用率 发生背景 问题伊始 根因分析 解决方案: 发生背景 pod内存使用率过高需要自动重启pod防止被kill影响线上业务 制定计算规则 首先制定的规则:(cont ...

  8. 由一个stack OOM引发的血案

    转自: https://blog.csdn.net/oscaryue/article/details/72967448 近期在App监测平台上发现如下错误信息: java.lang.OutOfMemo ...

  9. 一个缓存穿透引发的血案

    2010年9月23日,Facebook遭遇了迄今为止最严重的宕机事件之一,网站关闭了4个小时.为恢复工作,不得不让FB下线,影响了10亿用户. 在事后的故障报告中提到: 今天,我们修改了一个错误的配置 ...

最新文章

  1. python一节课多久_第一节课 python简介
  2. 复制文本朗读_原创:昭明文选配乐朗读 卷第五十一 论一 东方曼倩 非有先生论 王子渊 四子讲德论 并序...
  3. 深圳python工程师 vue_Laravel 招聘:[深圳] [15K-25K] 明源云招聘 PYTHON [SAAS] [研发基地] | Laravel China 社区...
  4. npm、cnpm、yarn的安装与常用命令
  5. 函数指针 回调函数 面向对象风格的C语言
  6. 软件开发过程中的Visio使用
  7. 数据采集的大致流程(离线和实时)
  8. 笔记本键盘失灵的修复方法
  9. 洛谷P4767 [IOI2000]邮局(决策单调DP,四边形不等式优化)
  10. 在Textview中获取指定文字位置(兼顾网址链接和emoji表情),并在其附近展示图片
  11. 华为天才少年稚晖君自制「电子」机器人!应用OpenPose,项目已开源!
  12. ContexIoT: Towards Providing Contextual Integrity to Appified IoT Platforms
  13. 若依路由刷新缓存页面 + router.push
  14. Django项目中浏览器显示127.0.0.1拒绝我们的连接请求
  15. pages文稿进去显示服务器,pages要提供服务器地址
  16. 阮一峰 React Router 教程
  17. C语言每日一练——第83天:爱因斯坦的数学题
  18. java编写数独_求用java写一个数独游戏
  19. Quorum Raft
  20. 柯杰下赢机器人_柯洁再战人工智能 大胜“一秒识人机器人”

热门文章

  1. Java自定义异常、全局捕获异常、拦截器 实现动态控制登录超时
  2. 【java】 java 高并发解决方案和高负载优化方法
  3. 【开发环境】 irun(ncverilog)无法dump fsdb波形问题解决方法
  4. 解决 screen 连接不上,提示“There is no screen to be resumed matching 18352.” 的问题
  5. VS2010无法调试问题解决
  6. 在github上托管Maven存储库
  7. 识别和非识别关系之间有什么区别?
  8. win11HDMI端口无法使用怎么办 windows11HDMI端口无法使用的解决方法
  9. Netty框架入门案例,代码示例
  10. mysql安装出现中文乱码_MySQL安装以及中文乱码问题