sqrt函数实现

        转自：http://blog.csdn.net/gqtcgq/article/details/48392255

        Implement int sqrt(int x).

        Compute and return the square root of x.

 

1：二分查找

        思路：要实现一个sqrt函数，可以使用二分法，首先确定一个范围[begin, end]，这个范围的中间数mid，看mid的平方是否等于x，如果相等，则返回mid，如果不等则缩小[begin,end]的范围，为原来的一半。这里的初始范围可以是[1, x]，也可以是更精确一些的[1, (x/2) + 1]。（因 (x/2) + 1 的平方等于 x+1+(x^2/4)，它一定大于x，所以，x的平方根一定在[1, (x/2) + 1]范围内）

 

        题目中给出的函数原型是int mySqrt(int x)。参数和返回值都是整数。这里稍微扩展一下，将函数原型改为double mySqrt(int x)。解题思路还是一样的，但是浮点数因精度的原因，无法判断两个浮点数是否完全相等，只能说两者的差值绝对值小于某个精度，所以在二分查找时，需要一定的技巧，具体的代码如下：

[cpp]  view plain  copy

1. double mySqrt_binarysearch(int x)   
2. {  
3.     if(x <= 0)   return 0;  
4.   
5.     double begin = 1;  
6.     double end = x/2+1;  
7.     double mid, lastmid;  
8.   
9.     mid = begin + (end-begin)/2;  
10.     do{  
11.         if(mid < x/mid) begin = mid;  
12.         else    end = mid;  
13.   
14.         lastmid = mid;  
15.         mid = begin + (end-begin)/2;  
16.     }  
17.     while(ABS(lastmid-mid) > FLT_MIN);  
18.       
19.     return mid;  
20. }  

        上面的代码中，逐步缩小[begin,end]的范围，通过判断上次的lastmid与本次的mid的差值绝对值是否在精度之内，来决定是否继续寻找下去。

 

2：牛顿迭代法

        上面的实现方法只能说是中规中矩，但是实现sqrt有更牛逼的方法，就是牛顿迭代法。该方法就是由我们熟知的牛顿提出的。具体思想可以自行搜索。简而言之，如下图：

 

        x^2 = a的解，也就是函数f(x) = x^2 – a与x轴的交点。可以在x轴上先任选一点x0，则点（x0, f(x0)）在f(x)上的切线，与x轴的交点为x1，它们满足切线的方程：f(x0)=(x0-x1)f’(x0)，可得x1更接近最终的结果，解方程得到：

x1 = (x0 + (a/x0))/2。以x1为新的x0，按照切线的方法依次迭代下去，最终求得符合精确度要求的结果值。它的实现代码如下：

[cpp]  view plain  copy

1. double mySqrt_newton(int x)   
2. {  
3.     if(x <= 0)   return 0;  
4.   
5.     double res, lastres;  
6.   
7.     res = x;    //初始值，可以为任意非0的值  
8.       
9.     do{  
10.         lastres = res;  
11.         res = (res + x/res)/2;  
12.     }  
13.     while(ABS(lastres-res) > FLT_MIN);  
14.   
15.     return res;  
16. }  

       使用牛顿法解决sqrt的效率非常高，关于效率比较可参见本文最后一节。牛顿法的效率很大程度上取决于初始值的选取，这就引出了下一节。

 

3：神迹

       下面这段代码出自《雷神之锤》，至今尚未找到该代码的真正作者，代码如下：

[cpp]  view plain  copy

1. float InvSqrt(float x)  
2. {  
3.     float xhalf = 0.5f * x;  
4.     int i = *(int*)&x;   
5.     i = 0x5f375a86 - (i>>1);   
6.     x = *(float*)&i;  
7.     x = x*(1.5f-xhalf*x*x);   
8.     x = x*(1.5f-xhalf*x*x);   
9.     x = x*(1.5f-xhalf*x*x);  
10.   
11.     return 1/x;  
12. }  

       它本质上还是使用的牛顿迭代法，真正牛逼的地方在于它初始值的选择，0x5f375a86这个魔法数字的由来尚不可知，该算法的具体原理及其背景可以参见维基百科，不再赘述。

       要注意的是，上面算法使用的是float和int类型，实验可知他们不能替换为double和long。

 

4：效率

       使用下面的代码，测试上述三种方法，以及系统默认方法的效率：

[cpp]  view plain  copy

1. int main(int argc, char **argv)  
2. {  
3.     clock_t begin, end;  
4.     int num = atoi(argv[1]);  
5.     double res;  
6.       
7.     int i;  
8.     int loopcnts = 1000000;  
9.       
10.   
11.     begin = clock();  
12.     for(i = 0; i < loopcnts; i++)  
13.         res = mySqrt_binarysearch(num);  
14.     end = clock();  
15.     printf("mySqrt_binarysearch(%d) = %f, spent time is %f\n", num, res, (double)(end-begin));  
16.   
17.   
18.     begin = clock();  
19.     for(i = 0; i < loopcnts; i++)  
20.         res = mySqrt_newton(num);  
21.     end = clock();  
22.     printf("mySqrt_newton(%d) = %f, spent time is %f\n", num, res, (double)(end-begin));  
23.   
24.   
25.     begin = clock();  
26.     for(i = 0; i < loopcnts; i++)  
27.         res = InvSqrt(num);  
28.     end = clock();  
29.     printf("InvSqrt(%d) = %f, spent time is %f\n", num, res, (double)(end-begin));  
30.   
31.   
32.   
33.     begin = clock();  
34.     for(i = 0; i < loopcnts; i++)  
35.         res = sqrt(num);  
36.     end = clock();  
37.     printf("system sqrt(%d) = %f, spent time is %f\n", num, res, (double)(end-begin));  
38.   
39. }  

 

       测试结果如下：

[cpp]  view plain  copy

1. mySqrt_binarysearch(65535) = 255.998047, spent time is 3437535.000000  
2. mySqrt_newton(65535) = 255.998047, spent time is 659694.000000  
3. InvSqrt(65535) = 255.998047, spent time is 65902.000000  
4. system sqrt(65535) = 255.998047, spent time is 82605.000000  

       可见，二分法最慢，普通的牛顿迭代法次之，神迹代码要比系统库的还要快一些。