拉格朗日插值學習筆記

逛知乎時偶然看到了一個很經典的找規律填數問題，然后下面的回答基本都是 114514惡臭，突然想知道大伙是如何構造出這種能填入惡臭數字的函數的，於是就去了解了一波插值，於是就學了一波拉格朗日插值，於是就有了這篇博客。

引入

眾所周知，\(n+1\)​​ 個點 \((x_i,y_i)\)​​ 可以唯一地確定一個 \(n\)​​ 次多項式，給你 \(n+1\)​​ 個點，讓你求出所對應的那個多項式，怎么做？

常規的套路是待定系數法列出方程組，然后用高斯消元 \(O(n^3)\)​​​​​ 求解，顯然時間復雜度過大，還容易出精度問題，還有什么更優的方法嗎？有，拉格朗日插值，時間復雜度 \(O(n^2)\)​​​​​，且代碼實現起來基本不會有精度問題。

說到拉格朗日插值，就不得不先提一下到底什么是插值了。簡單來說，插值就是給你幾個離散的點 \(\big(x,f(x)\big)\)​​​​​​​​​​​​​，讓你根據這些點求出一個函數，滿足這些點都在這個函數上，稱為插值函數；然后一般會給你幾個新的 \(x\)​​​​​​，讓你求出這些 \(x\)​​​​ 所對應的 \(f(x)\)​​​​​​​​​​​​​ 是多少，稱為內插，一般多用於數值分析領域。

顯然，插值的方法有很多種，根據不同的插值法，能求到不同的插值函數，一般把插值法按照它們所求得的插值函數類型分類。比如你用某個插值法求出的插值函數是一個多項式，那么就稱該插值法為多項式插值法，拉格朗日插值就是一個經典的多項式插值法，它可以在給定 \(n+1\) 個點的條件下求出一個 \(n\) 次多項式函數，滿足這些點都在這個函數上。

舉個栗子，給出 \((1,1),~(2,3),~(3,5),~(4,7)\)，問你 \(f(5)\) 是多少。

答案是 \(114514\)​​。

因為當這個函數是 \(\displaystyle f(x)=\frac {18111}2x^4-90555~x^3+\frac{633885}2x^2-452773~x+217331\) 時，有：

\[f(1)=1\\f(2)=3\\f(3)=5\\f(4)=7\\f(5)=114514 \]

這插值怎么這么臭啊（惱）

拉格朗日插值

目的是構造函數 \(F(x)\)​​​​​​，使得對於給定的點 \(x_i\)​​​​​​，都有 \(F(x_i)=y_i\)​​​​​​。怎么構造呢？我們可以讓 \(F(x)=y_1f_1(x)+y_2f_2(x)+\dots+y_{n+1}f_{n+1}(x)\)​​，​且 \(f_i(x)\) 只有當 \(x=x_i\) 時取值為 1，否則取值為 0，顯然這樣構造出來的 \(F(x)\) 是滿足條件的。

那如何構造 \(f_i(x)\)​​ 呢？簡單，只需讓 \(f_i(x)=\prod\limits_{j\neq i}\dfrac{x-x_j}{x_i-x_j}\)​​​ 即可。就拿上面的很臭的例子來說吧，我們有：

\[f_1(x)=\frac{(x-2)(x-3)(x-4)}{(1-2)(1-3)(1-4)}\\ f_2(x)=\frac{(x-1)(x-3)(x-4)}{(2-1)(2-3)(2-4)}\\ f_3(x)=\frac{(x-1)(x-1)(x-4)}{(3-1)(3-2)(3-4)}\\ f_4(x)=\frac{(x-1)(x-2)(x-3)}{(4-1)(4-2)(4-3)}\\ \]

顯然當 \(x\neq x_i\)​​ 時，\(f_i(x)\)​​ 的分子肯定有一項為 0，值為 0；當 \(x=x_i\) 時，\(f_i(x)\) 的分子分母相等，值為 1，即 \(f_i(x)\) 剛好滿足上文所說的條件。

