偶然看见一段求根的神代码,于是就有了这篇博客;
对于求根问题,通常我们可以调用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);中的位运算大大提高了其速度;
然而我并没有看懂,待以后继续研究;