拉格朗日插值学习笔记

众所周知，若已知一 $n$ 次多项式 $f(x)$ 在 $x=\{x_0,x_1,x_2,\dots,x_n\}$ 这 $n+1$ 个点的值 $\{y_0,y_1,y_2,\dots,y_n\}$，那么就可以唯一确定这个多项式。而拉格朗日插值就是一种根据 $\{x_0,x_1,x_2,\dots,x_n\}$ 和 $\{y_0,y_1,y_2,\dots,y_n\}$ 求出 $f(x)$ 的算法。

不难发现，若我们构造 $n+1$ 个函数 $f_{[0,n]}(x)$，满足 $g_i(x)=\begin{cases}1&x=x_i\\0&x\not=x_i\\??&otherwise\end{cases}$，那么有 $f(x)=\sum\limits_{i=0}^n y_ig_i(x)$。

考虑怎么构造 $g_i(x)$，当 $x\not=x_i$ 时 $g_i(x)=0$ 显然可以通过让 $g_i(x)=\prod\limits_{j\not=i}x-x_j$，但是当 $x=x_i$ 时 $g_i(x)=\prod\limits_{j\not=i}x_i-x_j$，所以可以让这个东西除掉它，令 $g_i(x)=\frac{\prod\limits_{j\not=i}x-x_j}{\prod\limits_{j\not=i}x_i-x_j}$ 即可，这样就可以满足以上两个条件。

综上，有：

$\begin{aligned} f(x)&=\sum\limits_{i=0}^n y_ig_i(x)\\ &=\sum\limits_{i=0}^n y_i\frac{\prod\limits_{j\not=i}x-x_j}{\prod\limits_{j\not=i}x_i-x_j}\\ \end{aligned}$

按照这个公式直接插值的时间复杂度是 $O(n^2)$ 的。

看一道例题：

给 $n+1$ 个点，请你求出经过这 $n+1$ 个点的 $n$ 次多项式 $f(x)$ 在 $x=k$ 时的取值，对 $998244353$ 取模。

$1\le n\le 2\times 10^3$。

直接 $O(n^2)$ 插值然后计算即可。

代码如下：

#include <iostream>
#include <cstdio>

using namespace std;

const int S=2005,p=998244353;

int n,k;
int x[S],y[S];

inline int qpow(int x,int y=p-2)
{
	int res=1;
	for(;y>0;y>>=1,x=1ll*x*x%p) res=(y&1)?1ll*res*x%p:res;
	return res;
}

int main()
{
	scanf("%d%d",&n,&k);
	n--;
	for(int i=0;i<=n;i++) scanf("%d%d",&x[i],&y[i]);
	int ans=0;
	for(int i=0;i<=n;i++)
	{
		int s1=y[i],s2=1;
		for(int j=0;j<=n;j++) if(j!=i) s1=1ll*s1*(k-x[j])%p,s2=1ll*s2*(x[i]-x[j])%p;
		ans=(ans+1ll*s1*qpow(s2)%p+p)%p;
	}
	printf("%d\n",ans);
	return 0;
}

实际应用：求自然数幂和 $f_k(x)=\sum\limits_{i=1}^x i^k$

首先需要证明这个东西是多项式：

$\begin{aligned} f_k(x)&=\sum\limits_{i=1}^x i^k\\ &=1+\sum\limits_{i=1}^{x-1} (i+1)^k\\ &=1+\sum\limits_{i=1}^{x-1} \sum\limits_{j=0}^{k}\binom{k}{j}i^j\\ &=1+\sum\limits_{j=0}^{k}\binom{k}{j} \sum\limits_{i=1}^{x-1}i^j\\ &=\sum\limits_{i=1}^{x-1}i^k+1+\sum\limits_{j=0}^{k-1}\binom{k}{j} f_{j}(x-1)\\ \end{aligned}$

$\sum\limits_{i=1}^x i^k=\sum\limits_{i=1}^{x-1}i^k+1+\sum\limits_{j=0}^{k-1}\binom{k}{j} f_{j}(x-1)\\ x^k=1+\sum\limits_{j=0}^{k-1}\binom{k}{j} f_{j}(x-1)\\ x^k-\binom{k}{k-1}f_{k-1}(x-1)=1+\sum\limits_{j=0}^{k-2}\binom{k}{j} f_{j}(x-1)\\ \binom{k}{k-1}f_{k-1}(x-1)=x^k-1-\sum\limits_{j=0}^{k-2}\binom{k}{j} f_{j}(x-1)\\ f_{k-1}(x-1)=\frac{x^k-1-\sum\limits_{j=0}^{k-2}\binom{k}{j} f_{j}(x-1)}{\binom{k}{k-1}}\\$

