计算平方根的算法

总结一下一些常用的计算平方根的方法

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的来历比较复杂，在维基百科中也有介绍。最后的参考中也有介绍。

最后三种方法的测试代码如下：

#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 (&ltime1);               // start time s

    *y=func(x);

    time (&ltime2);               // 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

看来这三种方法时间都挺快，都小于毫秒

参考：http://www.cnblogs.com/vagerent/archive/2007/06/25/794695.html