一个Sqrt谋杀触发功能

我们平时常常会有一些数据运算的操作，须要调用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.. 
....
这样的算法的原理非常easy，我们不过不断用(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 y0
	x = *(float*)&i; // convert bits BACK to float
	x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
	x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
	x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy

	return 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 )

{

        longi;

        floatx2, y;

        constfloat threehalfs = 1.5F;

 

        x2= number * 0.5F;

        y   = number;

        i   = * ( long * ) &y;   // evil floating point bit level hacking

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

        y   = * ( float * ) &i;

        y   = y * ( threehalfs - ( x2 * y * y ) ); //1st iteration

        //y   = y * ( threehalfs - ( x2 * y * y )); // 2nd iteration, this can be removed

 

        returny;

} 

函数返回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

最后，给出最精简的1/sqrt()函数：

float InvSqrt(float x)

{

        floatxhalf = 0.5f*x;

        inti = *(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

        returnx;

}

 

前两天有一则新闻，大意是说 Ryszard Sommefeldt 非常久曾经看到这么样的一段code (可能出自 Quake III 的 source code)：

float InvSqrt (float x)

{

        floatxhalf = 0.5f*x;

        inti = *(int*)&x;

        i= 0x5f3759df - (i>>1);

        x= *(float*)&i;

        x= x*(1.5f - xhalf*x*x);

        returnx;

}
PS. 这个 function 之所以重要，是由于求开根号倒数这个动作在 3D 运算 (向量运算的部份) 里面经常会用到，假设你用最原始的 sqrt() 然后再倒数的话。速度比上面的这个版本号大概慢了四倍吧… XD 
源代码下载地址：http://diducoder.com/sotry-about-sqrt.html

 

 

 

 

问题出在我签入的来自卡马克的求平方根函数代码。

double InvSqrt(double number)
{
      __int64 i;
      double x2, y;
      const double threehalfs = 1.5F;
    
      x2 = number * 0.5F;
      y = number;
      i = *(__int64 *)&y;
      i = 0x5fe6ec85e7de30da - (i >> 1);
      y = *( double *)&i;
      y = y * (threehalfs - (x2 * y * y)); //1stiteration
      y = y * (threehalfs - (x2 * y * y)); //2nditeration, this can be removed
      return y;
}

红色部分代码在gcc开启-fstrict-aliasing选项后将得到错误的代码。因为使用了type-punned pointer将打破strict-aliasing规则。

因为-fstrict-aliasing选项在-O2, -O3, -Os等优化模式下都将开启（眼下dev不带优化。main带-O3所以该问题仅仅在main上出现）所以建议对linux编译中产生
warning:dereferencing type-punned pointer will break strict-aliasing rule
警告的情况作为编译失败。以便防止出现类似问题。
上述代码应当使用联合体重写为：

double InvSqrt(double number)
{
      double x2, y;
      const double threehalfs = 1.5F;
      union
      {
            double d;
            __int64 i;
      }d;
      x2 = number * 0.5F;
      y = number;
      d.d =  y;
      d.i = 0x5fe6ec85e7de30da - (d.i >> 1);
      y = d.d;
      y = y * (threehalfs - (x2 * y * y)); //1stiteration
      y = y * (threehalfs - (x2 * y * y)); //2nditeration, this can be removed
      return y;
}

这样就不会打破该规则。

 

 

 

平方根倒数速算法

平方根倒数速算法（英语：Fast Inverse Square Root，亦常以“FastInvSqrt()”或其使用的十六进制常数0x5f3759df代称）是用于高速计算（即的平方根的倒数，在此需取符合IEEE 754标准格式的32位浮点数）的一种算法。

此算法最早可能是于90年代前期由SGI所发明。后来则于1999年在《雷神之锤III竞技场》的源码中应用。但直到2002－2003年间才在Usenet一类的公共论坛上出现。这一算法的优势在于降低了求平方根倒数时浮点运算操作带来的巨大的运算耗费，而在计算机图形学领域，若要求取照明和投影的波动角度与反射效果，就常需计算平方根倒数。

此算法首先接收一个32位带符浮点数，然后将之作为一个32位整数看待，以将其向右进行一次逻辑移位的方式将之取半，并用十六进制“魔术数字”0x5f3759df减之，如此就可以得对输入的浮点数的平方根倒数的首次近似值。而后又一次将其作为浮点数，以牛顿法重复迭代。以求出更精确的近似值，直至求出符合准确度要求的近似值。

在计算浮点数的平方根倒数的同一精度的近似值时，此算法比直接使用浮点数除法要快四倍。

1介绍

平方根倒数速算法最早被觉得是由约翰·卡马克所发明。但后来的调查显示，该算法在这之前就于计算机图形学的硬件与软件领域有所应用，如SGI和3dfx就曾在产品中应用此算法。而就如今所知，此算法最早由Gary Tarolli在SGIIndigo的开发中使用。虽说在随后的相关研究中也提出了一些可能的来源。但至今为止仍未能确切知晓此常数的起源。

2算法的切入点

浮点数的平方根倒数经常使用于计算正规化矢量。^[1] 3D图形程序须要使用正规化矢量来实现光照和投影效果。因此每秒都需做上百万次平方根倒数运算，而在处理坐标转换与光源的专用硬件设备出现前。这些计算都由软件完毕，计算速度亦相当之慢；在1990年代这段代码开发出来之时，多数浮点数操作的速度更是远远滞后于整数操作。因而针对正规化矢量算法的优化就显得尤为重要。

以下陈述计算正规化矢量的原理：

要将一个矢量标准化，就必须计算其欧几里得范数以求得矢量长度。而这时就需对矢量的各分量的平方和求平方根；而当求取到其长度并以之除该矢量的每一个分量后。所得的新矢量就是与原矢量同向的单位矢量，若以公式表示：

·    可求得矢量v的欧几里得范数。此算法正类如对欧几里得空间的两点求取其欧几里得距离。

·    而求得的就是标准化的矢量，若以代表。则有，

可见标准化矢量时须要用到对矢量分量的平方根倒数计算。所以对平方根倒数计算算法的优化对计算正规化矢量也大有裨益。

为了加速图像处理单元计算，《雷神之锤III竞技场》使用了平方根倒数速算法，而后来採用现场可编程逻辑门阵列的顶点着色器也应用了此算法。^[2] 

3将浮点数转化为整数

要理解这段代码。首先需了解浮点数的存储格式。

一个浮点数以32个二进制位表示一个有理数，而这32位由其意义分为三段：首先首位为符号位，如若是0则为正数，反之为负数。接下来的8位表示经过偏移处理（这是为了使之能表示-127－128）后的指数。最后23位表示的则是有效数字中除最高位以外的其余数字。

将上述结构表示成公式即为，当中表示有效数字的尾数（此处，偏移量，而指数的值决定了有效数字（在Lomont和McEniry的论文中称为“尾数”（mantissa））代表的是小数还是整数。以上图为例。将描写叙述带入有），且，则可得其表示的浮点数为。

符号位
0	1	1	1	1	1	1	1	=	127
0	0	0	0	0	0	1	0	=	2
0	0	0	0	0	0	0	1	=	1
0	0	0	0	0	0	0	0	=	0
1	1	1	1	1	1	1	1	=	−1
1	1	1	1	1	1	1	0	=	−2
1	0	0	0	0	0	0	1	=	−127
1	0	0	0	0	0	0	0	=	−128

8位二进制整数补码演示样例

如上所述。一个有符号正整数在二进制补码系统中的表示中首位为0，而后面的各位则用于表示其数值。将浮点数取别名存储为整数时。该整数的数值即为，当中E表示指数，M表示有效数字；若以上图为例，图中例子若作为浮点数看待有。，则易知其转化而得的整数型号数值为。

因为平方根倒数函数仅能处理正数，因此浮点数的符号位（即如上的Si）必为0，而这就保证了转换所得的有符号整数也必为正数。

以上转换就为后面的计算带来了可行性，之后的第一步操作（逻辑右移一位）即是使该数的长整形式被2所除。^[3] 

4历史与考究

id Software创始人约翰·卡马克

《雷神之锤III》的代码直到QuakeCon 2005才正式放出，但早在2002年（或2003年）时平方根倒数速算法的代码就已经出如今Usenet与其它论坛上了。最初人们推測是卡马克写下了这段代码，但他在询问邮件的回复中否定了这个观点，并推測可能是先前曾帮id Software优化雷神之锤的资深汇编程序猿Terje Mathisen写下了这段代码；而在Mathisen的邮件里他表示在1990年代初他仅仅曾作过类似的实现。确切来说这段代码亦非他所作。如今所知的最早实现是由Gary Tarilli在SGIIndigo中实现的，但他亦坦承他仅对常数R的取值做了一定的改进，实际上他也不是作者。

RysSommefeldt则在向以发明MATLAB而闻名的CleveMoler查证后觉得原始的算法是Ardent Computer公司的GregWalsh所发明，但他也没有不论什么决定性的证据能证明这一点。

眼下不仅该算法的原作者不明，人们也仍无法明白当初选择这个“魔术数字”的方法。Chris Lomont在研究中曾做了个试验：他编写了一个函数，以在一个范围内遍历选取R值的方式将逼近误差降到最小，以此方法他计算出了线性近似的最优R值0x5f37642f（与代码中使用的0x5f3759df相当接近）。但以之代入算法计算并进行一次牛顿迭代后，所得近似值与代入0x5f3759df的结果相比精度却仍稍微更低；而后Lomont将目标改为遍历选取在进行1－2次牛顿迭代后能得到最大精度的R值。并由此算出最优R值为0x5f375a86，以此值代入算法并进行牛顿迭代后所得的结果都比代入原始值（0x5f3759df）更精确，于是他的论文最后以“原始常数是以数学推导还是以重复试错的方式求得”的问题作结。

在论文中Lomont亦指出64位的IEEE754浮点数（即双精度类型）所相应的魔术数字是0x5fe6ec85e7de30da，但后来的研究表明代入0x5fe6eb50c7aa19f9的结果准确度更高（McEniry得出的结果则是0x5FE6EB50C7B537AA，精度介于两者之间）。

在Charles McEniry的论文中，他使用了一种类似Lomont但更复杂的方法来优化R值：他最開始使用穷举搜索法。所得结果与Lomont同样。而后他尝试用带权二分法寻找最优值。所得结果恰是代码中所使用的魔术数字0x5f3759df，因此McEniry确信这一常数也许最初便是以“在可容忍误差范围内使用二分法”的方式求得。^[4] 

5凝视

1. 因为现代计算机系统对长整型的定义有所差异，使用长整型会减少此段代码的可移植性。

详细来说。由此段浮点转换为长整型的定义可知，如若这段代码正常执行，所在系统的长整型长度应为4字节（32位），否则又一次转为浮点数时可能会变成负数；而因为C99标准的广泛应用，在现今多数64位计算机系统（除使用LLP64数据模型的Windows外）中。长整型的长度都是8字节。^[5] 

2. 此处“浮点数”所指为标准化浮点数，也即有效数字部分必须满足，可參见DavidGoldberg. What Every Computer Scientist Should Know About Floating-PointArithmetic. ACM Computing Surveys. 1991.March, 23 (1): 5–48. doi:10.1145/103162.103163.

3. Lomont 2003确定R的方式则有所不同，他先将R分解为与，当中与分别代表R的有效数字域和指数域。^[6] 

 

 

 

 

 

卡马克高速平方根---平方根倒数算法 [转]

--------------------------------------------------------------------------------
 高速平方根（平方根倒数）算法

日前在书上看到一段使用多项式逼近计算平方根的代码，至今都没搞明确作者是如何推算出那个公式的。

但在尝试解决这个问题的过程中，学到了不少东西，于是便有了这篇心得，写出来和大家共享。

当中有错漏的地方，还请大家多多不吝赐教。

的确，正如很多人所说的那样，如今有有FPU，有3DNow。有SIMD，讨论软件算法好像不合时宜。

关于sqrt的话题事实上早在2003年便已在 GameDev.net上得到了广泛的讨论（可见我实在很火星了，当然不排除还有其它尚在冥王星的人，嘿嘿）。而尝试探究该话题则全然是出于本人的兴趣和好奇心（换句话说就是无知）。

我仅仅是个beginner，所以这样的大是大非的问题我也说不清楚（在GameDev.net上也有非常多类似的争论）。

但不管怎样，Carmack在DOOM3中还是使用了软件算法。而多知道一点数学知识对3D编程来说也仅仅有优点没坏处。3D图形编程事实上就是数学，数学。还是数学。

=========================================================

在3D图形编程中，常常要求平方根或平方根的倒数，比如：求向量的长度或将向量归一化。C数学函数库中的sqrt具有理想的精度，但对于3D游戏程式来说速度太慢。我们希望可以在保证足够的精度的同一时候。进一步提快速度。

Carmack在QUAKE3中使用了以下的算法。它第一次在公众场合出现的时候，差点儿震住了全部的人。据说该算法事实上并非Carmack发明的，它真正的作者是Nvidia的GaryTarolli（未经证实）。

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
//
// 计算參数x的平方根的倒数
//
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;
}
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
该算法的本质事实上就是牛顿迭代法（Newton-RaphsonMethod。简称NR）。而NR的基础则是泰勒级数（Taylor Series）。

