一、牛顿迭代法
以下算法在GameDev.net 上有人做过测试,该函数的相对误差约为0.177585%,速度比C标准库的sqrt提高超过20%。如果增加一次迭代过程,相对误差可以降低到e- 004 的级数,但速度也会降到和sqrt差不多。
/* 计算参数x的平方根的倒数
利用牛顿法来解决求平方根的倒数,实际就是求方程1/(x^2)-a=0的解。将该方程按牛顿迭代法的公式展开为:
x[n+1]=1/2*x[n]*(3-a*x[n]*x[n])
*/
float InvSqrt (float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f3759df - (i >> 1); // 计算第一个近似根
x = *(float*)&i;
x = x*(1.5f - xhalf*x*x); // 牛顿迭代法
return x;
}
本质其实就是牛顿迭代法(Newton-Raphson Method,简称NR),而NR的基础则是泰勒级数(Taylor Series)。 NR是一种求方程的近似根的方法。首先要估计一个与方程的根比较靠近的数值,然后根据公式推算下一个更加近似的数值,不断重复直到可以获得满意的精度。NR最关键的地方在于估计第一个近似根。如果该近似根与真根足够靠近的话,那么只需要少数几次迭代,就可以得到满意的解。
公式如下:
函数:y=f(x)
其一阶导数为:y'=f'(x)
则方程:f(x)=0 的第n+1个近似根为
x[n+1] = x[n] - f(x[n]) / f'(x[n])
注意:
1.0x5f3759df 是泰勒级数展开式的第一项 ,理论上最优秀的常数(精度最高)是 0x5f37642f
2.如果换成64位的double版本的话, 算法还是一样的,而最优常数则为0x5fe6ec85e7de30da
二阶展开
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 );
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed ,
//对于上面这句可以屏蔽,不需要那么高的精度,反倒是速度要紧
#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}
double InvSqrt(double number)
{
__int64 i;
double x2, y;
const double threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( __int64 * ) &y;
i = 0x5fe6ec85e7de30da - ( i >> 1 );
y = * ( double * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
return y;
}
参考:https://www.cnblogs.com/topameng/archive/2010/10/08/1845852.html
二、方法2
转载:https://wenku.baidu.com/view/a0f7ad255f0e7cd18525361d.html