拉格朗日插值学习笔记
引入
众所周知，\(n+1\)?? 个点 \((x_i,y_i)\)?? 可以唯一地确定一个 \(n\)?? 次多项式，给你 \(n+1\)?? 个点，让你求出所对应的那个多项式，怎么做？
常规的套路是待定系数法列出方程组，然后用高斯消元 \(O(n^3)\)????? 求解，显然时间复杂度过大，且往往会有精度问题，还有什么更优的方法吗？有，拉格朗日插值，时间复杂度 \(O(n^2)\)?????，而且完全不会有精度问题。
说到拉格朗日插值，就不得不先提一下到底什么是插值了。简单来说，插值就是给你几个离散的点 \(\big(x,f(x)\big)\)?????????????，让你根据这些点求出一个函数，满足这些点都在这个函数上，称为插值函数；然后一般会给你几个新的 \(x\)??????，让你求出这些 \(x\)???? 所对应的 \(f(x)\)????????????? 是多少，称为内插。显然，根据不同的插值方法，插出来的结果也不尽相同，一般多应用于数值分析领域。
举个栗子，给出 \((1,1),~(2,3),~(3,5),~(4,7)\)，问你 \(f(5)\) 是多少。
答案是 \(114514\)??。
因为当这个函数是 \(\displaystyle f(x)=\frac {18111}2x^4-90555~x^3+\frac{633885}2x^2-452773~x+217331\) 时，有：
\[f(1)=1\\f(2)=3\\f(3)=5\\f(4)=7\\f(5)=114514 \]这插值怎么这么臭啊（恼）拉格朗日插值
目的是构造函数 \(F(x)\)??????，使得对于给定的点 \(x_i\)??????，都有 \(F(x_i)=y_i\)??????。怎么构造呢？我们可以让 \(F(x)=y_1f_1(x)+y_2f_2(x)+\dots+y_{n+1}f_{n+1}(x)\)??，?且 \(f_i(x)\) 只有当 \(x=x_i\) 时取值为 1，否则取值为 0，显然这样构造出来的 \(F(x)\) 是满足条件的。
那如何构造 \(f_i(x)\)?? 呢？简单，只需让 \(f_i(x)=\prod\limits_{j\neq i}\dfrac{x-x_j}{x_i-x_j}\)??? 即可。就拿上面的很臭的例子来说吧，我们有：
显然当 \(x\neq x_i\)?? 时，\(f_i(x)\)?? 的分子肯定有一项为 0，值为 0；当 \(x=x_i\) 时，\(f_i(x)\) 的分子分母相等，值为 1，即 \(f_i(x)\) 刚好满足上文所说的条件。
此时 \(F(x)=f_1(x)+3f_2(x)+5f_3(x)+7f_4(x)\)，试着把几个点值都带进去，发现都对的上114514无视即可，于是满足条件的插值函数 \(F(x)\)???? 就成功构造出来了~
应用&&例题
在算法竞赛中，往往会遇到这样一类问题，它们会直接或间接地给你几个点，然后让你求另一个点的值，这类问题往往会跟多项式问题挂上钩。现在我们知道了拉格朗日插值的公式，就可以利用它解决一些经典的插值问题了，这里上几个例题帮助说明。
【模板】拉格朗日插值 - 洛谷
题目大意
给出 \(n\)??? 个点和一个数 \(k\)???，让你求出根据这 \(n\)??? 个点确定的多项式 \(y=F(x)\)??? 在 \(k\)??? 上的值，即 \(F(k)\)???，答案对 998244353 取模。
\((1\leq2e3,~1\leq x_i,y_i,k<998244353),~x_i~两两不同\)
题目分析
模板题，根据拉格朗日插值，我们有：
\[F(x)=\sum_{i=1}^n\left(y_i\prod_{j=1,j\neq i}^n\frac{x-x_j}{x_i-x_j}\right) \]直接把 \(x=k\)???? 带进去暴力算就行，时间复杂度 \(O(n^2)\)????，记得分母先做连乘最后再求逆元，不然复杂度多个 \(\log\)???，只能得 60 分。我就是每次都对分母求逆元然后过不了
代码实现
#include
using namespace std;
const long long mod=998244353;
const int maxn=2e3+10;
long long x[maxn],y[maxn];
int n,k;
void exgcd(long long a,long long b,long long &x,long long &y){
if(b==0){
x=1,y=0;
return ;
}
exgcd(b,a%b,y,x);
y-=a/b*x;
}
long long inv(long long a){
long long x,y;
exgcd(a,mod,x,y);
return x;
}
int main(){
cin >> n >> k;
for(int i=1;i<=n;i++){
cin >> x[i] >> y[i];
}
long long ans=0;
for(int i=1;i<=n;i++){
long long fz=y[i];
long long fm=1;
for(int j=1;j<=n;j++){
if(i==j) continue;
long long t=k-x[j];
fm*=x[i]-x[j];
fm%=mod;
fz=fz*t%mod;
}
fm<0?fm+=mod:0;
fm=inv(fm);
ans+=fz*fm%mod;
ans%=mod;
}
ans<0?ans+=mod:0;
cout << ans;
}
CF622F The Sum of the k-th Powers - 洛谷
题目大意
给定 \(n\) 和 \(k\)，让你求 \(\sum_{i=1}^ni^k\mod {(1e9+7)}\)。
\((1\leq n\leq 1e9,~0\leq k\leq1e6)\)??
题目分析
乍一看好像跟拉格朗日插值没啥关系，但众所周知，次方数的前缀和是有公式的，而 \(k\)???? 次方数的前缀和公式一定是个 \(k+1\)????? 次多项式。证明的话，当然方法不止一种，有兴趣的同学可以自行去了解一下。
假设 \(k\)????????? 次方数的前缀和公式 \(F_k(n)=\sum_{i=1}^ni^k\)?????????，则我们可以计算出 \(k+2\)??????? 个点 \(\big(x,F_k(x)\big)\)???????，然后利用拉格朗日插值求得 \(F_k(n)\)????。
计算 \(k+2\)???????? 个点可以直接用快速幂 \(O(n\log n)\)???????? 求，也可以线性筛 \(x^k\)???????? 然后做前缀和 \(O(n)\)????? 求，问题不大，主要问题是拉格朗日插值的时间复杂度是 \(O(n^2)\)????? 的，\(\bf n\)????? 在这里是点的数量，即 \(k+2\)?????，直接套板子必 T。
注意到横坐标 \(x\)? 的值是连续的，所以可以简化运算至 \(O(n)\) 时间复杂度，推导如下：
\[F(n)=\sum_{i=1}^{k+2}\Bigg({\large y_i*f_i(x)}\Bigg)=\sum_{i=1}^{k+2}\left({\large y_i}\prod_{j=1,j\neq i}^{k+2}\frac{n-x_j}{x_i-x_j}\right)=\sum_{i=1}^{k+2}\left({\large y_i}\prod_{j=1,j\neq i}^{k+2}\frac{n-j}{i-j}\right) \]我们可以 \(O(1)\)? 求 \(f_i(x)=\displaystyle\prod\limits_{j=1,j\neq i}^{k+2}\frac{n-j}{i-j}\)?? ：
\[分子：\prod_{j=1,j\neq i}^{k+2}(n-j)=\prod_{j=1,j\neq i}^{i-1}(n-j)*\prod_{j=i+1,j\neq i}^{k+2}(n-j)\\ 分母：\prod_{j=1,j\neq i}^{k+2}(i-j)=\prod_{j=1,j\neq i}^{i-1}(i-j)*\prod_{j=i+1,j\neq i}^{k+2}(i-j)=(-1)^{k+2-i}*(i-1)!*(k+2-i)! \]显然可以 \(O(n)\)??? 预处理分子和分母，然后 \(O(1)\)??? 求解 \(y_i*f_i(x)\)???，从而将总时间复杂度降为了 \(O(n)\)??。
代码实现
#include
using namespace std;
const int mod=1e9+7;
const int maxn=1e6+6;
long long n,k;
long long prime[maxn];
bool isp[maxn];
long long y[maxn],f[maxn],invf[maxn],pre[maxn],suf[maxn];
long long qsm(long long a,long long b){
long long ret=1;
while(b){
if(b&1){
ret*=a,ret%=mod;
}
b>>=1,a*=a,a%=mod;
}
return ret;
}
void inip(){
f[0]=invf[0]=f[1]=invf[1]=1;
y[1]=1;
for(int i=2;i> n >> k;
inip();
pre[0]=1,suf[k+3]=1,invf[k+3]=qsm(f[k+3],mod-2);
for(int i=1;i<=k+2;i++){
pre[i]=pre[i-1]*(n-i)%mod;
y[i]+=y[i-1];
y[i]%=mod;
}
for(int i=k+2;i>=1;i--){
invf[i]=invf[i+1]*(i+1)%mod;
suf[i]=suf[i+1]*(n-i)%mod;
}
if(n<=k+2){
cout << y[n];
}
else {
long long ans=0;
for(int i=1;i<=k+2;i++){
long long fz=pre[i-1]*suf[i+1]%mod;
long long fm=invf[i-1]*((k+2-i&1)?mod-invf[k+2-i]:invf[k+2-i])%mod;
ans+=y[i]*fz%mod*fm%mod;
ans%=mod;
}
ans<0?ans+=mod:ans;
cout << ans;
}
}
参考资料
- 插值（离散数学名词）_百度百科 (baidu.com)
- 拉格朗日插值法(图文详解) - Angel_Kitty - 博客园 (cnblogs.com)
- 从拉插到快速插值求值 - command_block 的博客 - 洛谷博客 (luogu.com.cn)
- 题解 P4781 拉格朗日插值 - attack 的博客 - 洛谷博客 (luogu.com.cn)