NR是一种求方程的近似根的方法。首先要预计一个与方程的根比較靠近的数值。然后依据公式推算下一个更加近似的数值，不断反复直到能够获得惬意的精度。

其公式例如以下：
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
函数：y=f(x)
其一阶导数为：y'=f'(x)
则方程：f(x)=0 的第n+1个近似根为
x[n+1] = x[n] - f(x[n]) / f'(x[n])
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
 NR最关键的地方在于预计第一个近似根。

假设该近似根与真根足够靠近的话，那么仅仅须要少数几次迭代，就能够得到惬意的解。

如今回过头来看看怎样利用牛顿法来解决我们的问题。求平方根的倒数，实际就是求方程1/(x^2)-a=0的解。将该方程按牛顿迭代法的公式展开为：
x[n+1]=1/2*x[n]*(3-a*x[n]*x[n])
 将1/2放到括号中面，就得到了上面那个函数的倒数第二行。

接着，我们要设法预计第一个近似根。这也是上面的函数最奇妙的地方。

它通过某种方法算出了一个与真根很接近的近似根，因此它仅仅须要使用一次迭代过程就获得了较惬意的解。它是如何做到的呢？全部的奥妙就在于这一行：
i = 0x5f3759df - (i >> 1); // 计算第一个近似根

超级莫名其妙的语句。不是吗？但细致想一下的话。还是能够理解的。我们知道，IEEE标准下，float类型的数据在32位系统上是这样表示的（大体来说就是这样，但省略了非常多细节，有兴趣能够GOOGLE）：
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
bits：31 30... 0
31：符号位
30-23：共8位，保存指数（E）
22-0：共23位，保存尾数（M）
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
所以。32位的浮点数用十进制实数表示就是：M*2^E。开根然后倒数就是：M^(-1/2)*2^(-E/2)。如今就十分清晰了。语句i> >1其工作就是将指数除以2。实现2^(E/2)的部分。

