拉格朗日插值法学习笔记

STEP1 问题引入

现在给定在平面直角坐标系中的 $k+1$ 个点，记 $P_i(x_i,y_i)$。现在要求一个 $k$ 次多项式函数 $y={\sum }_{i=0}^k{a}_i{x}^i$，使得它经过这 $k+1$ 个点，即 ${\forall }_{i|0\le i\le k},y{|}_{x=x_i}=y_i$。输出 $x=n$ 时 $y$ 的取值。

我们首先想到的是把 $k+1$ 个点的坐标分别代入 $y={\sum }_{i=0}^k{a}_i{x}^i$，则可以得到 $k+1$ 个方程。于是便可以用高斯消元法求出所有的 $a_i$，最后把 $n$ 代进去算就好了时间复杂度为 $O(k^3)$。这种方法效率非常低，需要用更高效的方法代替。

STEP2 拉格朗日插值公式

存在性

先考察多项式 $l_i={\prod }_{j\ne i}\frac{x-x_j}{x_i-x_j}$。不难发现，如果 $x=x_i$，那么每个分数分子与分母相等，即每个分数都为 $1$，则此时 $l_i=1$；如果 $x\ne x_i$且 $x$ 为某个点 $P_j$ 的横坐标，那么分子有一项为 $x-x_j=0$，那么此时 $l_i=0$。像这样的多项式，我们称之为拉格朗日基本多项式。
那么如果将 $l_i$ 乘上 $y_i$ ，就可以得到在 $x_i$ 上取值为 $y_i$，在 $x_j$ （$j\ne i$）上取值为 $0$ 的一系列多项式。我们把所有这种多项式累加起来，就可以得到要求的多项式$$L(x)=\sum _{i=0}^k{y}_i{l}_i=\sum _{i=0}^k{y}_i\sum _{0\le j\le k,j\ne i}\frac{x-x_j}{x_i-x_j}$$
时间复杂度是 $O(k^2)$。
由于每个点的横坐标不相等，$i\ne j$，故 $x_i\ne x_j$，即 $x_i-x_j\ne 0$。也就是说，上式中的所有分母不等于零，故这样的多项式是存在的。

唯一性

假设有两个不同的 $k$ 次多项式 $L_1\ne L_2$，记 $L=L_1-L_2$，则 $L$ 的次数不大于 $k$。又 $L_1$ 和 $L_2$ 在 $x_0,x_1,\dots ,x_k$ 上的取值相同，则 $L$ 在这些点上取值都为零，也就是说 $L=a{\prod }_{i=0}^k(x-x_i)$。若 $a\ne 0$，则这是一个 $k+1$次多项式，与它次数不大于 $k$ 相矛盾，故 $a=0$，那么 $L=0$ 即 $L_1=L_2$。这又与 $L_1\ne L_2$ 相矛盾。故有且只有一个 $k$ 次多项式经过平面上的 $k+1$ 个点。

代码模板

#include<iostream>
#include<cstdio>
#define ll long long
using namespace std;
int n;
ll x[2001],y[2001],mod=1e9+7;
ll k,ans,w,l;
inline ll pw(ll a,int b){  //快速幂
	if(b==0) return 1ll;
	if(b==1) return a;
	ll res=pw(a,b>>1);
	(res*=res)%=mod;
	if(b&1) (res*=a)%=mod;
	return res;
}
inline ll sub(ll a,ll b){  //模下的减法
	return (a-b>0)?(a-b):(a-b+mod);
}
inline ll inv(ll a){  //乘法逆元
	return pw(a,998244351);
}
inline void lagrange(){
	for(register int i=0;i<=k;i+=1){
		l=1ll; w=1ll;
		for(register int j=0;j<=k;j+=1){
			if(i==j) continue;
			(l*=sub(n,x[j]))%=mod;  //计算分子
			(w*=sub(x[i],x[j]))%=mod;  //计算分母
		}
		(l*=inv(w))%=mod;  //拉格朗日基本多项式的值
		(ans+=l*y[i])%=mod;  //累加答案
	}
	return;
}
int main(){
	scanf("%d%lld",&n,&k);
	for(register int i=0;i<=k;i+=1){  //下标从0开始
		scanf("%lld%lld",&x[i],&y[i]);
	}
	lagrange();
	printf("%lld\n",ans);
	return 0;
}

