拉格朗日插值学习笔记
拉格朗日插值学习笔记
简介
在 2022 年联合省选 D1T2 吃了这个的亏，来补一补。虽说我压根就没看出来跟插值有啥关系。
我们都知道，\(n\) 个点的点值 \((x_i,y_i)\) 可以唯一确定一个 \(n-1\) 次多项式（为了叙述方便，本文中所有“\(k\) 次多项式”“\(k\) 次函数”的最高次项系数可以为 \(0\)）。
拉格朗日插值就是用来求这个多项式的。
例子
例如，我们已知四个点值 \((-1,1)(0,2)(0.5,1.375)(1,1)\)，要求过这四个点的三次函数 \(f\)。
当然，你可以直接待定系数用高斯消元解方程，但那是 \(\mathcal{O}(n^3)\) 的，拉格朗日插值可以在 \(\mathcal{O}(n^2)\) 内求解。
约瑟夫·拉格朗日认为这个函数可以用四个三次函数线性组合出来。
首先构造一个三次函数 \(f_1\)，在 \(x=-1\) 的取值为 \(1\)，但在其他三个点的取值为 \(0\)。
类似地，构造 \(f_{2,3,4}\) 依次在每个点取值为 \(1\)，在其他三个点取值为 \(0\)。
画到一张图里就是这样：
那么这几个函数有啥用呢？
容易发现，\(f(x)=y_1f_1(x)+y_2f_2(x)+y_3f_3(x)+y_4f_4(x)\)，把那四个点点值代入进去就可以知道。
现在问题就转化为了怎么求 \(f_{1,2,3,4}\)。
类似于一次函数的两点式：
我们可以构造函数：
于是就做完了。
推导
是时候写一下推导了。
对于 \(n\) 个点值求 \(n-1\) 次多项式的问题，我们先有点值 \((x_j,y_j)\)（\(1\le j\le n\)），设函数 \(f_i\)（\(1\le i\le n\)），它们是 \(n-1\) 次函数，且满足：
则 \(f_i\) 可以写成：
最终得到：
代码
题目是 洛谷 P4781 拉格朗日插值
套上面公式直接算即可。
如果你在计算每个 \(\dfrac{x-x_j}{x_i-x_j}\) 的时候都求一遍逆元，会导致复杂度变为 \(\mathcal{O}(n^2\log n)\)，多带一个逆元的 \(\log\)。为了让复杂度瓶颈不在逆元上，我们通常分开算分子分母，在每个函数算完后再进行有理数取模，这样的复杂度为 \(\mathcal{O}(n^2)\)。
//By: Luogu@rui_er(122461)
#include <bits/stdc++.h>
#define rep(x,y,z) for(int x=y;x<=z;x++)
#define per(x,y,z) for(int x=y;x>=z;x--)
#define debug printf("Running %s on line %d...\n",__FUNCTION__,__LINE__)
#define fileIO(s) do{freopen(s".in","r",stdin);freopen(s".out","w",stdout);}while(false)
using namespace std;
typedef long long ll;
const int N = 2e3+5, mod = 998244353;
int n, k, x[N], y[N], ans;
template<typename T> void chkmin(T& x, T y) {if(x > y) x = y;}
template<typename T> void chkmax(T& x, T y) {if(x < y) x = y;}
int qpow(int x, int y) {
int ans = 1;
for(;y;y>>=1,x=1LL*x*x%mod) if(y & 1) ans = 1LL * ans * x % mod;
return ans;
}
int inv(int x) {return qpow(x, mod-2);}
int main() {
scanf("%d%d", &n, &k);
rep(i, 1, n) scanf("%d%d", &x[i], &y[i]);
rep(i, 1, n) {
int p = y[i], q = 1;
rep(j, 1, n) {
if(i != j) {
p = 1LL * p * (k - x[j]) % mod;
q = 1LL * q * (x[i] - x[j]) % mod;
}
}
ans = (ans + 1LL * p * inv(q) % mod + mod) % mod;
}
printf("%d\n", ans);
return 0;
}
连续点值的插值
如果已知的点值是连续点的点值，我们可以做到 \(\mathcal{O}(n)\) 的插值。
有时候发现了一个函数是 \(n\) 次多项式，就求 \(n+1\) 个点值进去插值。为了省事，这里我们令 \(x_i=i\)（\(1\le i\le n+1\)）。注意这里 \(n\) 与上面意义不一样，是次数而不是点数。
我们有拉插公式：
代入 \(x_i=i\)：
考虑怎么快速求 \(\displaystyle\prod\limits_{j=1\\j\ne i}^{n+1}\frac{x-j}{i-j}\)。
这个东西的分子是:
分母的话把 \(i-j\) 累乘拆成两段阶乘，就是：
于是连续点值的插值公式：
预处理前后缀积、阶乘阶乘逆然后代这个式子的复杂度为 \(\mathcal{O}(n)\)。
例子是 CodeForces 622F The Sum of the k-th Powers。
可以发现答案是 \(k+1\) 次多项式，因此代 \(k+2\) 个点进去拉插。
本题的 \(i^k\) 可以线性筛，因此复杂度可以做到 \(\mathcal{O}(n)\)。
//By: Luogu@rui_er(122461)
#include <bits/stdc++.h>
#define rep(x,y,z) for(int x=y;x<=z;x++)
#define per(x,y,z) for(int x=y;x>=z;x--)
#define debug printf("Running %s on line %d...\n",__FUNCTION__,__LINE__)
#define fileIO(s) do{freopen(s".in","r",stdin);freopen(s".out","w",stdout);}while(false)
using namespace std;
typedef long long ll;
const int N = 1e6+5, mod = 1e9+7;
int n, k, tab[N], p[N], pcnt, f[N], pre[N], suf[N], fac[N], inv[N], ans;
template<typename T> void chkmin(T& x, T y) {if(x > y) x = y;}
template<typename T> void chkmax(T& x, T y) {if(x < y) x = y;}
int qpow(int x, int y) {
int ans = 1;
for(;y;y>>=1,x=1LL*x*x%mod) if(y & 1) ans = 1LL * ans * x % mod;
return ans;
}
void sieve(int lim) {
f[1] = 1;
rep(i, 2, lim) {
if(!tab[i]) {
p[++pcnt] = i;
f[i] = qpow(i, k);
}
for(int j=1;j<=pcnt&&1LL*i*p[j]<=lim;j++) {
tab[i*p[j]] = 1;
f[i*p[j]] = 1LL * f[i] * f[p[j]] % mod;
if(!(i % p[j])) break;
}
}
rep(i, 2, lim) f[i] = (f[i-1] + f[i]) % mod;
}
int main() {
scanf("%d%d", &n, &k);
sieve(k+2);
if(n <= k + 2) return printf("%d\n", f[n])&0;
pre[0] = suf[k+3] = 1;
rep(i, 1, k+2) pre[i] = 1LL * pre[i-1] * (n - i) % mod;
per(i, k+2, 1) suf[i] = 1LL * suf[i+1] * (n - i) % mod;
fac[0] = inv[0] = fac[1] = inv[1] = 1;
rep(i, 2, k+2) {
fac[i] = 1LL * fac[i-1] * i % mod;
inv[i] = 1LL * (mod - mod / i) * inv[mod%i] % mod;
}
rep(i, 2, k+2) inv[i] = 1LL * inv[i-1] * inv[i] % mod;
rep(i, 1, k+2) {
int P = 1LL * pre[i-1] * suf[i+1] % mod;
int Q = 1LL * inv[i-1] * inv[k+2-i] % mod;
int mul = ((k + 2 - i) & 1) ? -1 : 1;
ans = (ans + 1LL * (Q * mul + mod) % mod * P % mod * f[i] % mod) % mod;
}
printf("%d\n", ans);
return 0;
}
重心拉格朗日插值
不会。
