淺談拉格朗日插值
在數值分析中,拉格朗日插值法是以法國十八世紀數學家約瑟夫·拉格朗日命名的一種多項式插值方法。許多實際問題中都用函數來表示某種內在聯系或規律,而不少函數都只能通過實驗和觀測來了解。拉格朗日插值法可以找到一個多項式,其恰好在各個觀測的點取到觀測到的值。這樣的多項式稱為拉格朗日(插值)多項式。
——百度百科
通俗地說,拉格朗日插值法可以找出一個恰好經過直角坐標系內\(n\)個給定點的函數
眾所周知,\(n\)個點能唯一確定的多項式最高次數是\(n-1\)次
這個可以兩點確定一次函數,三點確定二次函數來推出,或者我們由方程組有唯一解的充要條件也能得出
知道了這個之后我們就容易想到直接用高斯消元來搞出多項式的系數,但是復雜度消耗太大,是\(O(n^3)\)
而拉格朗日插值法就是一種一般可以在\(O(n^2)\)的復雜度下求出多項式的方法(不過通常用來求一個點值,以下的講述一般以此為參考)
一般的拉格朗日插值法
我們設要求的多項式為\(f(x)\),點的坐標為\((x_i,y_i)\),我們要找多項式在\(x_0\)處的取值
先上公式,由拉格朗日插值法:
看起來不是很好理解,其實很簡單,我們把原來的某個給定點\(x_k\)代入以下有:
容易發現,當\(k\not=i\)時,后面的\(\frac{x_k-x_j}{x_i-x_j}\)的分子總有一項是\(0\),此時\(\prod_{i\not= j} \frac{x_k-x_j}{x_i-x_j}=0\)
當\(k=i\)時,后面的\(\frac{x_k-x_j}{x_i-x_j}\)上下完全相同,此時\(\prod_{i\not= j} \frac{x_k-x_j}{x_i-x_j}=1\)
即對於\(f(x_k)\)來說,這個多項式的確給出了對應的\(y_k\)的值
不難發現這個方法對所有點都適用,因此它是正確的
從上面的式子可以看出每次計算要枚舉兩次,因此復雜度很簡單,就是\(O(n^2)\)
在\(x\)取值連續時的插值法
因為很多時候我們做題都是先發現某個函數是多少次的多項式,然后自己隨意取一些值代入插值
這樣的話為了省事橫坐標的取值完全可以從\(1\)開始連續取,那么我們把上面的式子中的\(x_i\)換成\(i\)就有:
考慮怎么快速求\(\prod_{i\not= j} \frac{k-j}{i-j}\),我們分別考慮:
分子的話容易發現就是\(\frac{\prod_{t=1}^n x_0-t}{x_0-i}\)
分母比較復雜,\(i-j\)的累乘可以分成兩個階乘部分,因此推導一下就是\((-1)^{n-i}\cdot i!\cdot(n-i)!\)
這樣我們一般就可以\(O(n\log n)\)來算了(\(\log\)主要是有求逆元的過程,當然你預處理\(O(n)\)求也沒有問題)
重心拉格朗日插值法
我們考慮到朴素的拉格朗日插值每次多加入一個點時就要整個重新算過,很浪費時間
那么能不能把重復算的一些東西利用起來?
我們對於
把分子提取出來,設為\(g=\prod_{i=1}^n x_0-x_i\),則此時:
設一個\(t_i=\frac{y_i}{\prod_{i\not= j} x_i-x_j}\),則:
因此每次多加入一個點只需要重新\(O(n)\)算它的\(t_i\)就好了
拉格朗日插值的應用以及常用解題思路
一個經典例子:自然數冪和,即求\(\sum_{i=1}^n i^k\)
之前也提到過用第二類斯特林數做的方法,但那種方法是\(O(k^2)\)的,不夠優秀
但是現在我們觀察這個式子,如果不看求和的話\(i^k\)就是一個\(k\)次多項式
那么前綴和(其實就是差分)之后,次數要\(+1\),即此時答案為一個關於\(n\)的\(k+1\)多項式
那我們直接代入\(k\)個值(取值連續,反正自己定)之后插值算即可,復雜度\(O(k\log k)\)
那么很多具體題目怎么辦呢,一般就是先推出某個式子,然后證明它是關於某個自變量的多少次函數(一定要判斷出次數!),然后自己選一些點(一般是連續的)代入插值即可
其實做了一些題目之后就很套路了
拉格朗日插值的模板
模板看Luogu P4781 【模板】拉格朗日插值 ,就是用一般的\(O(n^2)\)來做就好了
CODE
#include<cstdio>
#define RI register int
#define CI const int&
using namespace std;
const int N=2005,mod=998244353;
int n,x[N],y[N],k;
inline int sub(CI x,CI y)
{
int t=x-y; return t<0?t+mod:t;
}
inline int inv(int x,int p=mod-2,int mul=1)
{
for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
inline int Lagrange(CI n,int *x,int *y,CI k)
{
int ret=0; for (RI i=1;i<=n;++i)
{
int s1=1,s2=1; for (RI j=1;j<=n;++j) if (i!=j)
s1=1LL*s1*sub(k,x[j])%mod,s2=1LL*s2*sub(x[i],x[j])%mod;
(ret+=1LL*y[i]*s1%mod*inv(s2)%mod)%=mod;
}
return ret;
}
int main()
{
scanf("%d%d",&n,&k); for (RI i=1;i<=n;++i) scanf("%d%d",&x[i],&y[i]);
return printf("%d",Lagrange(n,x,y,k)),0;
}
一些入門例題
- CF622F The Sum of the k-th Powers 提到過的自然數冪和
- Luogu P3270 [JLOI2016]成績比較 需要容斥+組合數學+廣義容斥化式子的自然數冪和
- Luogu P4463 [國家集訓隊] calc DP+多項式差分系數推導之后用拉格朗日插值算答案