拉格朗日插值法学习笔记

拉格朗日插值法

无关内容

• 无关前置：今天学习$m Lagrange$，然后我一直疑惑为什么音译过来是拉格朗日（不应该是拉格阮籍吗？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 eq 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} imes frac{y_i-y_j}{y-y_j}$，那么比例就大概为$frac{x-x_j}{x_i-x_j}$（因为y并不知道）

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

所以，我们正式推出拉格朗日插值公式：
$sumlimits_{i=0}^{n}y_iprodlimits_{i eq j}(frac{x-x_j}{x_i-x_j})$

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

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

此公式还有另外一个证明，就是将给出的$x_i$带入该公式，我们会发现，出除了当你枚举的$i$等于带入的$x_i$的$i$时，其他情况都有一个$x_i-x_i=0$，而当相等时，则有$y_iprodlimits_{i eq j}(frac{x_i-x_j}{x_i-x_j})$，我们会发现左边的值一定为$1$（分子分母相同）,那么最后插值的答案确实就等于$y_i$，所以是正确的。
推广来说，我们令$eta_i(x)=prodlimits_{i eq j}frac{x-x_j}{x_i-x_j}$，我们就可以得知函数$eta_i(x_i)=1$，而$eta_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$为$1sim n$的连续正整数，这时再观察原式子，则变成了：
$sumlimits_{i=0}^ny_i prod limits_{i eq j} frac{x-j}{i-j}$

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

此时插入$x=k$:
前缀积: $pre[i]=prodlimits_{j=1}^i(k-i)$
后缀积: $suf[i]=prodlimits_{j=i}^n(k-i)$
那么，显然分子就等于$pre[i-1] imes suf[i+1]$（刚好取出第$i$个）
阶乘: $fac[i]=prodlimits_{j=1}^ij$
逆元可以这样处理:(这里使用费马小定理，模数为质数，你也可以使用其他定理求）)
$invfac[n]=pow(fac[n],mod-2)$
然后线性递推出所有的逆元，$invfac[i]=invfac[i+1] imes (i+1)$
那么，显然分母就等于$invfac[i-1] imes 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[i+1]%Mod;
		ll flag=((n-i)&1)?-1:1;
		ans=(ans+flag*s1*s2%Mod*y[i]%Mod)%Mod;
	}
	return (ans+Mod)%Mod;
}

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

自然数的幂的前缀和

---

求下列式子的值：
$sumlimits_{i=0}^ni^k$

数据范围：$1leq nleq10^{15},0leq kleq 10^6$，$Mod=998244353$

---

这是拉格朗日插值法的经典题目$m CF622F The Sum of the k^{th} Powers$，根据大佬说，也可以使用伯努利数 + 任意模数NTT（大佬文章【真正讲解伯努利数的文章】），但是蒟蒻并不会QAQ（要用到生成函数和NTT的知识），所以我们还是老老实实用拉格朗日插值法吧（而且复杂度还要小一些）。

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

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

证明（方法来自$m Imagine$大佬的blog）：

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

• 再将$n^{k+1}$减去$(n-1)^{k+1}$，同理得到：
  $n^{k+1}-(n-1)^{k+1}=sumlimits_{i=0}^k binom{k+1}{i}(n-1)^i$

• 我们令$zeta(x)=sumlimits_{i=0}^k binom{k+1}{i}(n-x)^i$，然后我们求出下列式子：
  $sum_{x=-1}^nzeta(x)$

化简得到$(n+1)^{k+1}=sumlimits_{i=0}^k binom{k+1}{i}i^k$，那么显然$sum_{i=0}^ni^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)$，然后带入下列式子:
$sumlimits_{i=1}^ni^k=sumlimits_{i=1}^{k+2}y_iprodlimits_{j=1,j eq i}^{k+2}frac{n-j}{i-j}$
最后的复杂度为$O(klogk)$的预处理和$O(k)$的插值。

---

下面提供几道例题：

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

---

其他的较好的文章 ：

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