而前面用一个常数减去它，目的就是得到M^(1/2)同一时候反转全部指数的符号。

至于那个0x5f3759df，呃，我仅仅能说，的确是一个超级的MagicNumber。

那个Magic Number是能够推导出来的，但我并不打算在这里讨论，由于实在太繁琐了。简单来说。其原理例如以下：由于IEEE的浮点数中，尾数M省略了最前面的1。所以实际的尾数是1+M。

假设你在大学上数学课没有打瞌睡的话，那么当你看到(1+M)^(-1/2)这种形式时。应该会立即联想的到它的泰勒级数展开，而该展开式的第一项就是常数。

以下给出简单的推导过程：
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
对于实数R>0，如果其在IEEE的浮点表示中，
指数为E，尾数为M，则：

R^(-1/2)
= (1+M)^(-1/2) * 2^(-E/2)

将(1+M)^(-1/2)按泰勒级数展开，取第一项，得：

原式
= (1-M/2) * 2^(-E/2)
= 2^(-E/2) - (M/2) * 2^(-E/2)

假设不考虑指数的符号的话。
(M/2)*2^(E/2)正是(R>>1)。
而在IEEE表示中。指数的符号仅仅需简单地加上一个偏移就可以，
而式子的前半部分刚好是个常数，所以原式能够转化为：

