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


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

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

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

    代码:

     1 #include<bits/stdc++.h>
     2 #define MAXN 100000+10
     3 #define MAX 100000000
     4 #define eps 1e-6
     5 #define ll long long
     6 using namespace std;
     7 
     8 float get_sqrt(float x)
     9 {
    10     float low=0, up=x, mid, now;
    11     mid=(low+up)/2;
    12     do
    13     {
    14         now=mid;        //**now保存上一次计算的值
    15         if(mid*mid<x)   //**mid偏小,右移
    16         {
    17             low=mid;
    18         }
    19         else       //**mid偏大,左移
    20         {
    21             up=mid;
    22         }
    23         mid=(low+up)/2;
    24     }while(abs(mid-now)>eps); //**两次计算的误差小于eps,mid即为所求值
    25     return mid;
    26 }
    27 
    28 int main(void)
    29 {
    30     std::ios::sync_with_stdio(false), cin.tie(0), cin.tie(0);
    31     float x;
    32     cin >> x;
    33     cout << get_sqrt(x) << endl;
    34     return 0;
    35 }



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

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

    假设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;

    代码:

     1 #include<bits/stdc++.h>
     2 #define MAXN 100000+10
     3 #define MAX 100000000
     4 #define eps 1e-6
     5 #define ll long long
     6 using namespace std;
     7 
     8 float get_sqrt(float a)
     9 {
    10     float x, now;
    11     x=a;
    12     do
    13     {
    14         now=x;     //**now保存上一次的x值
    15         x=(x+a/x)/2;  //**通过迭代更新x的值使其接近解
    16     }while(abs(now-x)>eps); //**两次计算的误差小于eps,x即为所求值
    17     return x;
    18 }
    19 
    20 int main(void)
    21 {
    22     std::ios::sync_with_stdio(false), cin.tie(0), cin.tie(0);
    23     float x;
    24     cin >> x;
    25     cout << get_sqrt(x) << endl;
    26     return 0;
    27 }



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

    神代码(效率为库函数的4倍):

     1 #include<bits/stdc++.h>
     2 #define MAXN 100000+10
     3 #define MAX 100000000
     4 #define eps 1e-6
     5 #define ll long long
     6 using namespace std;
     7 
     8 /*float InvSqrt( float number )
     9 {
    10     long i;
    11     float x2, y;
    12     const float threehalfs = 1.5F;
    13     x2 = number * 0.5F;
    14     y   = number;
    15     i   = * ( long * ) &y;   // evil floating point bit level hacking
    16     i   = 0x5f3759df - ( i >> 1 ); // what the fuck?
    17     y   = * ( float * ) &i;
    18     y   = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
    19     // y   = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed //**可以通过减少迭代次数来用精度换取时间
    20     #ifndef Q3_VM
    21     #ifdef __linux__
    22          assert( !isnan(y) ); // bk010122 - FPE?
    23     #endif
    24     #endif
    25     return 1/y;
    26 }*/
    27 
    28 float InvSqrt(float x)
    29 {
    30     float xhalf = 0.5f*x;
    31     int i = *(int*)&x; // get bits for floating VALUE
    32     i = 0x5f375a86- (i>>1); // gives initial guess y0
    33     x = *(float*)&i; // convert bits BACK to float
    34     x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
    35     x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
    36 //    x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy  //**可以通过减少迭代次数来用精度换取时间
    37     return 1/x;
    38 }
    39 
    40 int main(void)
    41 {
    42     std::ios::sync_with_stdio(false), cin.tie(0), cin.tie(0);
    43     float x;
    44     cin >> x;
    45     cout << InvSqrt(x) << endl;
    46     return 0;
    47 }



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

    然而我并没有看懂,待以后继续研究;

  • 相关阅读:
    Jenkins安装2
    11月14 vagrant 虚拟机连接
    11月13 ssh 和 netstat
    Jenkins安装
    11月12号 springboot1.5 引入redis
    11月12号 用户登录输入密码错误达到指定次数后,锁定账户 004
    js模式-观察者模式
    立即执行函数
    vue 获取数据
    在登陆退出时候使用Vuex
  • 原文地址:https://www.cnblogs.com/geloutingyu/p/5882716.html
Copyright © 2020-2023  润新知