拉格朗日插值學習筆記
簡介
在數值分析中，拉格朗日插值法是以法國十八世紀數學家約瑟夫·拉格朗日命名的一種多項式插值方法。許多實際問題中都用函數來表示某種內在聯系或規律，而不少函數都只能通過實驗和觀測來了解。如對實踐中的某個物理量進行觀測，在若干個不同的地方得到相應的觀測值，拉格朗日插值法可以找到一個多項式，其恰好在各個觀測的點取到觀測到的值。這樣的多項式稱為拉格朗日（插值）多項式。數學上來說，拉格朗日插值法可以給出一個恰好穿過二維平面上若干個已知點的多項式函數。拉格朗日插值法最早被英國數學家愛德華·華林於1779年發現，不久后（1783年）由萊昂哈德·歐拉再次發現。1795年，拉格朗日在其著作《師范學校數學基礎教程》中發表了這個插值方法，從此他的名字就和這個方法聯系在一起。
插值
對於學習拉格朗日插值之前，我們首先要弄懂，什么叫插值？
插值 數學領域數值分析中的通過已知的離散數據求未知數據的過程或方法。
概念可能不好理解，舉個栗子：
我們給出：
\(x_0=1,y_0=3\)
\(x_1=2,y_0=5\)
\(x_2=4,y_0=8\)
\(x_3=5,y_0=4\)
求當 \(x=3\) 時，\(y\) 的值。
首先，我們可以通過列方程的辦法，求出函數，將 \(x\) 帶入求解。
如，我們設 \(y=ax^3+bx^2+cx+d\)
\(\begin{cases} a+b+c+d=3 \\8a+4b+2c+d=5 \\ 64a+16b+4c+d=8 \\125a+25b+5c+d=4 \end{cases}\)
然后聯立求解，這里就不給出求解過程了，需要用到高斯消元。
這樣做的復雜度是 \(O(n^3)\) 的，有點慢，我們考慮用拉格朗日差值優化這一問題。
拉格朗日插值
我們首先考慮 \(x_0\) ,建立一個函數，使其經過 \((x_0,1),(x_1,0),(x_2,0),(x_3,0)\)。
為什么這么構造，為什么不構造一個經過 \((x_0,2),(x_1,1),(x_2,1),(x_3,1)\) 的函數或者其他的呢，這個我們一會再說。
我們先考慮一下，如何才能構造這樣一個函數呢？
設
\(y=\frac{(X-x_1)\times(X-x_2)\times(X-x_3)}{(x_0-x_1)\times(x_0-x_2)\times(x_0-x_3)}\)
我們不難發現，當 \(X=x_0\) 時，分子變成了 \((x_0-x_1)\times(x_0-x_2)\times(x_0-x_3)\)，這與分母相同，於是式子可化為 \(1\)。再看，當 \(X=x_1,X=x_2,X=x_3\) 時，分子會變為 \(0\),於是，我們完成了構造的一大半。
我們發現，當前的函數只是經過了 \((x_0,1)\)，我們想要的是它經過 \((x_0,y_0)\)，怎么辦呢。
很簡單，我們直接在原始的基礎上乘上一個 \(y_0\)。
\(y=\frac{(X-x_1)\times(X-x_2)\times(X-x_3)}{(x_0-x_1)\times(x_0-x_2)\times(x_0-x_3)}\times y_0\)
此時，由於當 \(X=x_1,X=x_2,X=x_3\) 時，函數值為 \(0\)，我們乘上一個 \(y_0\) 不會對其有影響，於是，我們就構造出一個經過 \((x_0,y_0)\) 的函數。
如圖：
此時，我們再回答之前的問題：為什么要使函數經過 \((x_0,1),(x_1,0),(x_2,0),(x_3,0)\)？
按照我自己的理解，這其實是為了之后乘 \(y_0\) 做鋪墊，因為這樣我們的 \(x_0\) 就對應了 \(y_0\) ,而其他值依舊為 \(0\) ,滿足了我們的構造要求：只經過 \((x_0,y_0)\) 。而如果是構造一個經過 \((x_0,2),(x_1,1),(x_2,1),(x_3,1)\) 的函數，還是應該有辦法讓它經過 \(y_0\) 的，我沒有嘗試過，但是應該會非常麻煩，所以我覺得這個構造的方式也是拉格朗日插值的一個精妙之處。
返回題目，這樣做還是不夠的，我們的目的是為了構造經過這四個點的函數，這才經過了一個呢。
我們根據上面的方法，依次構造出只經過 \((x_1,y_1)\) 的函數，只經過 \((x_2,y_2)\) 的函數，只經過 \((x_3,y_3)\) 的函數，
\(y0=\frac{(X-x_1)\times(X-x_2)\times(X-x_3)}{(x_0-x_1)\times(x_0-x_2)\times(x_0-x_3)}\times y_0\)
\(y1=\frac{(X-x_0)\times(X-x_2)\times(X-x_3)}{(x_1-x_0)\times(x_1-x_2)\times(x_1-x_3)}\times y_1\)
\(y2=\frac{(X-x_0)\times(X-x_1)\times(X-x_3)}{(x_2-x_0)\times(x_2-x_1)\times(x_2-x_3)}\times y_2\)
\(y3=\frac{(X-x_0)\times(X-x_1)\times(X-x_2)}{(x_3-x_0)\times(x_3-x_1)\times(x_3-x_2)}\times y_3\)
然后，我們把它們加起來。
\(f(x)=y0+y1+y2+y3\)
圖像如下：
由於這四個函數都只經過了各不相同的四個點，其余點的值都為 \(0\)，因此，這個解析式就滿足了我們的要求，我們此時只要將待求的 \(x\) 帶入式子，就可以求出對應的值。
如我們的栗子：
我們再回頭看我們的函數
\(f(x)=y0+y1+y2+y3\)
\(=\frac{(X-x_1)\times(X-x_2)\times(X-x_3)}{(x_0-x_1)\times(x_0-x_2)\times(x_0-x_3)}\times y_0+\frac{(X-x_0)\times(X-x_2)\times(X-x_3)}{(x_1-x_0)\times(x_1-x_2)\times(x_1-x_3)}\times y_1+\frac{(X-x_0)\times(X-x_1)\times(X-x_3)}{(x_2-x_0)\times(x_2-x_1)\times(x_2-x_3)}\times y_2+\frac{(X-x_0)\times(X-x_1)\times(X-x_2)}{(x_3-x_0)\times(x_3-x_1)\times(x_3-x_2)}\times y_3\)
\(=y_0\times\prod\limits_{j=1}^{3}\frac{x-x_j}{x_0-x_j}+y_1\times \prod\limits_{j=0,j\not=1}^{3}\frac{x-x_j}{x_1-x_j}+y_2\times\prod\limits_{j=0,j\not=2}^{3}\frac{x-x_j}{x_2-x_j}+y_3\times\prod\limits_{j=0,j\not=3}^{3}\frac{x-x_j}{x_3-x_j}\)
\(=\sum\limits_{i=0}^{3}y_i\prod\limits_{i\not=j}\frac{x-x_j}{x_i-x_j}\)
我們把界限設定為 \(n\)，就可以得到常見的拉格朗日插值的式子：
\(f(x)=\sum\limits_{i=0}^{n}y_i\prod\limits_{j\not=i}\frac{x-x_j}{x_i-x_j}\)
此時，這道模板題我們就可以直接利用這個公式完成了。
附上代碼：
/*
work by smyslenny
2021.06.27
P4781 【模板】拉格朗日插值
*/
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <cstring>
#include <cstdlib>
#include <queue>
#include <vector>
#define int long long
#define INF 0x3f3f3f3f
using namespace std;
const int mod=998244353;
const int M=2e3+5;
int x[M],y[M];
int n,k,Ans;
int read()
{
register int x=0,y=1;
register char c=getchar();
while(c<'0'||c>'9') {if(c=='-') y=0;c=getchar();}
while(c>='0'&&c<='9') { x=x*10+(c^48),c=getchar();}
return y?x:-x;
}
int ksm(int a,int b)
{
int res=1;
while(b)
{
if(b&1) res=res*a%mod;
a=a*a%mod;
b>>=1;
}
return res;
}
int ny(int a)
{
return ksm(a,mod-2);
}
signed main()
{
n=read(),k=read();
for(int i=1;i<=n;i++) x[i]=read(),y[i]=read();
for(int i=1;i<=n;i++)
{
int fz=y[i]%mod,fm=1;
for(int j=1;j<=n;j++)
{
if(i==j) continue;
fz=fz*(k-x[j])%mod;
fm=fm*(x[i]-x[j])%mod;
}
Ans=((Ans+fz*ny(fm)%mod)+mod)%mod;
}
printf("%lld\n",(Ans+mod)%mod);
return 0;
}
上面的方式的復雜度是 \(O(n^2)\) 的，對於一些題目中，給出的 \(x\) 的取值是連續的，那么此時我們可以將復雜度優化到 \(O(n)\)。
\(x\)取值連續時的優化
對於我們的式子，我們可以知道左邊的求和 \(O(n)\) 是不能動的，我們考慮使右邊的求積優化成 \(O(1)\) 的。
由於 \(x\) 是連續的，我們可以直接將式子轉化成 \(f(x)=\sum\limits_{i=0}^ n y_i\prod\limits_{j\not=i}\frac{x-j}{i-j}\)。
對於分子，我們將其展開 \((x-1)\times(x-2)\times(x-3)\times\dots\times(x-(i-1))\times(x-(i+1))\times\dots\times(x-(j-1))\times(x-j)\)
我們維護一個關於 \(x\) 的前綴積和后綴積（類比前綴和），\(pre[i]=\prod\limits_{j=0}^i(x-j),suc[i]=\prod\limits_{j=i}^n(x-j)\)
分子就可以寫成 \(pre[i-1]\times suc[i+1]\)。
再看分母，拆開
\(i\times(i-1)\times(i-2)\times(i-3)\times\dots\times (i-(i-1))\times(i-(i+1))\times\dots\times(i-n)\)
觀察一下，\((i-(i-1))=1\) ，在這一塊之前的可以看做是 \(1\) ~ \(i\) 的乘積，也就是 \(i!\) ，后面 \((i-(i+1))=-1\) 后面我們可以看做是 \(-1\backsim(i-n)\) 的階乘，為了好計算，我們可以看作是 \(1\backsim(n-i)\) 的階乘。
PS： 這里我們要注意正負號，當所乘的總數為奇數個時，需要帶上負號，由於公式中不知道個數，下面默認為是正數，請不要認為就是正的。
所以分母可以化成 \(i!\times(n-i)!\)
右邊也就可以化成 \(\frac{pre[i-1]\times suc[i+1]}{i!\times(n-i)!}\) ,分子和分母都可以 \(O(n)\) 預處理，查詢是 \(O(1)\) 的，復雜度就可以優化成 \(O(n)\)。
一道模板題
附代碼：
/*
work by smyslenny
2021.06.27
CF622F The Sum of the k-th Powers
*/
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <cstring>
#include <cstdlib>
#include <queue>
#include <vector>
#define int long long
#define INF 0x3f3f3f3f
using namespace std;
const int M=1e6+5;
const int mod=1e9+7;
int n,k,tmp,Ans,pre[M],suc[M],fac[M];
int read()
{
register int x=0,y=1;
register char c=getchar();
while(c<'0'||c>'9') {if(c=='-') y=0;c=getchar();}
while(c>='0'&&c<='9') { x=x*10+(c^48),c=getchar();}
return y?x:-x;
}
void init()
{
pre[0]=suc[k+3]=fac[0]=1;
for(int i=1;i<=k+2;i++) pre[i]=pre[i-1]*(n-i)%mod;
for(int i=k+2;i>=1;i--) suc[i]=suc[i+1]*(n-i)%mod;
for(int i=1;i<=k+2;i++) fac[i]=fac[i-1]* i %mod;
}
int ksm(int a,int b)
{
int res=1;
while(b)
{
if(b&1) res=res*a%mod;
a=a*a%mod;
b>>=1;
}
return res;
}
signed main()
{
n=read(),k=read();
init();
for(int i=1;i<=k+2;i++)
{
tmp=(tmp+ksm(i,k))%mod;
int fz=pre[i-1]*suc[i+1]%mod,fm=fac[i-1]*((k-i)&1?-1:1)*fac[k+2-i]%mod;
Ans=(Ans+tmp*fz%mod*ksm(fm,mod-2)%mod)%mod;
}
printf("%lld\n",(Ans%mod+mod)%mod);
return 0;
}
重心拉格朗日插值法
不難發現，當我們每插入一個值的時候，都需要推倒重算，這里可以用重心拉格朗日插值法來解決。
\(f(x)=\sum\limits_{i=0}^{n}y_i\prod\limits_{j\not=i}\frac{x-x_j}{x_i-x_j}\)
\(=\frac{\prod\limits_{i\not=j}x-x_j}{\prod\limits_{i!=j}x_i-x_j}\)
上下乘 \(x-x_i\)
\(=\sum\limits_{i=1}^{n}y_i\frac{\prod\limits_{i=1}^n x-x_j}{\prod\limits_{i\not=j}(x_i-x_j)\times(x-x_i)}\)
\(=\prod\limits_{i=1}^n (x-x_j)\sum\limits_{i=1}^{n}\frac{y_i}{\prod\limits_{i\not=j}(x_i-x_j)\times(x-x_i)}\)
設 \(\lambda=\prod\limits_{i=1}^n (x-x_j)，\mu(x)=\prod\limits_{i\not=j}(x_i-x_j)\)
\(f(x)=\lambda\sum\limits_{i=1}^n\frac{y_i}{\mu(x)\times(x-x_i)}\)
插入一個值，只用 \(O(n)\) 更新 \(\mu(x)\) ，因為分子我們已經乘上 \((x-x_i)\) 了
詢問一個值，\(O(n)\) 更新 \(\lambda\) ，再套公式即可。
參考資料：