拉格朗日插值法学习笔记

拉格朗日插值法

无关内容

• 无关前置：今天学习$\mathrm{L}\mathrm{a}\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}$，然后我一直疑惑为什么音译过来是拉格朗日（不应该是拉格阮籍吗？QAQ），最后在同学的帮助下，在知乎上找到了原因：原来真相居然是，法语，emmmm
  

进入正题

---

• 百度一下

• 拉格朗日公式

• 问题引入：

现在给你一个$n$次多项式的点值表示法（说简单点就是给你一个多项式函数图像上的$n+1$个点对$(x_i,y_i)$），请你求出当$x=k$时的函数值（也就是将$k$带入多项式求值）。

---

• 暴力想法：我们有了$n+1$个点对，那么将其分别带入可以得到$n+1$个$n$元一次（本来是一元$n$次方程，但是将$x^r$看作不同的未知数就变成了$n+1$元一次方程），然后高斯消元即可，复杂度$O(n^3)$。

• 进一步思考：但是对于$n$达到几千的情况，$O(n^3)$显然会超时，那么有没有更快的方法呢？

• 我们可以用斜率来拟合原函数，我们可以通过先枚举一个点，然后求给定的$x=k$点与其他点的$x_i$组成的两条直线的斜率之比:

假设插入的$k$求出点对为$(x,y)$，现在枚举的 $i,j(i̸=j)$ , 那么斜率之比为$\frac{\frac{x-x_j}{y-y_j}}{\frac{x_i-x_j}{y_i-y_j}}$，整理得$\frac{x-x_j}{x_i-x_j}\times \frac{y_i-y_j}{y-y_j}$，那么比例就大概为$\frac{x-x_j}{x_i-x_j}$（因为y并不知道）

然后用直线的斜率之比的乘积来当做新的斜率之比，然后我们就拟合出了一条$\prod _{i̸=j}(\frac{x-x_j}{x_i-x_j})$斜率的直线，那么该点对于给定的$x=k$的贡献（也就是对于的$y$值）可以由比例得知为$y_i\prod _{i̸=j}(\frac{x-x_j}{x_i-x_j})$，那么将$n+1$个点的所有贡献累加起来便得知了当$x=k$的时候$y=\sum _{i=0}^n{y}_i\prod _{i̸=j}(\frac{x-x_j}{x_i-x_j})$（$i̸=j$是因为自己和自己一个点是不能求斜率的）

所以，我们正式推出拉格朗日插值公式：
$$\sum _{i=0}^n{y}_i\prod _{i̸=j}(\frac{x-x_j}{x_i-x_j})$$

拉格朗日插值法可能会有一定的误差，因为是拟合出来的。

（证明与推导可能有些不妥或错误的地方，可以提出或者参考其他的blog或者百科）

此公式还有另外一个证明，就是将给出的$x_i$带入该公式，我们会发现，出除了当你枚举的$i$等于带入的$x_i$的$i$时，其他情况都有一个$x_i-x_i=0$，而当相等时，则有$y_i\prod _{i̸=j}(\frac{x_i-x_j}{x_i-x_j})$，我们会发现左边的值一定为$1$（分子分母相同）,那么最后插值的答案确实就等于$y_i$，所以是正确的。
推广来说，我们令${\beta }_i(x)=\prod _{i̸=j}\frac{x-x_j}{x_i-x_j}$，我们就可以得知函数${\beta }_i(x_i)=1$，而${\beta }_i(x_j)=0$

---

那么，我们通过公式可以看出，这个方法的复杂度是$O(n^2)$的，如果涉及取模和逆元，那么使用快速幂，复杂度则为$O(n_{log}val+n^2)$还是$O(n^2)$。

代码实现就非常简单啦！

丑陋代码如下:

点对为x[i],y[i]。
n-1次多项式，求插入k的值
#define ll long long
#define fpow(a,b) pow(a,b)
ll Lagrange(ll *x,ll *y,int n,ll k){
	ll s1,s2,ans=0;
	for(RG int i=1;i<=n;++i){
		s1=1;s2=1;
		for(RG int j=1;j<=n;++j){
			if(i==j) continue;
			s1=(s1*(k-x[j]+Mod)%Mod)%Mod;
			s2=(s2*(x[i]-x[j]+Mod)%Mod)%Mod;
		}
		ans=(ans+s1*fpow(s2,Mod-2)%Mod*y[i]%Mod)%Mod;
	}
	return ans;
}

• 特殊情况：

有时我们的$x_i$会是连续的取值，如$x_i$为$1\sim n$的连续正整数，这时再观察原式子，则变成了：
$$\sum _{i=0}^n{y}_i\prod _{i̸=j}\frac{x-j}{i-j}$$

于是我们对于分子预处理前缀积和后缀积，对于分母再预处理阶乘（取模还有阶乘的逆元），那么就$O(n+logMaxval+n)$相当于$O(n)$求出来了。

此时插入$x=k$:
前缀积: $pre[i]=\prod _{j=1}^i(k-j)$
后缀积: $suf[i]=\prod _{j=i}^n(k-j)$
那么，显然分子就等于$pre[i-1]\times suf[i+1]$（刚好取出第$i$个）
阶乘: $fac[i]=\prod _{j=1}^{i}j$
逆元可以这样处理:(这里使用费马小定理，模数为质数，你也可以使用其他定理求）)
$invfac[n]=pow(fac[n],mod-2)$
然后线性递推出所有的逆元，$invfac[i]=invfac[i+1]\times (i+1)$
那么，显然分母就等于$invfac[i-1]\times invfac[n-i]$

