largrange 插值学习笔记

引入

给出 $n$ 个点 $(x_i,y_i)$，求一个 $n-1$ 次的多项式 $G(x)$，满足 $G(x_i)=y_i$，求出 $G(k)$ 的值。

largrange 插值

根据待定系数法的推论，若已知 $n$ 个点，就可确定一个 $n-1$ 次的多项式。而通过点的坐标来求多项式的过程叫做 插值。

显然，问题就是让我们求过点 $(x_i,y_i)$ 的 $n-1$ 次多项式。

当然可以考虑高斯消元，也就是计算机模拟待定系数法。但是这样的复杂度是 $\Theta (n^3)$。
当待定系数的时间复杂度无法满足需求时，我们就可以使用 $\mathbf{l}\mathbf{a}\mathbf{r}\mathbf{g}\mathbf{r}\mathbf{a}\mathbf{n}\mathbf{g}\mathbf{e}$ 插值（拉格朗日插值）

我们考虑对单点进行求解，构造一组多项式 $F_i(x),i\in [1,n]$，使得 $F_i(x)$ 满足以下条件：

1. $F_i(x_i)=y_i$
2. $F_i(x_j)=0,j\ne i$

若拥有这组多项式，那么最终的 $G(x)=\sum _{i=1}^n{F}_i(x)$，显然满足所要求的条件。

现在考虑构造 $F_i(x)$,若要满足条件 $2$，则可考虑利用因式定理，使 $F_i(x)=\prod _{j=1,j\ne i}^n(x-x_j)$。

此时，对于 $F_i(x_k)=\prod _{j=1,j\ne i}^n(x_k-x_j),k\ne i$，显然都存在一个 $x_j$ 使得 $(x_k-x_j)=0$，即 $F_i(x_k)=0$，满足条件 $2$。

此时我们将 $x_i$ 代入 $F_i(x)$，即 $F_i(x_i)=\prod _{j=1,j\ne i}^n(x_i-x_j)$，显然这是一个非 $0$ 的常数。（由于 $x_i$ 肯定互不相等，所以不可能为 $0$）

我们考虑使 $F_i(x_i)$ 的值变为 $y_i$，那么可以有这么一种思路，将多项式的值除以 $F_i(x_i)$ ，使 $F_i(x_i)=1$，然后乘上 $y_i$，即 $
F_i(x)=\dfrac{y_i\cdot\prod\limits_{j=1,j\ne i}^{n}(x-x_j)}{\prod\limits_{j=1,j\ne i}^{n}(x_i-x_j)}$。

显然，这样的 $F_i(x)$ 是同时满足上述两个条件的，因此

$$\begin{array}{rl}G(x)& =\sum _{i=1}^n{F}_i(x)\\ & =\sum _{i=1}^n(y_i\cdot \frac{\prod _{j=1,j\ne i}^n(x-x_j)}{\prod _{j=1,j\ne i}^n(x_i-x_j)})\\ & \\ & =\sum _{i=1}^n(y_i\cdot \prod _{j=1,j\ne i}^n\frac{x-x_j}{x_i-x_j})\end{array}$$

此时，我们就得到了所求的 $G(x)$，通过这种方法求得的它，有个名字，叫做拉格朗日多项式

或许你会觉得这个多项式长得有点怪，那是因为他带个 $\sum$ 的原因，如果将求和符号打开，还是可以转换为有 $n$ 个系数，我们所熟知的那种多项式，这种方法我们后面也会提到，但是这里仅需求出 $G(k)$ 的值，这里就不做过多讨论了。

代码实现

在普通代码中，显然，对于每一层 $i$ 的循环，我们都要处理出分子 $\prod _{j=1,j\ne i}^n(k-x_j)$ 和分母 $\prod _{j=1,j\ne i}^n(x_i-x_j)$，然后两式相除再照着结论打，就可以求出 $G(k)$ 的值。