此時 \(F(x)=f_1(x)+3f_2(x)+5f_3(x)+7f_4(x)\)，試着把幾個點值都帶進去，發現都對的上114514無視即可，於是滿足條件的插值函數 \(F(x)\)​​​​ 就成功構造出來了~~

應用&&例題

在算法競賽中，往往會遇到這樣一類問題，它們會直接或間接地給你幾個點，然后讓你求另一個點的值，這類問題往往會跟多項式問題掛上鈎。現在我們知道了拉格朗日插值的公式，就可以利用它解決一些經典的插值問題了，這里上幾個例題幫助說明。

【模板】拉格朗日插值 - 洛谷

題目大意

給出 \(n\)​​​ 個點和一個數 \(k\)​​​，讓你求出根據這 \(n\)​​​ 個點確定的多項式 \(y=F(x)\)​​​ 在 \(k\)​​​ 上的值，即 \(F(k)\)​​​，答案對 998244353 取模。

\((1\leq2e3,~1\leq x_i,y_i,k<998244353),~x_i~\)兩兩不同。

題目分析

模板題，根據拉格朗日插值，我們有：

\[F(x)=\sum_{i=1}^n\left(y_i\prod_{j=1,j\neq i}^n\frac{x-x_j}{x_i-x_j}\right) \]

直接把 \(x=k\)​​​​ 帶進去暴力算就行，時間復雜度 \(O(n^2)\)​​​​，記得分母先做連乘最后再求逆元，不然復雜度多個 \(\log\)​​​，只能得 60 分。我就是每次都對分母求逆元然后過不了

代碼實現

#include<bits/stdc++.h>
using namespace std;
const long long mod=998244353;
const int maxn=2e3+10;
long long x[maxn],y[maxn];
int n,k;
void exgcd(long long a,long long b,long long &x,long long &y){
	if(b==0){
		x=1,y=0;
		return ;
	}
	exgcd(b,a%b,y,x);
	y-=a/b*x;
}
long long inv(long long a){
	long long x,y;
	exgcd(a,mod,x,y);
	return x;
}
int main(){
	cin >> n >> k;
	for(int i=1;i<=n;i++){
		cin >> x[i] >> y[i];
	}
	long long ans=0;
	for(int i=1;i<=n;i++){
		long long fz=y[i];
		long long fm=1;
		for(int j=1;j<=n;j++){
			if(i==j) continue;
			long long t=k-x[j];
			fm*=x[i]-x[j];
			fm%=mod;
			fz=fz*t%mod;
		}
		fm<0?fm+=mod:0;
		fm=inv(fm);
		ans+=fz*fm%mod;
		ans%=mod;
	}
	ans<0?ans+=mod:0;
	cout << ans;
} 

---

CF622F The Sum of the k-th Powers - 洛谷

題目大意

給定 \(n\) 和 \(k\)，讓你求 \(\sum_{i=1}^ni^k\mod {(1e9+7)}\)。

\((1\leq n\leq 1e9,~0\leq k\leq1e6)\)​​

題目分析

乍一看好像跟拉格朗日插值沒啥關系，但眾所周知，次方數的前綴和是有公式的，而 \(k\)​​​​ 次方數的前綴和公式一定是個 \(k+1\)​​​​​ 次多項式。簡單的證明的話請看這里，當然方法不止一種，有興趣的同學可以自行去了解一下。

假設 \(k\)​​​​​​​​​ 次方數的前綴和公式 \(F_k(n)=\sum_{i=1}^ni^k\)​​​​​​​​​，則我們可以計算出 \(k+2\)​​​​​​​ 個點 \(\big(x,F_k(x)\big)\)​​​​​​​，然后利用拉格朗日插值求得 \(F_k(n)\)​​​​。

計算 \(k+2\)​​​​​​​​ 個點可以直接用快速冪 \(O(n\log n)\)​​​​​​​​ 求，也可以線性篩 \(x^k\)​​​​​​​​ 然后做前綴和 \(O(n)\)​​​​​ 求，問題不大，主要問題是拉格朗日插值的時間復雜度是 \(O(n^2)\)​​​​​ 的，\(\bf n\)​​​​​ 在這里是點的數量，即 \(k+2\)​​​​​，直接套板子必 T。