原式 = C - (M/2)*2^(E/2) = C -(R>>1)，当中C为常数

所以仅仅须要解方程：
R^(-1/2)
= (1+M)^(-1/2) * 2^(-E/2)
= C - (R>>1)
求出令到相对误差最小的C值就能够了
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
上面的推导过程仅仅是我个人的理解，并未得到证实。而ChrisLomont则在他的论文中具体讨论了最后那个方程的解法。并尝试在实际的机器上寻找最佳的常数C。有兴趣的朋友能够在文末找到他的论文的链接。

所以。所谓的Magic Number，并非从N元宇宙的某个星系因为时空扭曲而掉到地球上的，而是几百年前就有的数学理论。仅仅要熟悉NR和泰勒级数，你我相同有能力作出类似的优化。

在GameDev.net上有人做过測试，该函数的相对误差约为0.177585%，速度比C标准库的sqrt提高超过20%。假设添加一次迭代过程，相对误差能够减少到e-004 的级数。但速度也会降到和sqrt差点儿相同。

据说在DOOM3中。Carmack通过查找表进一步优化了该算法，精度近乎完美。并且速度也比原版提高了一截（正在努力弄源代码。谁有发我一份）。

值得注意的是。在Chris Lomont的演算中。理论上最棒的常数（精度最高）是0x5f37642f，而且在实际測试中，假设仅仅使用一次迭代的话，其效果也是最好的。但奇怪的是。经过两次NR后，在该常数下解的精度将减少得很厉害（天知道是怎么回事！

）。经过实际的測试，Chris Lomont觉得。最棒的常数是0x5f375a86。

假设换成64位的double版本号的话。算法还是一样的，而最优常数则为0x5fe6ec85e7de30da（又一个令人冒汗的Magic Number - -b）。

这个算法依赖于浮点数的内部表示和字节顺序。所以是不具移植性的。假设放到Mac上跑就会挂掉。

假设想具备可移植性。还是乖乖用sqrt好了。但算法思想是通用的。大家能够尝试推算一下对应的平方根算法。

以下给出Carmack在QUAKE3中使用的平方根算法。

Carmack已经将QUAKE3的全部源码捐给开源了。所以大家能够放心使用，不用操心会受到律师信。
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
//
// Carmack在QUAKE3中使用的计算平方根的函数
//
float CarmSqrt(float x){
union{
int intPart;
float floatPart;
} convertor;
union{
int intPart;
float floatPart;
} convertor2;
convertor.floatPart = x;
convertor2.floatPart = x;
convertor.intPart = 0x1FBCF800 + (convertor.intPart >> 1);
convertor2.intPart = 0x5f3759df - (convertor2.intPart >> 1);
return 0.5f*(convertor.floatPart + (x * convertor2.floatPart));
}

 

 

 

日前在书上看到一段使用多项式逼近计算平方根的代码，至今都没搞明确作者是如何推算出那个公式的。

但在尝试解决这个问题的过程中，学到了不少东西，于是便有了这篇心得，写出来和大家共享。

当中有错漏的地方，还请大家多多不吝赐教。 

