算法笔记--牛顿迭代法

牛顿迭代法：

牛顿迭代法（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
",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