自学教程：详解Python牛顿插值法

/* * @Author: csc * @Date: 2021-04-30 08:52:45 * @LastEditTime: 2021-04-30 11:57:46 * @LastEditors: Please set LastEditors * @Description: In User Settings Edit * @FilePath: /code_formal/course/cal/newton_quo.cpp */#include <bits/stdc++.h>#define pr printf#define sc scanf#define for0(i, n) for (i = 0; i < n; i++)#define for1n(i, n) for (i = 1; i <= n; i++)#define forab(i, a, b) for (i = a; i <= b; i++)#define forba(i, a, b) for (i = b; i >= a; i--)#define pb push_back#define eb emplace_back#define fi first#define se second#define int long long#define pii pair<int, int>#define vi vector<int>#define vii vector<vector<int>>#define vt3 vector<tuple<int, int, int>>#define mem(ara, n) memset(ara, n, sizeof(ara))#define memb(ara) memset(ara, false, sizeof(ara))#define all(x) (x).begin(), (x).end()#define sq(x) ((x) * (x))#define sz(x) x.size()const int N = 2e5 + 100;const int mod = 1e9 + 7;namespace often{    inline void input(int &res)    {        char c = getchar();        res = 0;        int f = 1;        while (!isdigit(c))        {            f ^= c == '-';            c = getchar();        }        while (isdigit(c))        {            res = (res << 3) + (res << 1) + (c ^ 48);            c = getchar();        }        res = f ? res : -res;    }    inline int qpow(int a, int b)    {        int ans = 1, base = a;        while (b)        {            if (b & 1)                ans = (ans * base % mod + mod) % mod;            base = (base * base % mod + mod) % mod;            b >>= 1;        }        return ans;    }    int fact(int n)    {        int res = 1;        for (int i = 1; i <= n; i++)            res = res * 1ll * i % mod;        return res;    }    int C(int n, int k)    {        return fact(n) * 1ll * qpow(fact(k), mod - 2) % mod * 1ll * qpow(fact(n - k), mod - 2) % mod;    }    int exgcd(int a, int b, int &x, int &y)    {        if (b == 0)        {            x = 1, y = 0;            return a;        }        int res = exgcd(b, a % b, x, y);        int t = y;        y = x - (a / b) * y;        x = t;        return res;    }    int invmod(int a, int mod)    {        int x, y;        exgcd(a, mod, x, y);        x %= mod;        if (x < 0)            x += mod;        return x;    }}using namespace often;using namespace std;int n;signed main(){    auto polymul = [&](vector<double> &v, double er) {        v.emplace_back(0);        vector<double> _ = v;        int m = sz(v);        for (int i = 1; i < m; i++)            v[i] += er * _[i - 1];    };    auto polyval = [&](vector<double> const &c, double const &_x) -> double {        double res = 0.0;        int m = sz(c);        for (int ii = 0; ii < m; ii++)            res += c[ii] * pow(_x, (double)(m - ii - 1));        return res;    };    int __ = 1;    //input(_);    while (__--)    {        double _x, temp;        cin >> _x;        vector<double> x, y;        while (cin >> temp)            x.emplace_back(temp), y.emplace_back(cos(temp));        n = x.size();        vector<vector<double>> a(n, vector<double>(n));        int i, j;        for0(i, n) a[i][0] = y[i];        forab(j, 1, n - 1) forab(i, j, n - 1) a[i][j] = (a[i][j - 1] - a[i - 1][j - 1]) / (x[i] - x[i - j]);        vector<double> v;        v.emplace_back(a[n - 1][n - 1]);        forba(i, 0, n - 2)        {            polymul(v, -x[i]);            int l = sz(v);            v[l - 1] += a[i][i];        }        for0(i, n)            pr("%g/n", v[i]);        for0(i, n)        {            for0(j, n)                pr("%g ", a[i][j]);            puts("");        }        double _y =  polyval(v, _x);        pr("%g/n", _y);        pr("%.6e",fabs(_y-cos(_x)));    }    return 0;}

'''Author: cscDate: 2021-04-29 23:00:57LastEditTime: 2021-04-30 09:58:07LastEditors: Please set LastEditorsDescription: In User Settings EditFilePath: /code_py/newton_.py'''import numpy as npdef difference_quotient(x, y):    n = len(x)    a = np.zeros([n, n], dtype=float)    for i in range(n):        a[i][0] = y[i]    for j in range(1, n):        for i in range(j, n):            a[i][j] = (a[i][j-1]-a[i-1][j-1])/(x[i]-x[i-j])    return adef newton(x, y, _x):    a = difference_quotient(x, y)    n = len(x)    s = a[n-1][n-1]    j = n-2    while j >= 0:        s = np.polyadd(np.polymul(s, np.poly1d(            [x[j]], True)), np.poly1d([a[j][j]]))        j -= 1    for i in range(n):        print('%g' % s[n-1-i])    for i in range(n):        for j in range(n):            print('%g' % a[i][j], end=' ')        print()    _y = np.polyval(s, _x)    print('%g' % _y)    # re_err    real_y = np.cos(_x)    err = abs(_y-real_y)    print('%.6e' % err)def main():    _x = float(input())    x = list(map(float, input().split()))    y = np.cos(x)    newton(x, y, _x)if __name__ == '__main__':    main()