的确，正如很多人所说的那样，如今有有FPU，有3DNow。有SIMD，讨论软件算法好像不合时宜。关于sqrt的话题事实上早在2003年便已在 GameDev.net上得到了广泛的讨论（可见我实在很火星了，当然不排除还有其它尚在冥王星的人。嘿嘿）。而尝试探究该话题则全然是出于本人的兴趣和好奇心（换句话说就是无知）。 

我仅仅是个beginner，所以这样的大是大非的问题我也说不清楚（在GameDev.net上也有非常多类似的争论）。

但不管怎样。Carmack在DOOM3中还是使用了软件算法，而多知道一点数学知识对3D编程来说也仅仅有优点没坏处。3D图形编程事实上就是数学，数学，还是数学。 

========================================================= 

在3D图形编程中，常常要求平方根或平方根的倒数。比如：求向量的长度或将向量归一化。

C数学函数库中的sqrt具有理想的精度。但对于3D游戏程式来说速度太慢。

我们希望可以在保证足够的精度的同一时候，进一步提快速度。 

Carmack在QUAKE3中使用了以下的算法。它第一次在公众场合出现的时候，差点儿震住了全部的人。据说该算法事实上并非Carmack发明的，它真正的作者是Nvidia的Gary Tarolli（未经证实）。 
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 
// 
// 计算參数x的平方根的倒数 
// 
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; 
} 
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 
该算法的本质事实上就是牛顿迭代法（Newton-Raphson Method，简称NR）。而NR的基础则是泰勒级数（TaylorSeries）。NR是一种求方程的近似根的方法。

首先要预计一个与方程的根比較靠近的数值，然后依据公式推算下一个更加近似的数值。不断反复直到能够获得惬意的精度。

其公式例如以下： 
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 
函数：y=f(x) 
其一阶导数为：y'=f'(x) 
则方程：f(x)=0 的第n+1个近似根为 
x[n+1] = x[n] - f(x[n]) / f'(x[n]) 
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 
NR最关键的地方在于预计第一个近似根。假设该近似根与真根足够靠近的话，那么仅仅须要少数几次迭代，就能够得到惬意的解。 

如今回过头来看看怎样利用牛顿法来解决我们的问题。求平方根的倒数，实际就是求方程1/(x^2)-a=0的解。将该方程按牛顿迭代法的公式展开为： 
x[n+1]=1/2*x[n]*(3-a*x[n]*x[n]) 
将1/2放到括号中面，就得到了上面那个函数的倒数第二行。

 

接着。我们要设法预计第一个近似根。

这也是上面的函数最奇妙的地方。它通过某种方法算出了一个与真根很接近的近似根。因此它仅仅须要使用一次迭代过程就获得了较惬意的解。

它是如何做到的呢？全部的奥妙就在于这一行： 
i = 0x5f3759df - (i >> 1); // 计算第一个近似根 

超级莫名其妙的语句。不是吗？但细致想一下的话，还是能够理解的。我们知道。IEEE标准下，float类型的数据在32位系统上是这样表示的（大体来说就是这样，但省略了非常多细节。有兴趣能够GOOGLE）： 
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 
bits：3130 ... 0 
31：符号位 
30-23：共8位。保存指数（E） 
22-0：共23位。保存尾数（M） 
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 
所以。32位的浮点数用十进制实数表示就是：M*2^E。开根然后倒数就是：M^(-1/2)*2^(-E/2)。如今就十分清晰了。

语句i> >1其工作就是将指数除以2。实现2^(E/2)的部分。而前面用一个常数减去它，目的就是得到M^(1/2)同一时候反转全部指数的符号。 

至于那个0x5f3759df。呃，我仅仅能说，的确是一个超级的Magic Number。 

那个MagicNumber是能够推导出来的。但我并不打算在这里讨论，由于实在太繁琐了。简单来说，其原理例如以下：由于IEEE的浮点数中。尾数M省略了最前面的1。所以实际的尾数是1+M。假设你在大学上数学课没有打瞌睡的话，那么当你看到(1+M)^(-1/2)这种形式时，应该会立即联想的到它的泰勒级数展开，而该展开式的第一项就是常数。以下给出简单的推导过程： 
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 
对于实数R>0。如果其在IEEE的浮点表示中。 
指数为E，尾数为M。则： 

R^(-1/2) 
= (1+M)^(-1/2) * 2^(-E/2) 

将(1+M)^(-1/2)按泰勒级数展开，取第一项，得： 

原式 
= (1-M/2) * 2^(-E/2) 
= 2^(-E/2) - (M/2) * 2^(-E/2) 

假设不考虑指数的符号的话， 
(M/2)*2^(E/2)正是(R>>1)， 
而在IEEE表示中。指数的符号仅仅需简单地加上一个偏移就可以。 
而式子的前半部分刚好是个常数，所以原式能够转化为： 

原式= C - (M/2)*2^(E/2) = C - (R>>1)，当中C为常数 

