[筆記] 拉格朗日插值法

拉格朗日插值法

問題：給你 \(n+1\) 個點值，求這 \(n+1\) 個點確定的 \(n\) 次多項式 \(f(x)\)（求出給定點 \(x_0\) 的值 \(f(x_0)\) 即可）。

我們可以直接高斯消元，\(\mathcal{O}(n^3)\)

一般的拉格朗日插值法

簡單來說，拉格朗日插值法可以找出一個恰好經過直角坐標系內 \(n\) 個給定點的函數。

我們設所求的多項式為 \(f(x)\) ，點的坐標為 \((x_i,y_i)\) ，那么我們有：

\[f(x_0)=\sum_{i=1}^n y_i \prod_{i\neq j} \frac{x_0-x_i}{x_i-x_j} \]

我們發現，代入任意一個 \(x_k,k\in[1,n]\) ，當 \(k\neq i\) 時，后面的 \(\prod_{j\neq i} \frac{x_0-x_i}{x_i-x_j}\) 為 \(0\) ；而 \(k==i\) 時，后面的 \(\prod_{j\neq i} \frac{x_0-x_i}{x_i-x_j}\) 為 \(1\) 。

時間復雜度 \(\mathcal{O}(n^2)\)

在 \(x\) 取值連續時的拉格朗日插值法

很多時候我們做題都是先發現某個函數是多少次的多項式，然后自己取一些值代入插值求 \(f(x_0)\)。

這時我們爆算的點的橫坐標可以從 \(1\) 開始連續取，這樣我們把上面的式子中的 \(x_i\) 換成 \(i\) 有：

\[f(x_0)=\sum_{i=1}^{n}y_i \prod_{i\neq j} \frac{x_0−j}{i−j} \]

我們發現：

對於分子我們可以與處理前綴積和后綴積：

\(pl[i]=\prod_{j=1}^i (x_0-j)\)

\(pr[i]=\prod_{j=i}^n(x_0-j)\)

這樣可以 \(\mathcal{O}(1)\) 求分子；

對於分母，我們可以預處理階乘，分母實際上就是：

\((-1)^{n-i}(i-1)!\cdot (n-i)!\)

這樣我們可以一個 \(\mathcal{O}(\log n)\) 求出分母。（如果你預處理逆元則可以 \(\mathcal{O}(1)\) ）

重心拉格朗日插值法

我們發現，一般的拉格朗日插值法每多加入一個點時就要整個重新計算，考慮如何利用已經計算的信息。

我們對於

\(f(x_0)=\sum_{i=1}^n y_i \prod_{i\neq j} \frac{x_0−x_j}{x_i−x_j}\)

把分子提取出來，設為 \(g=\prod_{i=1}^n x_0-x_i\) ，則此時：

\(f(x_0)=g \times \sum_{i=1}^n \prod_{i\neq j} \frac{y_i}{(x_0-x_i)(x_i-x_j)}\)

設 \(t_i=\frac{y_i}{\prod_{i\neq j} x_i-x_j}\) 則：

\(f(x_0)=g\times \sum_{i=1}^n \frac{t_i}{x_0-x_i}\)

因此每次多加入一個點只需要重新 \(\mathcal{O}(n)\)算它的 \(t_i\) 就好了。

拉格朗日插值法的應用

暴力算式子：確定題目中的某個函數的次數，就可以暴力代入，然后差值去算出這個式子。

或是說題目給定的式子不確定，你需要在程序中算出來（有回校內考試出現過

拉格朗日插值法的板子（？

Luogu P4781 【模板】拉格朗日插值

#include<cstdio>
#include<iostream>
#define ll long long
#define R register int
using namespace std;
namespace Luitaryi {
template<class I> inline I g(I& x) { x=0; register I f=1;
  register char ch; while(!isdigit(ch=getchar())) f=ch=='-'?-1:f;
  do x=x*10+(ch^48); while(isdigit(ch=getchar())); return x*=f;
} const int N=2010,M=998244353;
int n,k;
int x[N],y[N];
inline ll qpow(ll a,ll b) { register ll ret=1;
  for(;b;b>>=1,(a*=a)%=M) if(b&1) (ret*=a)%=M; return ret;
}
inline int Lagrange(const int& n,const int& k,int* x,int* y) { R ret=0;
  for(R i=1;i<=n;++i) { R up=1,dn=1;
    for(R j=1;j<=n;++j) if(i!=j) 
      up=1ll*up*(k-x[j]+M)%M,dn=1ll*dn*(x[i]-x[j]+M)%M;
    ret=(ret+1ll*y[i]*up%M*qpow(dn,M-2))%M;
  } return ret;
}
inline void main() {
  g(n),g(k); for(R i=1;i<=n;++i) g(x[i]),g(y[i]);
  printf("%d\n",Lagrange(n,k,x,y));
}
} signed main() {Luitaryi::main(); return 0;}

CF622F The Sum of the k-th Powers

#include<bits/stdc++.h>
#define R register int
using namespace std;
namespace Luitaryi {
inline int g() { R x=0,f=1;
  register char s; while(!isdigit(s=getchar())) f=s=='-'?-1:f;
  do x=x*10+(s^48); while(isdigit(s=getchar())); return x*f;
} const int N=1000010,M=1000000007;
int n,k,ans;
int fac[N],pl[N],pr[N];
inline int qpow(int a,int b) { R ret=1;
  for(;b;b>>=1,a=1ll*a*a%M) if(b&1) ret=1ll*ret*a%M; return ret;
}
inline void main() { n=g(),k=g();
  fac[0]=pl[0]=pr[k+3]=1;
  for(R i=1;i<=k+2;++i) fac[i]=1ll*fac[i-1]*i%M;
  for(R i=1;i<=k+2;++i) pl[i]=1ll*pl[i-1]*(n-i)%M;
  for(R i=k+2;i;--i) pr[i]=1ll*pr[i+1]*(n-i)%M;
  for(R i=1,y=0,up,dn;i<=k+2;++i) {
    y=(y+qpow(i,k))%M;
    up=1ll*pl[i-1]*pr[i+1]%M;
    dn=((k-i)&1?-1ll:1ll)*fac[i-1]*fac[k+2-i]%M;
    ans=(ans+1ll*y*up%M*qpow(dn,M-2))%M;
  } printf("%d\n",(ans+M)%M);
}
} signed main() {Luitaryi::main(); return 0;}