【数学】多项式插值 - 拉格朗日插值
多项式插值简介
这部分内容收录在这篇文章中：【数学】多项式插值。
原理
由这 \(n + 1\) 个点，可以构造对应的 \(n + 1\) 个多项式，其中第 \(i\) 的多项式的形式为：
把下标为 \(i\) 的点 \(x_i\) 代入，可以得到：
把给定的点中，下标不为 \(i\) 的某个点 \(x_{j0}\) 代入，可以得到：
因为分子部分总存在 $ (x_{j0} - x_{j0}) = 0 $ 所以上式为 \(0\) 。
把不在给定的点中的任意点 \(x\) 代入，显然有 \(n\) 项，所以，这是一个 \(n\) 次的多项式。
于是，把上面的 \(n + 1\) 个不同的 \(l_i(x)\) 和 \(y_i\) 相乘，再相加，即可构造下面的多项式：
由上面的推理可知， $ L(x_i) = y_i $ ，并且 \(L(x_i)\) 是一个不超过 \(n\) 次的多项式（由于最高次的系数可能会被抵消，导致次数降低）。
由唯一性定理可知，这个 \(L(x)\) 就是我们要找的 \(P(x)\) ，即
模板
在代码中，为了统一我所有的模板和我一直以来的习惯，也就是从1开始计数的写法以及用 \(n\) 表示点数，代码部分中的下标和其他文章（包括前述内容）中或许有差异。
本文中所有的代码都用 \(n\) 表示点数，他们分别为 \((x_1, y_1), (x_2, y_2), (x_3, y_3), \cdots, (x_{n}, y_{n})\)，也就是多项式的最高次数不超过 \(n - 1\) 次。对一个不超过 \(n - 1\) 次的多项式插值，至少需要 \(n\) 个点。
对于任意点求解
直接求解
复杂度
时间复杂度 \(O(n^2)\) ，空间复杂度 \(O(1)\)
使用说明
- 检查数据规模
MAXN，注意数据规模指的是**点数 \(n\) **（或者多项式的次数 \(n - 1\) ）和模数MOD。 - 输入点数 \(n\) ，和两个1-index的数组 \(x, y\) 表示给定的点。
- 无需初始化，直接调用
ll lagrange_interpolation (ll x0, ll *x, ll *y, int n)，返回值即 \(L(x_0)\)。
代码
namespace LagrangeInterpolation {
static const int MOD = 998244353;
ll qpow (ll x, ll n) {
ll res = 1;
while (n) {
if (n & 1) {
res = res * x % MOD;
}
x = x * x % MOD;
n >>= 1;
}
return res;
}
ll inv (ll x) {
return qpow (x, MOD - 2);
}
/**
* 用 n 个点 (x[1], y[1]), (x[2], y[2]), ..., (x[n], y[n])
* 插值得到 n - 1 次多项式 L，直接求解 L(x0) 的值
*
* 时间复杂度: O(n^2)
* 空间复杂度：O(1)
*/
ll lagrange_interpolation (ll x0, ll *x, ll *y, int n) {
for (int i = 1; i <= n; ++i) {
if (x0 == x[i]) {
return y[i];
}
}
ll Lx0 = 0;
ll P = 1;
for (int i = 1; i <= n; ++i) {
P = P * (x0 - x[i] + MOD) % MOD;
}
for (int i = 1; i <= n; ++i) {
ll q = 1;
for (int j = 1; j <= n; ++j) {
if (j == i) {
continue;
}
q = q * (x[i] - x[j] + MOD) % MOD;
}
ll p = P * inv (x0 - x[i] + MOD) % MOD;
ll lix0 = p * inv (q) % MOD;
Lx0 = (Lx0 + lix0 * y[i] % MOD) % MOD;
}
return Lx0;
}
};
转换成多项式的系数表达
原理
当 \(x \neq x_i\) 恒成立时：
记 \(P = \prod\limits_{j = 0}^{n} (x - x_j)\) 这是一个和 \(i\) 无关的量，结果是一个 \(n\) 次多项式
记 \(C = \prod\limits_{j = 0, \; j \neq i}^{n} (x_i - x_j)\) 这是一个常数
记 \(Q_i = (x - x_i)\) 这是一个二项式
得到
预处理多项式 \(P\) ，然后用多项式除法除以 \(Q_i\) 再每一项除以常数 \(C_i\) 即可得到所需的系数，最后全部加起来得到 \(L(x)\) 的系数。
复杂度
预处理系数：时间复杂度 \(O(n^2)\) ，空间复杂度 \(O(n)\)
单次求值：时间复杂度 \(O(n)\) ，空间复杂度 \(O(1)\)
对横坐标为连续整数的点求解
根据 \(x[] = 1, 2, 3, \cdots, k\) 时的值 \(y[]\) ，
插值计算不超过 \(k-1\) 次的多项式 \(L\) 在 \(x_0\) 处的取值 \(L(x_0)\) 。
时间复杂度 \(O(k)\) ，空间复杂度 \(O(MAXK)\) 。
\(P1[i]\) 是 \(x - i\) 的前缀积，
\(P2[i]\) 是 \(x - i\) 的后缀积，
\(Q[i]\) 是阶乘 \(i!\) 的逆元。
注意：
1. \(k\) 为点的数量， \(MAXK\) 是否足够大？
2. \(y[]\) 是否传入负数？
3. \(int\) 是否可能溢出？
4. \(MOD\) 必须是大于其他所有数字的质数，否则逆元会不存在。
卡常：
1. \(q\) 只与 \(k\) 有关，对于固定的 \(k\) ，可以预处理后多次使用，但是真没什么必要。
使用：
1. 直接调用 int lagrange (int x0, int *y, int k) 即可，其会自动初始化需要的数组，并且可以多次调用。
namespace Lagrange1 {
const int MAXK = 1e6 + 5;
static int P1[MAXK], P2[MAXK], Q[MAXK];
ll qpow (ll x, ll n) {
ll res = 1;
for (; n; n >>= 1) {
if (n & 1) res = res * x % MOD;
x = x * x % MOD;
}
return res;
}
void init_P (int x0, int k) {
P1[0] = 1;
for (int i = 1; i <= k; ++i)
P1[i] = 1LL * P1[i - 1] * (x0 - i) % MOD;
P2[k + 1] = 1;
for (int i = k; i >= 1; --i)
P2[i] = 1LL * P2[i + 1] * (x0 - i) % MOD;
}
void init_Q (int k) {
if (Q[k] != 0) return;
Q[k] = 1;
for (int i = 1; i <= k; ++i)
Q[k] = 1LL * Q[k] * i % MOD;
Q[k] = qpow (Q[k], MOD - 2);
for (int i = k; i >= 1; --i)
Q[i - 1] = 1LL * Q[i] * i % MOD;
}
int lagrange (int x0, int *y, int k) {
if (x0 <= k) return y[x0];
init_P (x0, k);
init_Q (k);
int Lx0 = 0;
for (int i = 1; i <= k; ++i) {
ll p = 1LL * P1[i - 1] * P2[i + 1] % MOD;
ll q = 1LL * Q[i - 1] * Q[k - i] % MOD;
if ( (k - i) & 1) q = MOD - q;
Lx0 = (Lx0 + p * q % MOD * y[i] % MOD) % MOD;
}
return Lx0;
}
};
using Lagrange1::lagrange;
对横坐标为连续整数的点求解（不使用额外空间）
时间复杂度 \(O(k\log{M})\) ，空间复杂度 \(O(1)\) 。
使用：
1. 直接调用 int lagrange (int x0, int *y, int k) 即可，并且可以多次调用。
namespace Lagrange2 {
ll qpow (ll x, ll n) {
ll res = 1;
for (; n; n >>= 1) {
if (n & 1) res = res * x % MOD;
x = x * x % MOD;
}
return res;
}
ll inv (ll x) {
return qpow (x, MOD - 2);
}
int lagrange (int x0, int *y, int k) {
if (x0 <= k)
return y[x0];
ll Lx0 = 0, P = 1, Q1 = 1, Q2 = 1;
for (int i = 1; i <= k; ++i) {
P = P * (x0 - i) % MOD;
if (i < k) Q2 = Q2 * (MOD - i) % MOD;
}
for (int i = 1; i <= k; ++i) {
ll p = P * inv (x0 - i) % MOD;
ll q = inv (Q1 * Q2 % MOD);
Lx0 = (Lx0 + p * q % MOD * y[i]) % MOD;
Q1 = Q1 * i % MOD;
Q2 = Q2 * inv (MOD - k + i) % MOD;
}
return Lx0;
}
};
using Lagrange2::lagrange;
拓展
选取等差数列进行插值
例如：求 \(\sum\limits_{i=1}^ni^{k}\) 的值，众所周知，这是一个 \(k+1\) 次多项式，至少需要 \(k+2\) 个点。
这个值有很多种求法，有 \(O(n\log k)\) 的朴素解法（枚举+快速幂），有 \(O(k^2)\) 的系数递推法，使用上述的适用范围广泛的拉格朗日插值法，也可以做到 \(O(k^2)\) ，算法的瓶颈在于计算插值多项式的分母部分。
这里可以选取一些特殊点来加速计算插值多项式的分母部分。例如通过选取一个等差数列，那么这个分母的临项会变得非常有规律。简单起见可以取正整数。
观察分子部分，对每个 \(i\) 来说，分子部分乘上 \((x-x_i)\) 后，都变为公共的 \(\prod\limits_{j=0}^{k}(x-x_j)\) ，那么在使用的时候再把这个 \((x-x_i)\) 除掉就行，复杂度为 \(O(k)\) 预处理，\(O(\log M)\) 单次使用，复杂度为 \(O(k\log M)\)。算法的复杂度瓶颈在于计算分母部分。
观察分母部分，在选取 \(x_i=i\) 后，分母变为两个自然数的前缀积之积，符号是正负交替。更准确地说，第 \(i\) 项的值是： \([\prod\limits_{j=0}^{i-1}(i-j)][\prod\limits_{j=i+1}^{k}(i-j)]=(\prod\limits_{j=1}^{i}j)(\prod\limits_{j=1}^{k-i}j)(-1)^{k-i}\)
由于算法竞赛一般是在模 \(M\) 意义下进行，算上求解逆元时间复杂度应为 \(O(k\log{M})\) 。
注： \(\sum_{i=1}^ni^{k}\) 的暴力求解其实可以是 \(O(n)\) 的，是使用线性筛去筛出每个数的 \(k\) 次方，然后线性递推出前缀和，不过这实际上是在画蛇添足，因为快速幂的常数实际上很小，而且 \(k\) 本身不大。
卡常：对同一个多项式多次求值时，插值完成后可以保存其分母（毒瘤）。
重心拉格朗日插值法
观察上面这个式子：
若记 \(g=\prod\limits_{j=1}^{k} (x_0-x_j)\) ，代入上式得到：
若记 \(t_i=\prod\limits_{j=1}^{k} [i\neq j]\frac{y_i}{(x_i-x_j)}\) ，代入上式得到
那么增加一个点或者删除一个点，只需要重新计算 \(g\) ，和所有的新的 \(t\) 即可，对于已有的值只需要乘上或者除以变动的差值，对于新的值也是直接求解即可。
TODO：给出代码实现并通过版本3能通过的题目，因为其理论复杂度接近。
验证
https://www.luogu.com.cn/problem/P4781 （使用版本3）
https://codeforces.com/contest/622/problem/F （使用版本1或者版本2）