所以仅仅须要解方程： 
R^(-1/2) 
= (1+M)^(-1/2) * 2^(-E/2) 
= C - (R>>1) 
求出令到相对误差最小的C值就能够了 
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 
上面的推导过程仅仅是我个人的理解，并未得到证实。而Chris Lomont则在他的论文中具体讨论了最后那个方程的解法，并尝试在实际的机器上寻找最佳的常数C。有兴趣的朋友能够在文末找到他的论文的链接。

 

所以，所谓的Magic Number，并非从N元宇宙的某个星系因为时空扭曲而掉到地球上的。而是几百年前就有的数学理论。

仅仅要熟悉NR和泰勒级数，你我相同有能力作出类似的优化。

 

在GameDev.net上有人做过測试，该函数的相对误差约为0.177585%。速度比C标准库的sqrt提高超过20%。

假设添加一次迭代过程，相对误差能够减少到e-004 的级数。但速度也会降到和sqrt差点儿相同。据说在DOOM3中。Carmack通过查找表进一步优化了该算法。精度近乎完美，并且速度也比原版提高了一截（正在努力弄源代码，谁有发我一份）。 

值得注意的是，在Chris Lomont的演算中，理论上最棒的常数（精度最高）是0x5f37642f，而且在实际測试中，假设仅仅使用一次迭代的话，其效果也是最好的。

但奇怪的是，经过两次NR后，在该常数下解的精度将减少得很厉害（天知道是怎么回事！）。经过实际的測试，Chris Lomont觉得。最棒的常数是0x5f375a86。假设换成64位的double版本号的话，算法还是一样的。而最优常数则为0x5fe6ec85e7de30da（又一个令人冒汗的Magic Number - -b）。 

这个算法依赖于浮点数的内部表示和字节顺序。所以是不具移植性的。假设放到Mac上跑就会挂掉。假设想具备可移植性，还是乖乖用sqrt好了。

但算法思想是通用的。大家能够尝试推算一下对应的平方根算法。 

以下给出Carmack在QUAKE3中使用的平方根算法。

Carmack已经将QUAKE3的全部源码捐给开源了，所以大家能够放心使用。不用操心会受到律师信。

 
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 
// 
// Carmack在QUAKE3中使用的计算平方根的函数 
// 
float CarmSqrt(float x){ 
union{ 
int intPart; 
float floatPart; 
} convertor; 
union{ 
int intPart; 
float floatPart; 
} convertor2; 
convertor.floatPart = x; 
convertor2.floatPart = x; 
convertor.intPart = 0x1FBCF800 + (convertor.intPart >> 1); 
convertor2.intPart = 0x5f3759df - (convertor2.intPart >> 1); 
return 0.5f*(convertor.floatPart + (x * convertor2.floatPart)); 
} 

 

 

 

 

1.     #include <iostream> #include <math.h>

2.     using namespace std;

3.     

4.     float InvSqrt (float x)

5.     {

6.             float xhalf = 0.5f*x;

7.             int i = *(int*)&x;

8.             i = 0x5f3759df - (i >> 1); // 计算第一个近似根

9.             x = *(float*)&i;

10.           x = x*(1.5f - xhalf*x*x); // 牛顿迭代法

11.           return x;

12.   }

13.   

14.   float CarmSqrt(float x){

15.           union{

16.                   int intPart;

17.                   float floatPart;

18.           } convertor;

19.           union{

20.                   int intPart;

21.                   float floatPart;

22.           } convertor2;

23.           convertor.floatPart = x;

24.           convertor2.floatPart = x;

25.           convertor.intPart = 0x1FBCF800 + (convertor.intPart >> 1);

26.           convertor2.intPart = 0x5f3759df - (convertor2.intPart >> 1);

27.           return 0.5f*(convertor.floatPart + (x * convertor2.floatPart));

28.   }

29.   

30.   int main()

31.   {

32.           float a = 0.25;

33.   

34.           cout << InvSqrt(a) << endl;

35.           cout << CarmSqrt(a) << endl;

36.           cout << sqrt(a) << endl;

37.   

38.           return 0;

39.   }

结果例如以下：

普通浏览复制代码

1.  1.99661     

2.  0.488594

3.  0.5

4.  请按随意键继续. . .

代码没错，仅仅是cout << InvSqrt(a) << endl;应该改成cout <<1.0f / InvSqrt(a) << endl;。所以我提到了就这个样例来说，精度已经非常不错了。

 

 

 

 

主要思想利用了浮点数的表示法，ieee的标准当中关键部分是

int i = *(int*)&x; 
i = 0x5f3759df - (i >> 1); // 计算第一个近似根 
x = *(float*)&i; 
对于0x5f3759df假设看成浮点数的话，其指数位为(0101)(1111)0->1011111=190

设原来的xIEEE表示的指数位为a。则新的结果的指数位变成了190-a/2

假设要转换成 m*2^n格式，须要将指数a-127。所以得到的新结果用数学表示的指数位为

190-a/2-127= 63-a/2 

而x的原来假设从ieee格式变成数学表示应为a-127。比較63-a/2与a-127。能够看出指数之间的关系为

63-a/2==（a-127)*(-1/2)。

