高效率嵌入式系统开平方根.doc

开平方根开平方根目录1. 开平方根22. 开平方根说明8181. 开平方根我们平时经常会有一些数据运算的操作，需要调用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) /小于0的按照你需要的处理 if(n 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)的切线来逼近方程x2-a=0的根。根号a实际上就是x2-a=0的一个正实根，这个函数的导数是2x。也就是说，函数上任一点(x,f(x)处的切线斜率是2x。那么，x-f(x)/(2x)就是一个比x更接近的近似值。代入 f(x)=x2-a得到x-(x2-a)/(2x)，也就是(x+a/x)/2。相关的代码如下：float SqrtByNewton(float x)/ 最终float val = x; / 保存上一个计算的值float last;dolast = 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- (i1); / 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原代码的下载地址： /file.x?fid=7547。我们知道，越底层的函数，调用越频繁。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。 论文下载地址：/clomont/Math/Papers/2003/InvSqrt.pdf ，/data/InvSqrt.pdf。最后，给出最精简的1/sqrt()函数：float InvSqrt(float x)float xhalf = 0.5f*x;int i = *(int*)&x; / get bits for floating VALUE i = 0x5f375a86- (i1); / 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 - (i1);x = *(float*)&i;x = x*(1.5f - xhalf*x*x);return x;他一看之下惊为天人，想要拜见这位前辈高人，但是一路追寻下去却一直找不到人；同时间也有其他人在找，虽然也没找到出处，但是 Chris Lomont 写了一篇论文 (in PDF) 解析这段 code 的算法 (用的是 Newtons Method，牛顿法；比较重要的是后半段讲到怎么找出神奇的 0x5f3759df 的)。 PS. 这个 function 之所以重要，是因为求 开根号倒数 这个动作在 3D 运算 (向量运算的部份) 里面常常会用到，如果你用最原始的 sqrt() 然后再倒数的话，速度比上面的这个版本大概慢了四倍吧 XD PS2. 在他们追寻的过程中，有人提到一份叫做 MIT HACKMEM 的文件，这是 1970 年代的 MIT 强者们做的一些笔记 (hack memo)，大部份是 algorithm，有些 code 是 PDP-10 asm 写的，另外有少数是 C code (有人整理了一份列表)。2. 开平方根说明人们很早就在Quake3源代码中发现了类似如下的C代码，它可以快速的求1/sqrt(x)，在3D图形向量计算方面应用很广float invSqrt(float x)float xhalf = 0.5 * x;int i = *(int*)&x; / get bits for floating valuei = 0x5f3759df - (i 1); / gives initial guessx = *(float*)&i; / convert bits back to floatx = x * (1.5 - xhalf * x * x); / Newton stepreturn x;在分析这段代码之前，先看看传统方法是怎么求一个数的平方根的倒数的，一般采用牛顿迭代法，为描述方便，假设输入数为a，显然需要满足a0，sqrt是C语言的求平方根函数，为方便起见，下面用sqrt(x)的形式代替x(1/2)求1/sqrt(a)，用迭代法，即求方程f(x)=x(-2)-a在f(x)=0时的解，选择适当的初始值x0，代入迭代式：x=x-f(x)/ f(x)化简此式得：x=3x/2-ax3/2这实际上就是上面函数倒数第二行，从函数注释也可以直接看出，这一步就是牛顿迭代，一般选择一个合适的初始值开始迭代后，迭代次数越多越接近解，换句话说就是精度越高，误差越小，当误差小于可接受值，即可获得近似结果了，就这个问题而言，初始值的选择一般要在区间(0, sqrt(3/a)，证明从略而invSqrt这个函数厉害的地方就在于，在正式迭代开始前的三行计算已经得到了一个非常接近于解的数，因此只需一次迭代，即可得到近似值，经测试，对于常用的浮点数范围，invSqrt(x)与标准解1/sqrt(x)的最大相对误差为1.75，平均相对误差为0.95，这个精度在很多时候已经满足QuakeIII的基本要求1了，而invSqrt(x)的速度则比直接计算1/sqrt(x)快4倍，这对于Quake和CS这类游戏的性能是非常重要的，而且如果需要更高的精度，则将迭代那一行再重复一次就可以了，相对误差会降到百万分之一的级别，只不过速度会慢一些，接下来我们来分析下这三行代码的原理最令人迷惑的是gives initial guess这行，即对i做了移位和减法运算，不过熟悉C语言的人应该能看出来，这个算法和float浮点数的内部表示有关，分析应该从这里入手正式开始前先轻松下，讲些历史故事，人们在QuakeIII源码发现了这个函数，于是很自然的认为这是卡马克（John Carmack）的杰作，其中0x5f3759df这个数被称为卡马克密码，我们在下面称这个数为magic，Beyond3D.com的Ryszard Sommefeldt一直在想到底是哪个家伙写了这些神奇的代码，于是就开始找作者，John Carmack在邮件回复中明确表示不是他，也不是Michael。Terje Mathisen说他写过类似的高效代码，但上面的不是。后来猜测这个来自于一些早期黑客的算法笔记，作者究竟是谁自然也难以追查了，可以肯定的是这个家伙对计算机和高数知识都有较好理解，很聪明2003年普渡大学的数学家Chris Lomont写了一篇文章对这段代码进行了分析。论文是英文的，地址在：/view/80b84d1fb7360b4c2e3f644b.html在这篇12页的论文中，Lomont对这个算法做了分析，并从推导出了一个理论上最优的magic 0x5f37642f，有意思的是，这个数居然没有invSqrt里的0x5f3759df效果好，最大相对误差达到1.78，Lomont一怒之下，用暴力搜索枚举了所有可能的magic，终于找到一个最优的magic 0x5f375a86，只比0x5f3759df效果好一点点，至于invSqrt的作者究竟如何找到0x5f3759df的，也就是个迷了开始正式分析，这三行代码是把float在内存中的表示作为一个整数i看待，然后对i进行一次移位和减法，然后再将i的值作为一个float看待，所以我们先看看float在内存中的表示，一般计算机的浮点数遵循IEEE754标准，采用以2为底数的科学计数法，例如二进制的11010.11001记为1.101011001*10100，float和int都是32位，占4个字节（注意，这个函数早期的代码中整数类型应该是long，因为那时候在dos下，int是16位的）：最高位d31：符号位，0表示非负，1表示非正，为什么不直接说正负呢，因为有数值0的存在，浮点数有+0和-0的区别，这个位用S表示d30d23：指数域，存放一个整数，表示127+E，E为指数，由于指数域的范围是0255，因此理论上可以表示的指数范围是-127128，不过0和255有特殊含义，所以范围实际要稍微小一点，这个先按下不表，我们认为常用浮点数不包括这两种极端情况d22d0：有效数字域，只是小数部分，由于科学计数法的规定，整数部分肯定是1，就省略了，这样可以避免不必要的精度浪费，为描述方便，这个域所表示的小数设为F。当然有人会问，那0怎么办，+0和-0有自身的特殊表示法，S位表示符号，其他位都为0的时候是+0或-0于是除去0和IEEE754规定的特殊值，一个常用浮点数的表示可以看做：(-1)S*(1+F)*2E，具体到我们需要分析的问题，由于输入是正数，S位肯定是0，就不做考虑了，简化为：(1+F)*2E好，现在我们需要求(1+F)*2E的平方根的倒数，即求1/sqrt(1+F)*2E)，求得的结果当然也要用这个浮点数表示法，有效数字必须在1, 2)，指数域为整数，则结果分两种情况推导出结果：E为奇数：sqrt(2/(1+F)*2(-(E+1)/2)E为偶数：2/sqrt(1+F)*2(-E/2-1)先看指数，如果我们需要通过计算机的整数运算（移位、位运算和加减法等）来逼近解，首先要在数量级上尽量靠近，因为只要数量级一样，两个数的误差范围是最小的，也就是说，我们需要：E为奇数时，将127+E变成127-(E+1)/2E为偶数时，将127+E变成127-E/2-1于是，通过右移一位来实现除以二，通过用一个数减去指数域来将E变成负的，这样invSqrt中的那句就很容易理解了：E为奇数时，127+E为偶数，右移等于除以二，190-(127+E)/2 = 190-63-(E+1)/2 = 127-(E+1)/2E为偶数时，127+E为奇数，右移等于先减一再除以二，189-(127+E-1)/2 = 189-63-E/2 = 127-E/2-1然后我们把invSqrt中的magic 0x5f3759df用float的形式展开，则其指数域为0x5f1这个操作将指数域最后一位的1（即上面说的“先减一”）也向右移了一位，i右移后d22位为1，这样一来只要被减数的d22这一位是0，就会因为不够减而产生借位，指数域被借了1，自然就变成189了细心的童鞋应该发现了，当E为偶数时，减法做完后，指数域一定是127-E/2-1，但是如果E为奇数，则不保证指数域是127-(E+1)/2，因为这时候被减数和减数的d22位都是0，但如果减数的d21到d0这个域的数字比被减数的大，就会产生借位而使得指数域比预期的127-(E+1)/2要小1，这个问题会导致一定的误差，但是在最后的迭代中误差会被缩小，这个误差具体有多大，这里就不详细讨论了，有兴趣的童鞋可以自己算算看我们还是先证明用这种方法得到的x0是落在上述区间(0, sqrt(3/a)的，这里a就是(1+F)*2E由于输入是一个非负数，则S为0，指数域肯定不为0（右移后不可能刚好为190），因此x0肯定大于0，我们将sqrt(3/a)展开成期望的解的指数的乘法形式：E为奇数时，展开为sqrt(6/(1+F)*2(-(E+1)/2)，如果按上面说的那种情况产生借位，则展开为sqrt(24/(1+F)*2(-(E+3)/2)E为偶数时，展开为sqrt(12/(1+F)*2(-E/2-1)可以看到，无论是哪种情况，在指数相同的情况下，有效数字都大于2，反过来说，用上述算法得到的x0是小于sqrt(3/a)的，而且还小了很多，非常接近解，这时候只需要一次迭代，就得到了误差很小的近似结果到这里，基本原理都清楚了，只要我们保证x0在这个区间中，再做迭代总是能进一步接近解的，现在的问题就在于magic中d21到d0这个域的值应该怎么取了，这个取值关系到每次迭代的误差，比如说，我们全取1，这样也避免了E为奇数时的借位情况，这样magic就是0x5f3fffff，用这个magic测试，结果最大误差超过1%，平均误差超过6，显然效果太差了我们假设将magic作为float看时，小数部分的值是M，由于d22位已确定是0，则0=M1)的减法运算中，d22到d0域的运算可看做是定点小数减法，分三种情况讨论：E为奇数时，减数的d22位为0，则小数部分的值为F/2，因为右移对小数来说也是除以二，d0位如果是1则会舍弃，这个因为太小而忽略，假设M=F/2，则不需要向指数域借位，计算结果的小数部分为M-F/2，指数域符合结果预期，此时的相对误差为：rd1(M,F)=|1-(1+M-F/2)/sqrt(2/(1+F)|，0=F=2M假设MF/2，则需要向指数域借位，计算结果小数部分为1+M-F/2，指数位比预期低1，计算误差的时候需要将指数域差值补回去，此时的相对误差：rd2(M,F)=|1-(2+M-F/2)/2/sqrt(2/(1+F)|，2MF1E为偶数时，减数的d22位为1，则小数部分的值为F/2+1/2，此时必定借位，相对误差为：rd3(M,F)=|1-(2+M-F/2-1/2)/2/sqrt(1+F)|，0=F 1;x = *(float *)&i;x = (x + a / x) * 0.5;return x;推导过程从略，有兴趣的童鞋可以自行研究，注意这里是先加再右移，所以用unsigned int防止负数右移，这个mySqrt计算平方根速度比invSqrt(a)*a快，而且最大误差只有0.6，那么，有没有更快的算法呢，如果单纯用计算，可能很难超越mySqrt了，想要更快，得从另外的方向想办法比较直接的想法是，造一个巨大的table，储存a到sqrt(a)的映射关系，这样不需要计算，只要查表就行，而且由于