但是这里分母还有一个正负的问题，其实不难发现，当$n-i$为奇数时就是负数，否则为偶数。

那么代码也就显然啦，下面上丑陋代码QAQ

ll fac_inv[M];
ll pre[M],suf[M];
ll SpLagrange(/*ll *x,(x为1~n连续整数)*/ll *y,int n,ll k){
	pre[0]=1;for(RG int i=1;i<=n;++i)pre[i]=pre[i-1]*((k-i+Mod)%Mod)%Mod;
	suf[n+1]=1;for(RG int i=n;i>=1;--i)suf[i]=suf[i+1]*((k-i+Mod)%Mod)%Mod;
	ll fac=1;for(RG int i=2;i<=n;++i)fac=fac*i%Mod;
	fac=fpow(fac,Mod-2);fac_inv[n]=fac;
	for(RG int i=n-1;i>=1;--i)fac_inv[i]=fac_inv[i+1]*(i+1)%Mod;
	for(RG int i=1;i<=n;++i){
		ll s1=pre[i-1]*suf[i+1]%Mod;
		ll s2=inv_fac[i-1]*inv_fac[n-i]%Mod;
		ll flag=((n-i)&1)?-1:1;
		ans=(ans+flag*s1*s2%Mod*y[i]%Mod)%Mod;
	}
	return (ans+Mod)%Mod;
}

• 既然已经有了一些基础知识，那么我们来做几道例题吧：

自然数的幂的前缀和

---

求下列式子的值：
$$\sum _{i=0}^n{i}^k$$

数据范围：$1\le n\le 10^{15},0\le k\le 10^6$，$Mod=998244353$

---

这是拉格朗日插值法的经典题目$\mathrm{C}\mathrm{F}622\mathrm{F}\mathrm{T}\mathrm{h}\mathrm{e}\mathrm{S}\mathrm{u}\mathrm{m}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}{\mathrm{k}}^{\mathrm{t}\mathrm{h}}\mathrm{P}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}\mathrm{s}$，根据大佬说，也可以使用伯努利数 + 任意模数NTT（大佬文章【真正讲解伯努利数的文章】），但是蒟蒻并不会QAQ（要用到生成函数和NTT的知识），所以我们还是老老实实用拉格朗日插值法吧（而且复杂度还要小一些）。

首先，显然复杂度不能有$n$，（最多也只能$logn$），所以我们要从较小的$k$上分析。

通过经验，我们知道：
当$k=0$时，求值公式为$n$，一个一次多项式。
当$k=1$时，求值公式为$\frac{n\times (n+1)}{2}$，一个二次多项式。
当$k=2$时，求值公式为$\frac{n\times (n+1)\times (2n+1)}{6}$，一个三次多项式。
当$k=3$时，求值公式为${\left(\frac{n\times (n+1)}{2}\right)}^2$，一个四次多项式。
当$k=4$时，求值公式为$\frac{n(n+1)(2n+1)(3n^2+3n-1)}{30}$，一个五次多项式。
$\cdots$
由此，我们大胆猜测，${\sum }_{i=0}^n{i}^k$为一个$k+1$次多项式。

证明（方法来自$\mathrm{I}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{i}\mathrm{n}\mathrm{e}$大佬的blog）：

• 我们将两个$k+1$次多项式$(n+1{)}^{k+1}$与$n^{k+1}$相减，结合二项式定理展开公式，得到（$\left(\genfrac{}{}{0ex}{}{n}{m}\right)$表示组合数）：
  $$\begin{gathered} (n+1{)}^{k+1}-n^{k+1}=\sum _{i=0}^{k+1}(\genfrac{}{}{0ex}{}{k+1}{i})n^i-n^{k+1} \\ =\sum _{i=0}^k(\genfrac{}{}{0ex}{}{k+1}{i})n^i \end{gathered}$$

• 再将$n^{k+1}$减去$(n-1{)}^{k+1}$，同理得到：
  $$n^{k+1}-(n-1{)}^{k+1}=\sum _{i=0}^k(\genfrac{}{}{0ex}{}{k+1}{i})(n-1{)}^i$$

• 我们令$\zeta (x)=\sum _{i=0}^k\left(\genfrac{}{}{0ex}{}{k+1}{i}\right)(n-x{)}^i$，然后我们求出下列式子：
  $$\sum _{x=-1}^n\zeta (x)$$

化简得到$(n+1{)}^{k+1}=\sum _{i=0}^k(\genfrac{}{}{0ex}{}{k+1}{i})i^k$，那么显然${\sum }_{i=0}^n{i}^k$为一个$k+1$次多项式。

---

那么原式是一个$k+1$次的多项式，所以我们取$k+2$个值$x_i$，并计算出对应取值$y_i={\sum }_{i=1}^{x_i}{i}^k$，显然，我们发现当你的$x_i$取$[1,k+2]$中的正整数值时最简单，此时这个$x_i$也符合前面讲的特殊情况，所以我们可以先预处理求出$(x_i,y_i)$，然后带入下列式子:
$$\sum _{i=1}^n{i}^k=\sum _{i=1}^{k+2}{y}_i\prod _{j=1,j̸=i}^{k+2}\frac{n-j}{i-j}$$
最后的复杂度为$O(klogk)$的预处理和$O(k)$的插值。

---

下面提供几道例题：

BZOJ4559 别人的讲解
HDU4059 别人的讲解
HDU6239 别人的讲解

---

其他的较好的文章 ：

• 1.in
• 2.in
• 3.in
• 4.in