拉格朗日插值法学习笔记

拉格朗日插值法

无关内容

无关前置：今天学习 $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 别人的讲解

其他的较好的文章：