总结一下一些常用的计算平方根的方法
1. 牛顿法
具体的做法如下：
具体的计算程序如下：
double sqrt_( double x) { double g=x; while(ABS(g*g-x)> 0.000001) { g=(g+x/g)/ 2; } return g; }
2. 利用级数进行逼近
微积分中的泰勒级数如下：
这样我们可以进行在一定精度内的逼近。
但是这儿存在一个问题，就是这个公式的收敛问题。它是存在收敛区间的。
所以可以得到最后的代码：
double Tsqrt( double x) // 计算[0,2)范围内数的平方根 { double sum,coffe,factorial,xpower,term; int i; sum= 0; coffe= 1; factorial= 1; xpower= 1; term= 1; i= 0; while(ABS(term)> 0.000001) { sum+=term; coffe*=( 0.5-i); factorial*=(i+ 1); xpower*=(x- 1); term=coffe*xpower/factorial; i++; } return sum; } double sqrt2( double x) // 大于2的数要转化为[0,2)区间上去 { double correction= 1; while(x>= 2) { x/= 4; correction*= 2; } return Tsqrt(x)*correction; }
3. 平方根倒数速算法
这是牛顿迭代法的应用。
牛顿迭代法是一种求方程的近似根的方法。首先要估计一个与方程的根比较靠近的数值，然后根据公式推算下一个更加近似的数值，不断重复直到可以获得满意的精度。其公式如下：
函数：y=f(x)
其一阶导数为：y'=f'(x) 则方程：f(x)=0 的第n+1个近似根为 x[n+1] = x[n] - f(x[n]) / f'(x[n])
牛顿迭代法最关键的地方在于估计第一个近似根。如果该近似根与真根足够靠近的话，那么只需要少数几次迭代，就可以得到满意的解。 现在回过头来看看如何利用牛顿法来解决我们的问题。求平方根的倒数，实际就是求方程1/(x^2)-a=0的解。将该方程按牛顿迭代法的公式展开为： x[n+1]=1/2*x[n]*(3-a*x[n]*x[n])
这个方法的一个关键地方还在于起始值的选取。
具体代码如下：
float SquareRootFloat( float number) { long i; float x, y; const float f = 1.5F; x = number * 0.5F; y = number; i = * ( long * ) &y; i = 0x5f3759df - ( i >> 1 ); // 注意这一行 y = * ( float * ) &i; y = y* ( f - ( x * y * y ) ); y = y * (f - ( x * y * y )); return number * y; }
该方法又叫卡马克反转。其中的0x5f3759df的来历比较复杂，在维基百科中也有介绍。最后的参考中也有介绍。
最后三种方法的测试代码如下：
View Code #include <stdio.h> #include <stdlib.h> #include <time.h> #include <math.h> #include <sys/timeb.h> #define TIMEB timeb #define ABS(x) ((x)>0?(x):-(x)) #define eps 0.0000000001 time_t time_test( double x, double *y, double(*func)( double)) { time_t ltime1, ltime2, tmp_time; struct TIMEB tstruct1, tstruct2; ftime (&tstruct1); // start time ms time (<ime1); // start time s *y=func(x); time (<ime2); // end time sec ftime (&tstruct2); // end time ms tmp_time = (ltime2 * 1000 + tstruct2.millitm) - (ltime1 * 1000 + tstruct1.millitm); return tmp_time; } double sqrt1( double x) { double g=x; while(ABS(g*g-x)>eps) { g=(g+x/g)/ 2; } return g; } double Tsqrt( double x) { double sum,coffe,factorial,xpower,term; int i; sum= 0; coffe= 1; factorial= 1; xpower= 1; term= 1; i= 0; while(ABS(term)>eps) { sum+=term; coffe*=( 0.5-i); factorial*=(i+ 1); xpower*=(x- 1); term=coffe*xpower/factorial; i++; } return sum; } double sqrt2( double x) { double correction= 1; while(x>= 2) { x/= 4; correction*= 2; } return Tsqrt(x)*correction; } double sqrt3( double xi) { long i; float number=( float)xi; float x, y; const float f = 1.5F; x = number * 0.5F; y = number; i = * ( long * ) &y; i = 0x5f3759df - ( i >> 1 ); // 注意这一行 y = * ( float* ) &i; y = y* ( f - ( x * y * y ) ); y = y * (f - ( x * y * y )); return number * y; } int main(){ int x; double y,y1,y2,y3; int i; time_t t,t1,t2,t3; srand(( int)time(NULL)); for(i= 0;i< 20;i++) { x=rand()% 1000; t=time_test(x,&y,sqrt); t1=time_test(x,&y1,sqrt1); t2=time_test(x,&y2,sqrt2); t3=time_test(x* 1.0,&y3,sqrt3); printf( " x=%3d y=%7.4f<%2d> y1=%7.4f<%2d> y2=%7.4f<%2d> y3=%7.4f<%2d>\n ",x,y,t,y1,t1,y2,t2,y3,t3); } return 0; }
结果如下：
第一列是随机数，第二列是库函数sqrt， <>中的是时间，y1是方法1，y2是方法2，y3是方法3
看来这三种方法时间都挺快，都小于毫秒
参考：