【算法竞赛 - 数学】拉格朗日插值法 · purinliang/MyCareer Wiki
模板
版本1： \(x\) 取从 \(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;
版本2： \(x\) 取从 \(1\) 开始的连续正整数，但是不预处理。
时间复杂度 \(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;
版本3：对于给定的取值 \(x[]\) 进行求解。
时间复杂度 \(O(k^2)\) ，空间复杂度 \(O(1)\) 。
使用：
1. 直接调用 int lagrange (int x0, int *x, int *y, int k) 即可，并且可以多次调用。
namespace Lagrange3 {
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 *x, int *y, int k) {
if (x0 <= k)
return y[x0];
ll Lx0 = 0, P = 1;
for (int i = 1; i <= k; ++i)
P = P * (x0 - x[i] + MOD) % MOD;
for (int i = 1; i <= k; ++i) {
ll p = P * inv (x0 - x[i] + MOD) % MOD, q = 1;
for (int j = 1; j <= k; ++j) {
if (j == i) continue;
q = q * (x[i] - x[j] + MOD) % MOD;
}
Lx0 = (Lx0 + p * inv (q) % MOD * y[i]) % MOD;
}
return Lx0;
}
};
using Lagrange3::lagrange;
原理
对某个多项式函数，已知有
在算法竞赛中，通常要求在模 \(MOD\) 的意义下进行，且 \(MOD\) 是某个较大的质数。
TODO
细节
拉格朗日插值法的 \(k\) 的取值很容易让人迷惑，本文中所有的公式、代码都用 \(k-1\) 表示多项式的最高次数，也就是不超过 \(k-1\) 次的多项式。
- 对一个不超过 \(k-1\) 次的多项式插值，至少需要 \(k\) 个点。
- 一个不超过 \(k-2\) 次多项式的部分和，通常是一个不超过 \(k-1\) 次的多项式。
上面的式子中，给出某个 \(x\) ，然后 \(O(k^2+k\log M)\) 计算出 \(f(x)\) 的值，慢得离谱，慎用。
拓展
选取等差数列进行插值
例如：求 \(\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）
引申阅读
重心拉格朗日插值法、牛顿插值法、差分法、高斯消元法、龙格现象