拉格朗日插值
1 引入
众所周知，给定二维平面上的任意 \(n\) 个点，一定会有一个不超过 \(n-1\) 次的多项式与之对应。那么假如给出 \(n\) 个点，需要找到一个多项式 \(f(x)\) 使其对应，这个过程就是多项式插值。
而拉格朗日插值就是求解多项式插值的一种插值法。
2 拉格朗日插值
2.1 朴素解法
假设我们现在给出的 \(n\) 个点坐标为 \((x_i,y_i)\)，考虑构造出如下函数：
其中 \(g_i(x)\) 需要满足当 \(x\) 取值为 \(x_i\) 时值为 \(1\)，取值为其他 \(x_j(i\ne j)\) 时值为 \(0\)。考虑因式定理，当 \(x=x_j(i\ne j)\) 时该函数值为 \(0\)，说明 \(g_i(x)\) 有因式 \((x-x_j)(i\ne j)\)。也就是说，\(g_i(x)\) 有因式 \(\prod\limits_{j\ne i}(x-x_j)\)。
这个时候我们已经可以满足取值为其他 \(x_j(i\ne j)\) 时值为 \(0\) 这个条件了，考虑满足 \(x\) 取值为 \(x_i\) 时值为 \(1\) 这个条件。实际上这个很简单，我们给上面那个因式除以 \(\prod\limits_{j\ne i}(x_i-x_j)\) 即可，所以最后构造出的 \(g_i(x)\) 如下：
此时假如我们要查询一个单点 \(x\) 的值 \(f(x)\)，直接套用上面公式求解显然是 \(O(n^2)\) 的。
模板题：【模板】拉格朗日插值，代码如下：
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int Maxn = 2e5 + 5;
const int Inf = 2e9;
const int Mod = 998244353;
int n, k, x[Maxn], y[Maxn];
int qpow(int a, int b) {
int res = 1;
while(b) {
if(b & 1) res = res * a % Mod;
a = a * a % Mod;
b >>= 1;
}
return res;
}
int lagrange(int x[], int y[], int X) {
int ans = 0;
for(int i = 1; i <= n; i++) {
int res = 1;
for(int j = 1; j <= n; j++) {
if(i == j) continue;
res = res * (X - x[j]) % Mod * qpow(x[i] - x[j], Mod - 2) % Mod;
}
ans = (ans + res * y[i] % Mod + Mod) % Mod;
}
return ans;
}
signed main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n >> k;
for(int i = 1; i <= n; i++) {
cin >> x[i] >> y[i];
}
cout << lagrange(x, y, k) << '\n';
return 0;
}
2.2 线性求值解法
上面的朴素解法单次查询的复杂度是 \(O(n^2)\) 的，显然不优，因为我们一定可以将构造出来的式子写成一个多项式，单次求解的复杂度应该是 \(O(n)\) 的。所以问题在于怎么去求出多项式中每一项的系数。
考虑到上面构造出来的多项式本质上是：
考虑令 \(t_i=\dfrac{y_i}{\prod\limits_{j\ne i}(x_i-x_j)}\)，此时多项式转化为 \(f(x)=\sum\limits_{i=1}^n t_i\times \prod\limits_{j\ne i}(x-x_j)\)。考虑先求出 \(\prod\limits_{j=1}^n (x-x_j)\) 这个多项式展开后每一项的系数，然后对于每一个 \(x_i\) 将 \((x-x_i)\) 除掉即可。这样求出的多项式系数乘上 \(t_i\) 之后就是其对应的贡献了。
这里说一下怎样直接除以多项式 \((x-x_i)\)。
我们先看常数项 \(a_0\)，显然它的值就是 \(\prod x_i\)，想要实现除掉多项式 \((x-x_i)\) 的效果直接给 \(a_0\) 除以 \(-x_i\) 即可。
然后看一次项 \(a_1\)，我们同样直接给它除以一个 \(-x_i\)，但是此时会发现一个问题，即 \(a_1\) 中本身有一项为 \(\prod\limits_{j\ne i} x_j\)，我们需要先将这个值减去才能除以 \(-x_i\)。不过好巧不巧，我们发现此时的 \(a_0\) 就是这个值，所以新的 \(a_1\) 就是 \(a_1\) 减去 \(a_0\) 后再除以 \(-x_i\)。
对于后面的每一项都是同理，因此可以得出 \({a_i}'=\dfrac{a_i-{a_{i-1}}'}{-x_i}\)。
代码如下：
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int Maxn = 2e5 + 5;
const int Inf = 2e9;
const int Mod = 998244353;
int n, k, x[Maxn], y[Maxn];
int qpow(int a, int b) {
int res = 1;
while(b) {
if(b & 1) res = res * a % Mod;
a = a * a % Mod;
b >>= 1;
}
return res;
}
int t[Maxn], f[Maxn], g[Maxn], fg[Maxn];
void init(int x[], int y[]) {
for(int i = 1; i <= n; i++) {
t[i] = 1;
for(int j = 1; j <= n; j++) {
if(i == j) continue;
t[i] = t[i] * (x[i] - x[j]) % Mod;
}
t[i] = qpow(t[i], Mod - 2) * y[i] % Mod;//求出 t[i]
}
fg[1] = 1;
for(int i = 1; i <= n; i++) {
for(int j = n; j >= 1; j--) fg[j] = (fg[j - 1] + fg[j] * (Mod - x[i]) % Mod) % Mod;//先求出所有 (x - x[i]) 的积
}
for(int i = 1; i <= n; i++) {
int inv = qpow(Mod - x[i], Mod - 2);
for(int j = 1; j <= n; j++) {
g[j] = (fg[j] - g[j - 1]) * inv % Mod;//利用上面的式子
f[j] = (f[j] + g[j] * t[i] % Mod + Mod) % Mod;//乘上 t[i] 之后贡献给对应系数
}
}
}
int lagrange(int X) {
int ans = 0, res = 1;
for(int i = 1; i <= n; i++) {//多项式求值
ans = (ans + res * f[i] % Mod) % Mod;
res = res * X % Mod;
}
return ans;
}
signed main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n >> k;
for(int i = 1; i <= n; i++) {
cin >> x[i] >> y[i];
}
init(x, y);
cout << lagrange(k) << '\n';
return 0;
}
2.3 横坐标连续的插值
在某些题目中，我们获得的 \(x_i\) 坐标是连续的整数，这个时候我们有一种 \(O(n)\) 预处理，\(O(n)\) 求解的做法。
考虑回到构造出的函数 \(f(x)=\sum\limits_{i=1}^ny_i\times\dfrac{\prod\limits_{j\ne i}(x-x_j)}{\prod\limits_{j\ne i}(x_i-x_j)}\)，此时它应该变成 \(f(x)=\sum\limits_{i=1}^ny_i\times\dfrac{\prod\limits_{j\ne i}(x-j)}{\prod\limits_{j\ne i}(i-j)}\)。考虑求出 \((x-i)\) 的前缀积和后缀积，分子就应该变成 \(pre_{i-1}\times suf_{i+1}\)。然后再看分母，发现它其实就是两个阶乘而已，所以分母可以写作 \((i-1)!\times(n-i)!\times(-1)^{n-i}\)。因此插值公式就可以写作：
此时即可 \(O(n)\) 求解。
代码如下：
int y[Maxn], f[Maxn], g[Maxn], pre[Maxn], suf[Maxn];
int qpow(int a, int b) {
int res = 1;
while(b) {
if(b & 1) res = res * a % Mod;
a = a * a % Mod;
b >>= 1;
}
return res;
}
void init(int n) {
f[0] = 1;
for(int i = 1; i <= n; i++) f[i] = f[i - 1] * i % Mod;
g[n] = qpow(f[n], Mod - 2);
for(int i = n - 1; i >= 0; i--) g[i] = g[i + 1] * (i + 1) % Mod;
}
int lagrange(int n, int X) {
pre[0] = suf[n + 1] = 1;
for(int i = 1; i <= n; i++) pre[i] = pre[i - 1] * (X - i) % Mod;
for(int i = n; i >= 1; i--) suf[i] = suf[i + 1] * (X - i) % Mod;
int res = 0;
for(int i = 1; i <= n; i++) {
int num = y[i] * pre[i - 1] % Mod * suf[i + 1] % Mod * g[i - 1] % Mod * g[n - i] % Mod;
if((n - i) & 1) res = (res - num + Mod) % Mod;
else res = (res + num) % Mod;
}
return res;
}
下面来看几道经典例题。
3 例题
例 1： [湖北省队互测2014] 一个人的数论
题意： 令 \(f_d(n)\) 表示所有小于 \(n\) 且与 \(n\) 互质的正整数的 \(d\) 次方之和，求 \(f_d(n)\)。
发现 \(n\) 过于巨大，因此我们想要构造出一个积性函数，这样就可以拆成每一个质因子去求解了。首先进行如下变换：
发现整个式子是一个狄利克雷卷积的形式，但是后面自然数的幂次方和并不是一个积性函数。这个时候运用拉格朗日插值，熟知自然数的 \(k\) 次方和可以写作一个 \(k+1\) 次多项式，所以只需要给定 \(k+2\) 个点值即可求出该多项式。
于是设该多项式为 \(g(x)=\sum\limits_{i=0}^{k+1} a_ix^i\)，则可得：
这个时候发现后面的和式就是一个积性函数了。由于 \(\mu(d)\) 在 \(d\) 有平方因子的时候值为 \(0\)，所以 \(d\) 的每种质因子都应只出现一次。于是令 \(x=\prod p_i\)，则后面和式的条件应该化为 \(d\mid x\)。设 \(F(n)=\sum\limits_{d\mid n}\mu(d)\cdot d^{k-i}\)，则有：
直接求解该式子即可，复杂度 \(O(k^2+kw)\)。
通过上例发现，拉格朗日插值的一个经典用途就是求出自然数的幂次方和。
例 2： [JLOI2016] 成绩比较
我们考虑 dp。由于最后一定会有 \(k\) 个人被碾压，所以剩下的 \(n-1-k\) 个人都至少有一科考的比 B 高。所以设 \(dp(i,j)\) 表示枚举到第 \(i\) 科，还有 \(j\) 个人没有满足 “至少有一科考的比 B 高” 的条件。枚举当前有多少个人考的比 B 高，则可得如下转移方程：
首先从 \(j+x\) 个人中任选 \(x\) 个，让他们考的比 B 高，则这些人就满足了要求。此时我们会有 \(n-1-k-j-x\) 已经满足了要求的人，从中任选几个填满 \(r_i-1\) 即可。
发现这个式子的后半部分无法直接求解，因为 \(u_i\le 10^9\)。但此时不难发现后面的和式中，每一项都是一个关于 \(u_i\) 的 \(n-1\) 次多项式，求总和后就是一个关于 \(u_i\) 的 \(n\) 次多项式。所以预处理出前 \(n+1\) 个点值，然后拉格朗日插值求解即可。复杂度 \(O(n^3)\)。
例 3： [集训队互测 2012] calc
首先将有序的 \(a\) 转化成选出一个无序的集合，算出总权值后乘上 \(n!\) 即为最后答案。
考虑对于无序集合怎么做，不难想到设 \(dp(i,j)\) 表示在 \([1,j]\) 中选了 \(i\) 个数的权值总和，则有转移方程：
发现此时直接暴力做是 \(O(kn)\) 的，肯定无法通过。此时我们将 \(i\) 看作常数，\(j\) 看作变量，\(dp(i,j)\) 实际上可以看作是一个关于 \(j\) 的多项式。那么现在的问题就在于这个多项式的次数，于是设 \(dp(i,j)\) 是一个 \(g(i)\) 次多项式。根据上面的转移式，可以得到 \(dp(i,j)-dp(i,j-1)=dp(i-1,j-1)\times i\)。
左边是多项式差分的形式，所以差分过后应该会得到一个 \(g(i)-1\) 次多项式，而右边本来是一个 \(g(i-1)\) 次多项式，再乘上 \(i\) 后就变为了 \(g(i-1)+1\) 次多项式，所以会有如下关系式：
所以可以得到 \(g(i)=g(i-1)+2\) 这个递推式。又由于 \(dp(0,0)=1\) 导出 \(g(0)=0\)，所以 \(g(i)=2i\)。
也就是说，\(dp(i,j)\) 是一个关于 \(j\) 的 \(2i\) 次多项式，只需要知道 \(2i+1\) 个点值即可拉格朗日插值求出 \(dp(n, k)\)。
