申明,本文非笔者原创,原文转载自:http://www.cnblogs.com/xkfz007/archive/2012/05/15/2502348.html
总结一下一些常用的计算平方根的方法
1. 牛顿法
具体的做法如下:
具体的计算程序如下:
double sqrt_(
double x)
{
double g=x;
while(ABS(g*g-x)> 0.000001)
{
g=(g+x/g)/ 2;
}
return g;
}
{
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;
}
{
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])
则方程: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])
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;
}
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
结果如下:
第一列是随机数,第二列是库函数sqrt, <>中的是时间,y1是方法1,y2是方法2,y3是方法3
看来这三种方法时间都挺快,都小于毫秒
参考:http://www.cnblogs.com/vagerent/archive/2007/06/25/794695.html