计算平方根的算法

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

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 (&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

 

posted @ 2012-05-15 22:34  Mr.Rico  阅读(...)  评论(... 编辑 收藏