STEP3 算法优化

重心拉格朗日插值

定义 $l(x)={\prod }_{i=0}^k(x-x_i)$，那么有 $l_i(x)=\frac{l}{x-x_i}$。
原式变成 $L(x)={\sum }_{i=0}^k\frac{y_{i}l(x)}{(x-x_i){\prod }_{0\le j\le k,j\ne i}(x_i-x_j)}$。定义重心权 $w_i=\frac{1}{{\prod }_{0\le j\le k,j\ne i}(x_i-x_j)}$，则$$L(x)=l(x)\sum _{i=0}^k\frac{y_i{w}_i}{x-x_i}$$
这样 $l_i(x)$ 就能以 $O(1)$ 的时间计算出来，总共耗时 $O(n)$。

$x$连续取值的情况

当所有 $x$ 的取值为连续整数时，基本多项式的分母被 $x_i$ 分成两段。因此我们与处理好阶乘的乘法逆元，就可以以 $O(n)$ 的复杂度计算分母了。
上述两种方法结合起来，就能把算法的时间复杂度从 $O(n^2)$ 变为 $O(n)$。

STEP4 经典例题

The Sum of the k-th Powers（CF622F）

There are well-known formulas: ${\sum }_{i=1}^{n}i=1+2+\cdots +n=\frac{n(n+1)}{2}$,${\sum }_{i=1}^n{i}^2=1^2+2^2+\cdots +n^2=\frac{n(2n+1)(n+1)}{6}$ ,${\sum }_{i=1}^n{i}^3=1^3+2^3+{\dots }_n^3=(\frac{n(n+1)}{2}{)}^2$ . Also mathematicians found similar formulas for higher degrees.
Find the value of the sum ${\sum }_{i=1}^n{i}^k=1^k+2^k+\cdots +n^k$ modulo $10^9+7$ (so you should find the remainder after dividing the answer by the value $10^9+7$).
题目大意： 求 ${\sum }_{i=1}^n{i}^k(mod10^9+7)$，其中 $1\le n\le 10^9$，$0\le k\le 10^6$。

记 $F(x)={\sum }_{i=1}^x{i}^k$ 跟据题目的提示，这是一个 $k+1$ 次多项式函数。因此只要计算 $F(1),F(2),\dots ,F(k+2)$，再用拉格朗日插值法就能求得 $F(n)$。注意到选取的 $x$ 为连续整数，可以用到上面提到的两个优化，总时间复杂度为 $O(k)$。

#include<iostream>
#include<cstdio>
#define ll long long
using namespace std;
int n,k;
ll y[1000005],jc[1000005],l=1ll,w,ans,mod=1e9+7;
ll qp(ll a,int b){
	if(b==0) return 1ll;
	if(b==1) return a;
	ll res=qp(a,b>>1);
	(res*=res)%=mod;
	if(b&1) (res*=a)%=mod;
	return res;
}
ll sub(ll a,ll b){
	return (a>b)?(a-b):(a-b+mod);
}
ll inv(ll a){
	return qp(a,mod-2);
}
void lagrange(){
	for(int i=1;i<=k+2;i+=1){
		(w*=sub(n,i))%=mod;
	}
	jc[0]=1ll;
	for(int i=1;i<=k+2;i+=1){
		jc[i]=(jc[i-1]*i*1ll)%mod;  //阶乘预处理
	}
	for(int i=1;i<=k+2;i+=1){  //拉格朗日插值
		l=(jc[i-1]*jc[k+2-i])%mod;
		if((k+2-i)&1) l=sub(0,w);
		l=inv(l);
		(l*=inv(n-i))%=mod;
		(l*=w)%=mod;
		(l*=y[i])%=mod;
		(ans+=l)%=mod;
	}
	return;
}
int main(){
	scanf("%d%d",&n,&k);
	for(int i=1;i<=k+2;i+=1){
		y[i]=qp(i*1ll,k);
		(y[i]+=y[i-1])%=mod;  //计算y
	}
	if(n<=k+2){
		ans=y[n];  //这种情况直接可以获得答案
	}
	else lagrange();
	printf("%lld\n",ans);
	return 0;
}

谢谢观看，记得点赞