求sqrt()底层效率问题（二分/牛顿迭代）

偶然看见一段求根的神代码，于是就有了这篇博客；

对于求根问题，通常我们可以调用sqrt库函数，不过知其然需知其所以然，我们看一下求根的方法；

比较简单方法就是二分咯：

代码：

#include<bits/stdc++.h>
#define MAXN 100000+10
#define MAX 100000000
#define eps 1e-6
#define ll long long
using namespace std;

float get_sqrt(float x)
{
    float low=0, up=x, mid, now;
    mid=(low+up)/2;
    do
    {
        now=mid;        //**now保存上一次计算的值
        if(mid*mid<x)   //**mid偏小，右移
        {
            low=mid;
        }
        else       //**mid偏大，左移
        {
            up=mid;
        }
        mid=(low+up)/2;
    }while(abs(mid-now)>eps); //**两次计算的误差小于eps，mid即为所求值
    return mid;
}

int main(void)
{
    std::ios::sync_with_stdio(false), cin.tie(0), cin.tie(0);
    float x;
    cin >> x;
    cout << get_sqrt(x) << endl;
    return 0;
}

 

然而虽然其计算的结果和库函数一样，然而其效率较之库函数差数百倍，当然不是我说的神代码咯；

 

我们再看一下牛顿迭代法如何；

假设a为需求根的数，x为其正根，则有a=x*x；即a的正根为函数f(x)=x*x-a与x轴的正交点；

由牛顿迭代法我们可以知道，可以通过（x，f(x)）的切线不断逼近解；

任取一点（x0，f(x0)），其切线方程为 y=f'(x0)*(x-x0)+f(x0)，其与x轴的交点为x1=x0-f(x0)/f'(x0)，x1是一个比x0更接近的近似解；

依次类推，可以求出x2，x2又比x1更接近；

可以求出x3......

由此我们得出迭代公式为：x'=x-f(x)/f'(x)，再带入f(x)=x*x-a得：x'=(x+a/x)/2；

 

代码：

#include<bits/stdc++.h>
#define MAXN 100000+10
#define MAX 100000000
#define eps 1e-6
#define ll long long
using namespace std;

float get_sqrt(float a)
{
    float x, now;
    x=a;
    do
    {
        now=x;     //**now保存上一次的x值
        x=(x+a/x)/2;  //**通过迭代更新x的值使其接近解
    }while(abs(now-x)>eps); //**两次计算的误差小于eps，x即为所求值
    return x;
}

int main(void)
{
    std::ios::sync_with_stdio(false), cin.tie(0), cin.tie(0);
    float x;
    cin >> x;
    cout << get_sqrt(x) << endl;
    return 0;
}

 

其效率较之二分高了很多，不过还是不如库函数；

 

神代码（效率为库函数的4倍）：

#include<bits/stdc++.h>
#define MAXN 100000+10
#define MAX 100000000
#define eps 1e-6
#define ll long long
using namespace std;

/*float InvSqrt( 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 ); // what the fuck?
    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 1/y;
}*/

float InvSqrt(float x)
{
    float xhalf = 0.5f*x;
    int i = *(int*)&x; // get bits for floating VALUE
    i = 0x5f375a86- (i>>1); // gives initial guess y0
    x = *(float*)&i; // convert bits BACK to float
    x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
    x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
//    x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy  //**可以通过减少迭代次数来用精度换取时间
    return 1/x;
}

int main(void)
{
    std::ios::sync_with_stdio(false), cin.tie(0), cin.tie(0);
    float x;
    cin >> x;
    cout << InvSqrt(x) << endl;
    return 0;
}

 

其本质还是迭代法，不过因为那个魔鬼般的常数0x5f375a86 和 i = 0x5f375a86- (i>>1);中的位运算大大提高了其速度；

然而我并没有看懂，待以后继续研究；