标签：c++   平方根   平方根倒数   

我们平时经常会有一些数据运算的操作，需要调用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 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））代表的是小数还是整数。以上图为例，将描述带入有），且，则可得其表示的浮点数为。

8位二进制整数补码示例

如上所述，一个有符号正整数在二进制补码系统中的表示中首位为0，而后面的各位则用于表示其数值。将浮点数取别名存储为整数时，该整数的数值即为，其中E表示指数，M表示有效数字；若以上图为例，图中样例若作为浮点数看待有，，则易知其转化而得的整数型号数值为。由于平方根倒数函数仅能处理正数，因此浮点数的符号位（即如上的Si）必为0，而这就保证了转换所得的有符号整数也必为正数。以上转换就为后面的计算带来了可行性，之后的第一步操作（逻辑右移一位）即是使该数的长整形式被2所除。^[3] 

4历史与考究

《雷神之锤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好了。但算法思想是通用的。大家可以尝试推算一下相应的平方根算法。

一个Sqrt函数引发的血案

标签：c++   平方根   平方根倒数   

原文地址：http://blog.csdn.net/u013467442/article/details/41147487