拉格朗日插值法學習筆記

拉格朗日插值法學習筆記

​ 在數值分析中，拉格朗日插值法是以法國十八世紀數學家約瑟夫·拉格朗日命名的一種多項式插值方法。許多實際問題中都用函數來表示某種內在聯繫或規律，而很多函數都只能經過實驗和觀測來了解。如對實踐中的某個物理量進行觀測，在若干個不一樣的地方獲得相應的觀測值，拉格朗日插值法能夠找到一個多項式，其剛好在各個觀測的點取到觀測到的值。這樣的多項式稱爲拉格朗日（插值）多項式。數學上來講，拉格朗日插值法能夠給出一個剛好穿過二維平面上若干個已知點的多項式函數。拉格朗日插值法最先被英國數學家愛德華·華林於1779年發現，不久後（1783年）由萊昂哈德·歐拉再次發現。1795年，拉格朗日在其著做《師範學校數學基礎教程》中發表了這個插值方法，今後他的名字就和這個方法聯繫在一塊兒。

​ 給定 \(n + 1​\)個橫座標不相同的點，能夠惟一肯定一個 \(n​\) 次的多項式，咱們能夠直接列方程利用高斯消元求解，時間複雜度\(O(n^3)​\)。

​ 而另外一種方法拉格朗日插值能夠在\(O(n^2)\)的時間複雜度內完美的求解出答案，能夠說是十分完美了。

拉格朗日插值法

爲了接下來的講解方便，咱們引進一些函數和變量來輔助講解。

拉格朗日插值法基函數（後文內簡稱爲基函數）\(g(x)​\)：

下面給出\(g(x)​\)的定義式：\(g(x)=\Pi_{i=1,i!=x}^{n}\frac{k-X[j]}{X[x]-X[j]}​\)

根據拉格朗日插值法，\(f(k)=\sum_{i=1}^{n}y_i*g(i)​\)

因而，咱們就能夠在\(O(n^2)​\)的複雜度內求出\(f(k)​\)。

咱們來看一個例子：

對於\(f(x)\),咱們有\(f(1)=1,f(2)=2,f(3)=3\)，求\(f(100)\)。

咱們先求出\(g(1,2,3)​\)：

\(g(1)=\frac{(100-2)*(100-3)}{(1-2)*(1-3)}=4753​\)，

\(g(2)=\frac{(100-1)*(100-3)}{(2-1)*(2-3)}=-9603​\)，

\(g(3)=\frac{(100-1)*(100-2)}{(3-1)*(3-2)}=4851​\)。

咱們有\(f(100)=g(1)*1+g(2)*2+g(3)*3=4753-9603*2+4851*3=100\)。

同時，咱們將\(f(k)​\)暴力展開有：

\(f(k)=1*\frac{(k-2)*(k-3)}{(1-2)*(1-3)}+2*\frac{(k-1)*(k-3)}{(2-1)*(2-3)}+3*\frac{(k-1)*(k-2)}{(3-1)*(3-2)}​\),

咱們發現若咱們將\(x_i\)帶入時，其餘項中都有\(x_i-x_i\)的項，一定爲\(0\)，因此這個正確性是很顯然的。

這是一道模板題：洛谷P4781

Ma code

#include <bits/stdc++.h>

using namespace std;

#define int long long
#define reg register
#define clr(a,b) memset(a,b,sizeof a)
#define Mod(x) (x>=mod)&&(x-=mod)
#define abs(a) ((a)<0?-(a):(a))
#define debug(x) cerr<<#x<<"="<<x<<endl;
#define debug2(x,y) cerr<<#x<<"="<<x<<" "<<#y<<"="<<y<<endl;
#define debug3(x,y,z) cerr<<#x<<"="<<x<<" "<<#y<<"="<<y<<" "<<#z<<"="<<z<<endl;
#define rep(a,b,c) for(reg int a=(b),a##_end_=(c); a<=a##_end_; ++a)
#define ret(a,b,c) for(reg int a=(b),a##_end_=(c); a<a##_end_; ++a)
#define drep(a,b,c) for(reg int a=(b),a##_end_=(c); a>=a##_end_; --a)
#define erep(i,G,x) for(int i=(G).Head[x]; i; i=(G).Nxt[i])
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize(3,"Ofast","inline")

