开平方实现

首先是最普通的CRT里自带的sqrt，只需要引用math.h就可以使用了：

[cpp] view plain copy

 print?

1. #include <math.h>  
2. result = sqrt(number);  

 

 

接下来是传统的牛顿迭代法，我们计算开方的时候就是手工不断尝试每一位最合适的数字，然后一步步收敛求得更精确的答案。牛顿迭代法就是通过折半来快速收敛，每迭代一次就得到更精确的结果，下面这段程序采用的是迭代10次的牛顿迭代法：

[cpp] view plain copy

 print?

1. float newtonSqrt(float a){  
2.   int i;  
3.   float x;  
4.   x=a;  
5.   for(i=1;i<=10;i++)  
6.     x=(x+a/x)/2;  
7.   return x;  
8. }  

 

 

 

既然每次计算都需要那么多步骤，为什么不预先计算好之后储存，使用的时候再通过查表来取得结果呢？

首先是建立一个表，把计算好的结果都存放进去，我这里建立的是0到10，精确到百分位总共有1001个条目的表：

[cpp] view plain copy

 print?

1. float *makeSqrtTable()  
2. {  
3.   int num;  
4.   float *tbl = (float *)malloc(sizeof(float) * (1000 + 1));  
5.   for (num = 0; num <= 1000; num += 1)  
6.   {  
7.     tbl[num] = sqrt((float)num / 100);  
8.   }  
9.   return tbl;  
10. }  

 

不需要用表的时候可以把表释放掉：

[cpp] view plain copy

 print?

1. free(tbl);  

 

然后是查表的函数，直接从对应的项目取出来就行了：

[cpp] view plain copy

 print?

1. float tabelSqrt(float *tbl, float num)  
2. {  
3.   return tbl[(int)(num * 100)];  
4. }  

 

 

最后，是出现在Quake3代码里的，顶级游戏引擎设计师约翰卡马克的开方方法，quake3-1.32b/code/game/q_math.c里可以找到实现：

[cpp] view plain copy

 print?

1. //============================================================================  
2. #if !idppc  
3. /* 
4. ** float q_rsqrt( float number ) 
5. */  
6. float Q_rsqrt( float number )  
7. {  
8.     long i;  
9.     float x2, y;  
10.     const float threehalfs = 1.5F;  
11.     x2 = number * 0.5F;  
12.     y  = number;  
13.     i  = * ( long * ) &y;                       // evil floating point bit level hacking  
14.     i  = 0x5f3759df - ( i >> 1 );               // what the f**k?  
15.     y  = * ( float * ) &i;  
16.     y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration  
17. //  y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed  
18. #ifndef Q3_VM  
19. #ifdef __linux__  
20.     assert( !isnan(y) ); // bk010122 - FPE?  
21. #endif  
22. #endif  
23.     return y;  
24. }  

 

还有另一个版本的实现是：

 

[cpp] view plain copy

 print?

1. float carmSqrt(float x) {  
2.   union{  
3.     int intPart;  
4.     float floatPart;  
5.   } convertor;  
6.   union{  
7.     int intPart;  
8.     float floatPart;  
9.   }convertor2;  
10.   convertor.floatPart = x;  
11.   convertor2.floatPart = x;  
12.   convertor.intPart = 0x1FBCF800 + (convertor.intPart >> 1);  
13.   convertor2.intPart = 0x5f3759df - (convertor2.intPart >> 1);  
14.   return 0.5f * (convertor.floatPart + (x * convertor2.floatPart));  
15. }  

 

 

卡马克进行开方时也进行了迭代，但是只迭代了一次，因为通过那个奇怪的数字0x5f3759df就使得“只用迭代一次就取到比较满意的结果”。

测试几个数开方：

4

Sqrt:2.000000

NewtonSqrt:2.000000

TabelSqrt:2.000000

CarmSqrt:1.954374

10

Sqrt:3.162278

NewtonSqrt:3.162278

TabelSqrt:3.162278

CarmSqrt:3.235606

除了卡马克开方之外，其他的方法得到的结果都是一样的，也是更精确的，因为卡马克开方虽然快速逼近最终值，但是只迭代了一次。

 

接下来是评测一下速度，我们从0.00到10.00每次递增0.01共1001个数字进行开方，统计时间，在我这台机器得到的结果如下，单位是毫秒：

Sqrt:0:71

Carm:0:28

Newton:0:472

Tabel:0:52

可以看到，卡马克开方惊人的快，甚至差不多比查表要少一半的时间。

如果在精度要求允许的范围，用卡马克的方法是最好的,要更精确的结果，可以多迭代几次，这都没问题；如果要实现任意精度开方，可以用牛顿迭代法计算；如果需要精度和速度之间做平衡，同时不在乎占用的空间，可以使用查表法；当然，CRT里自带的sqrt就已经做得相当好了，一般使用没什么问题。