#include <bits/stdc++.h>
#define int long long
using namespace std;
int mod=998244353;
int n,k,x,y;
struct node
{
	int x,y;
} pos[2010];
int ksm(int a,int b,int mod)
{
	int sum=1;
	while(b)
	{
		if(b&1) sum=sum*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return sum%mod;
}
int largrange(int s)
{
	int ans=0;
	for(int i=1;i<=n;i++)
	{
		int sum=pos[i].y,t=1;
		for(int j=1;j<=n;j++)
		{
			if(i==j) continue;
			sum=1ll*sum*(s-pos[j].x)%mod;  //处理分子
			t=1ll*t*(pos[i].x-pos[j].x)%mod;  //处理分母
		}
		ans=((ans%mod)+(sum*ksm(t,mod-2,mod)%mod))%mod;  //求和统计答案
	}
	return (ans%mod+mod)%mod;
}
signed main()
{
	ios::sync_with_stdio(false);
	cin.tie(0);
	cin>>n>>k;
	for(int i=1;i<=n;i++) cin>>pos[i].x>>pos[i].y;
	cout<<largrange(k);
	return 0;
}

经典例题

引入：CF622F

一句话题意：给定 $n,k$，求 $(\sum _{i=1}^n{i}^k)\mathrm{m}\mathrm{o}\mathrm{d}(1e9+7)$

解题思路

一般的拉格朗日插值问题的首要步骤就是先证明所求解问题最终结果是一个多项式，这样子我们就可以将这个多项式插出来了。

可以证明，数列 $S_n=\sum _{i=1}^n{i}^k$ 的通项是一个多项式。

证明：

差分可证，已经有大佬做过很好的证明了，这里直接贴贴（%%%）

---

largrange 插值的优化

知道了结果的通项一定是一个 $k+1$ 次的多项式，那么我们就可以直接用拉格朗日插值将这个通项插出来。

我们需要手动计算 $k+2$ 个点的值来确定最终答案的多项式。
考虑 $x_i$ 表示我们取得点，用 $y_i=\sum _{i=1}^{x_i}{i}^k$ 表示我们取得点的值，$G(x)$ 表示通项式，也就是求 $G(n)$。直接套入拉格朗日公式，此题的通项即为：

$$G(n)=\sum _{i=1}^{k+2}(y_i\prod _{j=1,j\ne i}^{k+2}\frac{n-x_j}{x_i-x_j})$$

但是，$k\le 10^6$，$\Theta (n^2)$ 的普通拉格朗日插值是跑不过去的。
考虑到在这道题中，我们可以自行取点来插值，即点对的横坐标 $x_i$ 是连续的，也就是说 $x_i=i$

我们就可以考虑将 $x_i=i$ 这个特殊条件代入，对基本的拉格朗日公式进行变形，即：

$$G(x)=\sum _{i=1}^{k+2}(y_i\prod _{j=1,j\ne i}^{k+2}\frac{n-j}{i-j})$$

现在我们需要考虑的就是，对于 $\prod _{j=1,j\ne i}^{k+2}\frac{n-j}{i-j}$，真的有必要每次都 $\Theta (n)$ 来求吗？

答案显然是否定的。

我们可以构造前缀积数组 $mul$，使得 $mu{l}_i=\prod _{j=1}^i(n-j)$，以及后缀积数组 $mul\alpha$ 使得 $mul{\alpha }_i=\prod _{j=i}^{k+2}(n-j)$

那么对于 $\prod _{j=1,j\ne i}^{k+2}\frac{n-j}{i-j}$ 我们就可以通过每次 $mu{l}_{i-1}\cdot mul{\alpha }_{i+1}$ 得到，复杂度 $\Theta (1)$

$mul,mul\alpha$ 可以 $\Theta (n)$ 预处理。

同时，我们再构造数组 $chu$，使得 $ch{u}_i=i!$,以及数组 $chu\alpha$，使得 $chu{\alpha }_i=\prod _{j=1}^i-j=i!(-1{)}^i$

那么对于模意义上的 $\frac{1}{\prod _{j=1,j\ne i}^{k+2}i-j}$，只需要求出 $ch{u}_{i-1}\cdot chu{\alpha }_{n-i}$ 再求逆元即可，时间总共 $\Theta (\log n)$，如果线性求逆元还可以优化到 $\Theta (1)$

$chu,chu\alpha$ 我们也可以 $\Theta (n)$ 预处理得到。

这样我们通过种种预处理，就将 $x_i$ 连续取值时的拉格朗日插值优化到了 $\Theta (n)$ 的时间复杂度。

奉上代码：

#include <bits/stdc++.h>
#define int long long
using namespace std;
int n,k,ans=0;
int mod=1e9+7;
struct node
{
	int x,y;
} pos[1000010];
int mul[1000010],chu[1000010],chu1[1000010],mul1[1000010];
int ksm(int a,int b)
{
	int sum=1;
	while(b)
	{
		if(b&1) sum=1ll*sum*a%mod;
		b>>=1;
		a=1ll*a*a%mod;
	}
	return (sum%mod+mod)%mod;
}
void largrange()
{
	for(int i=1;i<=k+2;i++)
	{
		int a=pos[i].y*mul[i-1]%mod*mul1[i+1]%mod;
		int b=chu[i-1]*chu1[k+2-i]%mod;
		b=ksm(b,mod-2);
		ans=1ll*(ans+(a*b%mod))%mod;
	}
	return ;
}
signed main()
{
	ios::sync_with_stdio(false);
	cin.tie(0);
	cin>>n>>k;
	for(int i=1;i<=k+2;i++)
	{
		pos[i].x=i;
		pos[i].y=(pos[i-1].y+ksm(i,k))%mod;
		if(i==n)
		{
			cout<<pos[i].y;
			return 0;
		}
	}
	mul[0]=mul1[k+3]=chu[0]=chu1[0]=1;
	for(int i=1;i<=k+2;i++) mul[i]=mul[i-1]*(n-pos[i].x)%mod;
	for(int i=k+2;i>=1;i--) mul1[i]=mul1[i+1]*(n-pos[i].x)%mod;
	for(int i=1;i<=k+2;i++) chu[i]=chu[i-1]*i%mod;
	for(int i=1;i<=k+2;i++) chu1[i]=chu1[i-1]*(-i)%mod;
	largrange();
	cout<<ans;
	return 0;
} 

ps:此份代码用了费马小定理求逆元，时间复杂度 $\Theta (n\log n)$，如果用线性求逆元可达到 $\Theta (n)$，都能通过此题。

---