inline int Read(void) {
    int res = 0, f = 1;
    char c;
    while (c = getchar(), c < 48 || c > 57)if (c == '-')f = 0;
    do res = (res << 3) + (res << 1) + (c ^ 48);
    while (c = getchar(), c >= 48 && c <= 57);
    return f ? res : -res;
}

template<class T>inline bool Min(T &a, T const&b) {
    return a > b ? a = b, 1 : 0;
}
template<class T>inline bool Max(T &a, T const&b) {
    return a < b ? a = b, 1 : 0;
}

const int N=2e3+5,M=1e5+5,K=10,mod=998244353,INF=2e9;

bool MOP1;

namespace LP {
    int n,k,A[N],B[N];
    inline int Pow(int x,int y) {
        int res=1;
        while(y) {
            if(y&1)res=res*x%mod;
            x=x*x%mod,y>>=1;
        }
        return res;
    }
    inline int Inv(int x) {
        return Pow(x,mod-2);
    }
    inline void solve(void) {
        n=Read(),k=Read();
        rep(i,1,n)A[i]=Read(),B[i]=Read();
        int Ans=0;
        rep(i,1,n) {
            int L1=1,L2=1;
            rep(j,1,n)if(i^j) {
                L1=L1*(k-A[j])%mod;
                L2=L2*(A[i]%mod-A[j]%mod)%mod;
            }
            Ans+=(L1*Inv(L2)%mod)*B[i]%mod,Mod(Ans);
            Ans+=mod,Mod(Ans);
        }
        printf("%lld\n",Ans);
    }
}

bool MOP2;

void _main(void) {
    LP::solve();
}

signed main() {
    _main();
    return 0;
}

在\(x_i\)連續時的作法

根據拉格朗日插值法，咱們有\(g(x)=\Pi_{i=1,i!=x}^{n}\frac{k-X[j]}{X[x]-X[j]}​\),\(f(k)=\sum_{i=1}^{n}y_i*g(i)​\)

關於\(f(k)\)的\(O(n)\)的複雜度已經沒法在進行優化了，複雜度主要堆積在\(g(x)\)的求值上，咱們須要優化\(g(x)\)的求值過程。

咱們先用\(i​\)來代替\(x_i​\)可得：\(g(x)=\Pi_{i=1,i!=x}^{n}\frac{k-j}{x-j}​\)，

仔細看一下，咱們發現了什麼呢？

先看一下分子，咱們能夠利用前綴和的思想，維護一個前綴積和後綴積。

即\(pre_i=\Pi_{j=0}^{i}k-i,nxt_i=\Pi_{j=i+1}^{n}k-j\)。

咱們發現分子就變成了\(pre_{i-1}*nxt_{i+1}​\)，能夠在\(O(1)​\)的複雜度內求出。

在來看一下分母，這不就是一個階乘的形式嗎？

咱們有分母\(=Frac[i]*Frac[n-i]*(-1)^{n-i}​\)。

因此咱們對於\(g(i)​\)有\(g(i)=\frac {pre_{i-1}*nxt_{i+1}}{Frac[i]*Frac[n-i]*(-1)^{n-i}}​\)

因而咱們就有了\(f(x)=\sum_{i=1}^{n}{\frac {pre_{i-1}*nxt_{i+1}}{Frac[i]*Frac[n-i]*(-1)^{n-i}}}\)。

因而咱們就能夠在\(O(n)\)的複雜度內解出了。。。

重心拉格朗日插值法

先咕咕咕了