注意到橫坐標 \(x\)​ 的值是連續的，所以可以簡化運算至 \(O(n)\) 時間復雜度，推導如下：

\[F(n)=\sum_{i=1}^{k+2}\Bigg({\large y_i*f_i(n)}\Bigg)=\sum_{i=1}^{k+2}\left({\large y_i}\prod_{j=1,j\neq i}^{k+2}\frac{n-x_j}{x_i-x_j}\right)=\sum_{i=1}^{k+2}\left({\large y_i}\prod_{j=1,j\neq i}^{k+2}\frac{n-j}{i-j}\right) \]

我們可以 \(O(1)\)​ 求 \(f_i(n)=\displaystyle\prod\limits_{j=1,j\neq i}^{k+2}\frac{n-j}{i-j}\)​​ ：

\[分子=\prod_{j=1,j\neq i}^{k+2}(n-j)=\prod_{j=1,j\neq i}^{i-1}(n-j)*\prod_{j=i+1,j\neq i}^{k+2}(n-j)\\ 分母=\prod_{j=1,j\neq i}^{k+2}(i-j)=\prod_{j=1,j\neq i}^{i-1}(i-j)*\prod_{j=i+1,j\neq i}^{k+2}(i-j)=(-1)^{k+2-i}*(i-1)!*(k+2-i)! \]

顯然可以 \(O(n)\)​​​ 預處理分子和分母，然后 \(O(1)\)​​​ 求解 \(y_i*f_i(n)\)​​​，從而將總時間復雜度降為了 \(O(n)\)​​。

代碼實現

#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
const int maxn=1e6+6;
long long n,k;
long long prime[maxn];
bool isp[maxn];
long long y[maxn],f[maxn],invf[maxn],pre[maxn],suf[maxn];
long long qsm(long long a,long long b){
	long long ret=1;
	while(b){
		if(b&1){
			ret*=a,ret%=mod;
		}
		b>>=1,a*=a,a%=mod;
	}
	return ret;
}
void inip(){
	f[0]=invf[0]=f[1]=invf[1]=1;
	y[1]=1;
	for(int i=2;i<maxn;i++){
		f[i]=f[i-1]*i%mod;
		if(!isp[i]){
			y[i]=qsm(i,k);
			prime[++prime[0]]=i;
		}
		for(int j=1;j<=prime[0]&&i*prime[j]<maxn;j++){
			isp[i*prime[j]]=true;
			y[i*prime[j]]=y[i]*y[prime[j]]%mod;
			if(i%prime[j]==0) break;
		}
	}
}
int main(){
	cin >> n >> k;
	inip();
	pre[0]=1,suf[k+3]=1,invf[k+3]=qsm(f[k+3],mod-2);
	for(int i=1;i<=k+2;i++){
		pre[i]=pre[i-1]*(n-i)%mod;
		y[i]+=y[i-1];
		y[i]%=mod;
	}
	for(int i=k+2;i>=1;i--){
		invf[i]=invf[i+1]*(i+1)%mod;
		suf[i]=suf[i+1]*(n-i)%mod;
	}
	if(n<=k+2){
		cout << y[n];
	}
	else {
		long long ans=0;
		for(int i=1;i<=k+2;i++){
			long long fz=pre[i-1]*suf[i+1]%mod;
			long long fm=invf[i-1]*((k+2-i&1)?mod-invf[k+2-i]:invf[k+2-i])%mod;
			ans+=y[i]*fz%mod*fm%mod;
			ans%=mod;
		}
		ans<0?ans+=mod:ans;
        cout << ans;
	}
} 

參考資料

• 插值（離散數學名詞）_百度百科 (baidu.com)
• 拉格朗日插值法(圖文詳解) - Angel_Kitty - 博客園 (cnblogs.com)
• 從拉插到快速插值求值 - command_block 的博客 - 洛谷博客 (luogu.com.cn)
• 題解 P4781 拉格朗日插值 - attack 的博客 - 洛谷博客 (luogu.com.cn)