观察到 $f_0(x)=x$ 是个 $1$ 次多项式，所以 $f_k(x)$ 是 $k+1$ 次多项式。

例题：

求 $\sum\limits_{i=1}^xi^k\mod 10^9+7$。

$1\le x\le 10^9,0\le k\le 10^6$。

直接拉插即可，但是直接拉插是 $O(k^2)$ 的，所以需要一些优化。

看回拉插式子：

$f(x)=\sum\limits_{i=0}^{k+1} y_i\frac{\prod\limits_{j\not=i}x-x_j}{\prod\limits_{j\not=i}x_i-x_j}$

观察到 $x=\{0,1,2,\dots,k+1\}$ 时，$\prod\limits_{j\not=i}x_i-x_j$ 是两个阶乘乘起来，$\prod\limits_{j\not=i}x-x_j$ 可以预处理前缀后缀积，所以可以优化到 $O(k\log k)$ 预处理，$O(k)$ 插值计算。

代码如下：

#include <iostream>
#include <cstdio>

using namespace std;

const int S=1000005,p=1000000007;

int n,k;
int y[S];
int fra[S],inv[S];
int pre[S],suf[S];

inline int qpow(int x,int y=p-2)
{
	int res=1;
	for(;y>0;y>>=1,x=1ll*x*x%p) res=(y&1)?1ll*res*x%p:res;
	return res;
}

inline void init(int k)
{
	y[0]=0;
	for(int i=1;i<=k+1;i++) y[i]=(y[i-1]+qpow(i,k))%p;
	fra[0]=1;
	for(int i=1;i<=k+1;i++) fra[i]=1ll*fra[i-1]*i%p;
	inv[k+1]=qpow(fra[k+1]);
	for(int i=k+1;i>=1;i--) inv[i-1]=1ll*inv[i]*i%p;
}

inline int calc(int n,int k)
{
	pre[0]=n;
	for(int i=1;i<=k+1;i++) pre[i]=1ll*pre[i-1]*(n-i+p)%p;
	suf[k+1]=(n-k-1+p)%p;
	for(int i=k;i>=0;i--) suf[i]=1ll*suf[i+1]*(n-i+p)%p;
	int res=0;
	for(int i=0;i<=k+1;i++)
	{
		int s1=1ll*(i==0?1:pre[i-1])*(i==k+1?1:suf[i+1])%p;
		int s2=1ll*inv[i]*((k+1-i&1)?p-inv[k+1-i]:inv[k+1-i])%p;
		res=(res+1ll*y[i]*s1%p*s2%p)%p;
	}
	return res;
}

int main()
{
	scanf("%d%d",&n,&k);
	init(k);
	printf("%d\n",calc(n,k));
	return 0;
}

例题2：（原题面极其垃圾，给出正确的题面）

小豆喜欢玩游戏，现在他在玩一个游戏遇到这样的场面，有 $n-m$ 个怪，血量为 $\{1,2,3,\dots,n\}-\{a_1,a_2,a_3,\dots,a_m\}$。小豆手里有无限张“亵渎”。亵渎的效果是对所有的怪造成 $1$ 点伤害，如果有怪死亡，则再次施放该法术（再次施放不消耗额外的“亵渎”）。我们认为血量为 $0$ 怪物死亡。

小豆使用一张“亵渎”会获得一定的分数，在使用一张“亵渎”之后，每一个被该“亵渎”造成伤害的怪会产生 $x^k$ 的分数，其中 $x$ 是造成伤害前怪的血量，$k$ 是杀死所有怪物所需的“亵渎”的张数。

请求出杀死所有怪物后小豆的分数，对 $10^9+7$ 取模。

$1\le m\le 50,1\le n\le 10^{13}$。

不难发现，血量为 $\{1,2,3,\dots,n\}$ 这样连续的怪只需要一张“亵渎”即可杀光，所以 $k=m+1$，并且杀光一段的得分为 $\sum\limits_{i=1}^n i^k=\sum\limits_{i=1}^n i^{m+1}$。

也就是说，若令 $a_0=0$，$a_1<a_2<a_3<\dots<a_m$ 即给它排个序，那么答案即为 $\sum\limits_{i=0}^m\sum\limits_{j=1}^{n-a_i}j^{m+1}-\sum\limits_{j=i+1}^{m} (a_j-a_i)^{m+1}$，那么直接拉插即可。

更强力的应用：拉格朗日插值优化 dp

例题：

对于给定的 $n$ 和 $k$，一个正整数序列 $a$ 是合法的当且仅当：

• $|a|=n$
• $1\le a_i\le k$
• $a_i$ 两两不同