x应该已经是1/sqrt(x)的近似值

 

 

 

首先看牛顿迭代法的原理

 

若函数

f（x）在点的某一临域内具有直到（n+1）阶导数，则在该邻域内f（x）的n阶泰勒公式为：

f(x) = f(x0)+(x－x0)f'(x0)+(x－x0)^2*f''(x0)/2! +… +...fn(x0)(x- x0)n/n!....

解非线性方程f(x)=0的牛顿法是把非线性方程线性化的一种近似方法。把f(x)在x0点附近展开成泰勒级数 f(x) = f(x0)+(x－x0)f'(x0)+(x－x0)^2*f''(x0)/2! +…

 

取其线性部分。作为非线性方程f(x) = 0的近似方程，即泰勒展开的前两项，则有

f(x0)+f'(x0)(x－x0)=f(x)=0 

设f'(x0)≠0则f(x0)+f'(x0)(x－x0)=0 

其解为x1=x0－f(x0)/f'(x0) 

这样，得到牛顿法的一个迭代序列：x(n+1)=x(n)－f(x(n))/f'(x(n))。

也就是说x(n+1)=x(n)－f(x(n))/f'(x(n))。是用来计算f(x)=0的根的迭代公式，

假设x(n)与真根足够靠近的话，那么仅仅须要少数几次迭代，就能够得到惬意的解。

 

我们如今要求的是一个数的平方根的倒数，设这个数为a，那么我们要计算的就是x= a^(-0.5)

。首先我们要把这个计算需求转换为f(x)=0的形式。

 

x = a^(-0.5) 

x^(-2) = a 

x^(-2)-a = 0 

这样就转化成了f(x) = 0的形式。当中x就是a的平方根的倒数。

 

f(x) = x^(-2)-a 则f'(x) = -2*x^(-3) 依据牛顿迭代公式能够得到

 

x(n+1) = x(n) - (x(n)^(-2) - a)/(-2*x(n)^(-3)) 

x(n+1) = x(n) + (x(n)^(-2) - a)/(2*x(n)^(-3)) 

x(n+1) = x(n) + (x(n)/2) - (a*x(n)^3)/2 

x(n+1) = 1.5*x(n) - (0.5*a*x(n)^3) 

这样我们就得到了迭代计算的公式。也就是源程序中的

 x=x*(1.5f-xhalf*x*x); 

 

以下关键的一步就是计算x(n)的一个真根足够靠近的近似值。

主要思想利用了浮点数的表示法，

ieee754标准当中关键部分是

 long i=*(long*)&x; 

 i=0x5f3759df - (i>>1); 

x=*(float *)&i; 

其它地方都比較简单，这三句的意思，应当是把x的浮点表示格式直接看出long类型的，

然后进行0x5f3759df - (i>>1)运算。之后再直接把long类型的看成float类型的。当中的数学道理，我大致看了下

对于浮点数的表示，在ieee里float类型是依照符号指数.尾数格式进行的。当中符号位1位，指数位8位，尾数23位。先仅仅考虑指数部分。

 

long i=*(long*)&x;是把浮点数看成了与他具有同样位格式的long型数。

i=0x5f3759df - (i>>1);

对于0x5f3759df,假设看成浮点数的话，二进制位为

.10111110.01101110101100111011111,

其指数位为10111110=190 。

 

设原来的x的IEEE表示的指数位为a则新的结果的指数位变成了190-a/2；注意假设要转换成

m*2^n格式，须要将指数减去127；所以得到的新结果用数学表示的指数位为190-a/2-127= 63-a/2 。而x的原来的指数a，假设从ieee格式变成数学表示应为a-127 。

 

比較63-a/2与a-127，能够看出指数之间的关系为63-a/2==（a-127)*(-1/2) 

所以经过上述三句。x应该已经是1/sqrt(x)的近似值。

当然上面仅仅是从指数位进行了简单分析，假设要细化到尾数，还须要看详细效果。

 

 

 

 

在3D图形编程中，常常要求平方根或平方根的倒数。比如：求向量的长度或将向量归一化。C数学函数库中的sqrt具有理想的精度，但对于3D游戏程式来说速度太慢。我们希望可以在保证足够的精度的同一时候。进一步提快速度。

Carmack在QUAKE3中使用了以下的算法，它第一次在公众场合出现的时候，差点儿震住了全部的人。

据说该算法事实上并非Carmack发明的，它真正的作者是Nvidia的Gary Tarolli（未经证实）。

 

//
// 计算參数x的平方根的倒数
//
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;
}

该算法的本质事实上就是牛顿迭代法（Newton-RaphsonMethod，简称NR）。而NR的基础则是泰勒级数（Taylor Series）。

NR是一种求方程的近似根的方法。首先要预计一个与方程的根比較靠近的数值，然后依据公式推算下一个更加近似的数值，不断反复直到能够获 得惬意的精度。其公式例如以下：

函数：y=f(x),其一阶导数为：y'=f'(x)

则方程：f(x)=0 的第n+1个近似根为

