算法筆記--牛頓迭代法

牛頓迭代法：

牛頓迭代法（Newton's method）又稱為牛頓-拉夫遜（拉弗森）方法（Newton-Raphson method），它是牛頓在17世紀提出的一種在實數域和復數域上近似求解方程的方法。

詳見：https://www.matongxue.com/madocs/205.html#/madoc

設x*為f(x) = 0 的根

計算公式（迭代公式）：x_n+1 = x_n - f(x_n) / f'(x_n) 

帶入一個初始點x₀ 即可啟動迭代

x_n ->x* (n -> ∞)

 

---

 

牛頓迭代法在acm中的應用：

1.求解方程的根

HUD 2899

枚舉初始迭代點求f'(x) = 0 的根

#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define pb push_back
#define mem(a, b) memset(a, b, sizeof(a))

const double eps = 1e-9;
double y;
double f(double x) {
    return 6 * x*x*x*x*x*x*x + 8 * x*x*x*x*x*x + 7 * x*x*x + 5 * x*x - y * x;
}
double _f(double x) {
    return 42 * x*x*x*x*x*x + 48 * x*x*x*x*x + 21 * x*x + 10 * x - y;
}
double __f(double x) {
    return 42 * 6 * x*x*x*x*x + 48 * 5 * x*x*x*x + 21 * 2 * x +10;
}
double newton(double x) {
    int tot = 0;
    while(fabs(_f(x))>eps) {
        x -= _f(x) / __f(x);
    }
    return x;
}
int main() {
    int T;
    scanf("%d", &T);
    while (T--) {
        scanf("%lf", &y);
        double ans = min(f(0), f(100));
        for (int i = 0; i <= 100; i++) {
            double x = newton(i);
            if(0<= x && x <= 100) {
                ans = min(ans, f(x));
            }
        }
        printf("%.4lf\n",ans);
    }
    return 0;
}

2.開根號

牛頓迭代法相比於二分法，迭代次數更少。

而且在雷神之錘3的源碼中還有一個更強大的求1/sqrt(a)的方法，比c++ 標准庫快好幾倍

其源碼如下所示

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 ); // 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 y;
}  

關於注釋為"what the fuck?"的一行，這里有一篇論文

http://www.matrix67.com/data/InvSqrt.pdf

求解1/sqrt(a) = x

1/(x^2) - a = 0

令f(x) = 1/(x^2) - a

f'(x) = -2/(x^3)

x_n+1 = x_n - f(x_n) / f'(x_n) = x*(1.5 -0.5a*x^2)

參考雷聲之錘源碼，我們得到一個速度更優的求sqrt(a)的方法

double Sqrt(float x){
    double xhalf = 0.5*(double)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.0/(double)x;
}

 參考：

https://www.cnblogs.com/ECJTUACM-873284962/p/6536576.html

https://blog.csdn.net/kenden23/article/details/14543521

https://blog.csdn.net/wubaizhe/article/details/75574798