$a$ 的权值定义为 $\prod\limits_{i=1}^n a_i$，请你求出所有不同的合法的 $a$ 的权值和，其中两个序列 $a,b$ 不同当且仅当有至少一个 $i$ 满足 $1\le i\le n$ 且 $a_i\not=b_i$。

$1\le n\le 500,1\le k\le 10^9$。

首先不难发现只要钦定序列是递增的，最后乘上 $n!$ 就行了。那么就有个朴素的 dp 做法，设 $dp_{i,j}$ 为填完 $a_{[1,i]}$，$1\le a_i\le j$ 的不同的 $a$ 的权值和，那么显然有转移 $dp_{i,j}=j\times dp_{i-1,j-1}+dp_{i,j-1}$。

但是这样做是 $O(nk)$ 的，所以考虑优化。

• 第一种方法：猜测法（无脑，感性）

不妨猜测 $dp_{i,j}$ 是一个关于 $j$ 的多项式，严谨证明见第二种方法。

那么直接拉插即可，插 $x$ 坐标连续的 $n$ 个点可以做到线性时间复杂度，所以尽可能多算几个 $dp_{n,\dots}$，插值即可。

毕竟拉插不难写，考场上遇到某一维特别大的 dp 都可以试试。

• 第二种方法：证明法（理性）

有一个结论：若 $f(x)$ 是 $n$ 次多项式，那么 $f(x)-f(x-1)$ 是 $n-1$ 次多项式。

证明如下：

设 $f(x)$ 的 $n$ 次项系数为 $a$，即 $f(x)=\dots+ax^n$。那么 $f(x)-f(x-1)$ 的 $n$ 次项为：

$\begin{aligned} &ax^n-a(x-1)^n\\ &=a(x^n-(x-1)^n)\\ &=a\left(x^n-\sum\limits_{i=0}^n\binom{n}{i}x^i(-1)^{n-i}\right)\\ &=a\left(x^n-x^n-\sum\limits_{i=0}^{n-1}\binom{n}{i}x^i(-1)^{n-i}\right)\\ &=a\sum\limits_{i=0}^{n-1}\binom{n}{i}x^i(-1)^{n-i} \end{aligned}$

这个东西是 $n-1$ 次的，得证。

设 $pd_{i,j}$ 为填完 $a_{[1,i]}$，$a_i=j$ 的所有合法序列的权值和，显然有转移 $pd_{i,j}=j\sum\limits_{k=1}^{j-1}pd_{i-1,k}$，那么因为 $pd_{1,j}$ 是关于 $j$ 的 $0$ 次多项式，并且以 $j$ 为上界的求和会使次数加一，乘上 $j$ 也会使次数加一，所以 $pd_{i,j}$ 是关于 $j$ 的 $2i$ 次多项式。

因为 $pd_{i,j}=dp_{i,j}-dp_{i,j-1}$，作差会使次数减 $1$，所以 $dp_{i,j}$ 是关于 $j$ 的 $2i+1$ 次多项式。

那么问题就变得简单起来了，$O(n^2)$ 求出 $dp_{n,[0,2n+1]}$ 后直接拉插即可。

代码如下：

#include <iostream>
#include <cstdio>

using namespace std;

const int S=1005;

int k,n,p;
int dp[S][S];

inline int qpow(int x,int y)
{
	int res=1;
	for(;y>0;y>>=1,x=1ll*x*x%p) res=(y&1)?1ll*res*x%p:res;
	return res;
}

int main()
{
	scanf("%d%d%d",&k,&n,&p);
	for(int i=0;i<=n*2+1;i++) dp[1][i]=(i+dp[1][i-1])%p;
	for(int i=2;i<=n;i++)
	{
		dp[i][0]=0;
		for(int j=1;j<=n*2+1;j++) dp[i][j]=(1ll*j*dp[i-1][j-1]%p+dp[i][j-1])%p;
	}
	int ans=0;
	for(int i=0;i<=n*2+1;i++)
	{
		int s1=1,s2=1;
		for(int j=0;j<=n*2+1;j++) if(j!=i) s1=1ll*s1*((k-j)%p+p)%p,s2=1ll*s2*((i-j)%p+p)%p;
		ans=(ans+1ll*dp[n][i]*s1%p*qpow(s2,p-2)%p)%p;
	}
	int fra=1;
	for(int i=1;i<=n;i++) fra=1ll*fra*i%p;
	printf("%d\n",1ll*ans*fra%p);
	return 0;
}

后记

场上遇到一些某维很大的式子，可以猜测它是关于这维的多项式，线性拉插尽可能多个点来尝试求解。