x[n+1] = x[n] - f(x[n]) / f'(x[n])NR最关键的地方在于预计第一个近似根。假设该近似根与真根足够靠近的话，那么仅仅须要少数几次迭代，就能够得到惬意的解。

如今回过头来看看怎样利用牛顿法来解决我们的问题。求平方根的倒数。实际就是求方程1/(x^2)-a=0的解。将该方程按牛顿迭代法的公式展开为：

x[n+1]=1/2*x[n]*(3-a*x[n]*x[n])将1/2放到括号中面，就得到了上面那个函数的倒数第二行。

接着，我们要设法预计第一个近似根。这也是上面的函数最奇妙的地方。

它通过某种方法算出了一个与真根很接近的近似根，因此它仅仅须要使用一次迭代过程就获得了较惬意的解。它是如何做到的呢？全部的奥妙就在于这一行：

i = 0x5f3759df - (i >> 1); // 计算第一个近似根超级莫名其妙的语句。不是吗？但细致想一下的话，还是能够理解的。

我们知道。IEEE标准下，float类型的数据在32位系统上是这样 表示的（大体来说就是这样，但省略了非常多细节，有兴趣能够GOOGLE）：

bits：31 30 ... 0
31：符号位
30-23：共8位，保存指数（E）
22-0：共23位，保存尾数（M）所 以，32位的浮点数用十进制实数表示就是：M*2^E。开根然后倒数就是：M^(-1/2)*2^(-E/2)。

如今就十分清晰了。

语句 i>>1其工作就是将指数除以2，实现2^(E/2)的部分。而前面用一个常数减去它，目的就是得到M^(1/2)同一时候反转全部指数的符号。

至于那个0x5f3759df。呃，我仅仅能说。的确是一个超级的Magic Number。

那个Magic Number是能够推导出来的，但我并不打算在这里讨论，由于实在太繁琐了。简单来说。其原理例如以下：由于IEEE的浮点数中。尾数M省略了最前面的1，所 以实际的尾数是1+M。假设你在大学上数学课没有打瞌睡的话。那么当你看到(1+M)^(-1/2)这种形式时，应该会立即联想的到它的泰勒级数展开。 而该展开式的第一项就是常数。以下给出简单的推导过程：

对于实数R>0。如果其在IEEE的浮点表示中。
指数为E，尾数为M，则：R^(-1/2) = (1+M)^(-1/2) * 2^(-E/2)

 将(1+M)^(-1/2)按泰勒级数展开，取第一项。得：

原式= (1-M/2) * 2^(-E/2)=2^(-E/2) - (M/2) * 2^(-E/2)

假设不考虑指数的符号的话。
(M/2)*2^(E/2)正是(R>>1)。
而在IEEE表示中，指数的符号仅仅需简单地加上一个偏移就可以，而式子的前半部分刚好是个常数，所以原式能够转化为：

原式 = C - (M/2)*2^(E/2) = C- (R>>1)，当中C为常数

所以仅仅须要解方程：
R^(-1/2)
= (1+M)^(-1/2) * 2^(-E/2)
= C - (R>>1)
求出令到相对误差最小的C值就能够了上面的推导过程仅仅是我个人的理解，并未得到证实。

而Chris Lomont则在他的论文中具体讨论了最后那个方程的解法。并尝试在实际的机器上寻找最佳的常数C。有兴趣的朋友能够在文末找到他的论文的链接。

所以，所谓的Magic Number，并非从N元宇宙的某个星系因为时空扭曲而掉到地球上的，而是几百年前就有的数学理论。

仅仅要熟悉NR和泰勒级数。你我相同有能力作出类似的优化。

在http://GameDev.net上有人做过測试，该函数的相对误差约为0.177585%，速度比C标准库的sqrt提高超过20%。

假设添加一次迭代过 程，相对误差能够减少到e-004的级数，但速度也会降到和sqrt差点儿相同。

据说在DOOM3中，Carmack通过查找表进一步优化了该算法。精度近乎完美。并且速度也比原版提高了一截（正在努力弄源代码。谁有发我一份）。

值得注意的是，在Chris Lomont的演算中，理论上最棒的常数（精度最高）是0x5f37642f。而且在实际測试中，假设仅仅使用一次迭代的话，其效果也是最好的。

但奇怪的 是。经过两次NR后。在该常数下解的精度将减少得很厉害（天知道是怎么回事！）。经过实际的測试，Chris Lomont觉得，最棒的常数是0x5f375a86。

假设换成64位的double版本号的话，算法还是一样的。而最优常数则为 0x5fe6ec85e7de30da（又一个令人冒汗的Magic Number- -b）。

这个算法依赖于浮点数的内部表示和字节顺序。所以是不具移植性的。假设放到Mac上跑就会挂掉。

假设想具备可移植性，还是乖乖用sqrt好了。但算法思想是通用的。大家庭可以尝试对相应的平方根算法预测。

 

版权声明：本文【借给你1秒】原创文章，转载请注明出处。