引子
拉格朗日插值可用于求一个用 \(n+1\) 个点确定了表达式的 \(n\) 次多项式某一个点的纵坐标。
算法本身
例题
\(\rm{luogu}\) 有模版。
现给定 \(n+1\) 个点,其中第 \(i\) 个点用 \((x_i, y_i)\) 表示。这 \(n+1\) 个点显然确定了一个 \(n\) 次多项式 \(f(x)\) ,现在给定 \(k\) ,求 \(f(k)\) 。答案对 \(998244353\) 取模。
方案 1:高斯消元
既然有了 \(n+1\) 个点,那么就有了 \(n+1\) 个 \(n+1\) 元方程,可以用高斯消元把所有系数都解出来。时间复杂度是 \(O(n^3)\) 的,显然超时。
方案 2: 拉格朗日插值法
这里先直接写出拉格朗日插值的计算式:
有两种证明方式,第一种是基于中国剩余定理的,第二种是基于构造的,这里讲第二种,因为比较通用。
设 \(P_i = (x_i, y_i), Q_i=(x_i, 0)\) 也就是 \(Q_i\) 是 \(P\) 在 \(x\) 轴上的投影,接下来构建 \(n+1\) 个函数,其中第 \(i\) 个函数 \(g_i(x)\) 经过 \(P_i\) 以及 \(Q_1, Q_2, \dots, Q_{i-1}, Q_{i+1}, \dots, Q_{n+1}\) 。这样的话只要所有 \(g_i(x)\) 合在一起,就是 \(f(x)\) 了。
接下来直接构造 \(g_i(x)\):
设一个常数 \(a\),令 \(g_i(x)=a\prod_{i\ne j} x-x_j\) ,显然对于所有 \(x_j\) ,这个函数值都为 \(0\) ,接下来代入 \(P_i\),即 \(y_i = a\prod_{i\ne j} x_i-x_j\) ,因此有:
因此有:
构造成功,那么就能得到 \(f(x)\) 了:
因此可以用 $O(n^2) $ 的时间算出一个位置的值。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int maxn = 2e3 + 5;
const LL MOD = 998244353;
LL qpow(LL a, LL b) {
LL ret = 1;
while (b) {
if (b & 1) ret = ret * a % MOD;
a = a * a % MOD, b >>= 1;
}
return ret;
}
inline LL inv(LL a) { return qpow(a, MOD - 2); }
struct Point {
LL x, y;
} p[maxn];
int n;
int main() {
// freopen("test.in", "r", stdin);
// freopen("test.out", "w", stdout);
LL px, ans = 0; scanf("%d%lld", &n, &px);
for (int i = 1; i <= n; i++) scanf("%lld%lld", &p[i].x, &p[i].y);
for (int i = 1; i <= n; i++) {
LL s1 = p[i].y, s2 = 1;
for (int j = 1; j <= n; j++) if (j ^ i)
s1 = s1 * (px - p[j].x + MOD) % MOD, s2 = s2 * (p[i].x - p[j].x + MOD) % MOD;
ans = (ans + s1 * inv(s2) % MOD) % MOD;
}
printf("%lld\n", ans);
return 0;
}
简单应用
计算自然数幂和
形式化的,就是给定两个非负整数 \(n, k\),求下面的式子:
其中 \(k \leq 10^6, n\leq 10^9\) 。
这里首先要证明 \(f(n,k)\) 是一个自变量为 \(n\) ,次数为 \(k+1\) 的多项式。
考虑归纳证明:
当 \(k=0\) ,显然有 \(f(n,0)=n\) ,为 \(1\) 次多项式。
现在证明当 \(k>0\) 的时候,
因为有:
接下来构造一个和:
这个东西明显提出一个 \(f(n,k)\) ,而上面这个和明显等于 \((n+1)^{k+1}-1\) ,因此有:
首先,\((n+1)^{k+1}\) 是一个关于 \(n\) 的 \(k+1\) 次多项式,而求和的部分根据归纳,可以证明出也是一个关于 \(n\) 的多项式,其次数一定比 \(k+1\) 小,因此可以得到 \(f(n,k)\) 是一个关于 \(n\) 的 \(k+1\) 次多项式。
当然,实际上有了上面这个式子,你已经可以在 \(O(k^2)\) 的时间内算出一个 \(f(n,k)\) 了,但是这还不够。
因为确定了次数,和自变量,我们直接构造 \(k+2\) 个点,其中第 \(i\) 个点 \(P_i = (i, f(i,k))\) 。接着优化拉格朗日插值的式子:
分式中的四块式子显然是可以用阶乘,前缀和和后缀和预处理的,预处理时间复杂度为 \(O(k)\) ,而计算的时间复杂度也变成了 \(O(k)\),因此总的时间复杂度是 \(O(k)\) 的。
代码
using namespace std;
typedef long long LL;
const int maxk = 1e6 + 5;
const LL MOD = 1e9 + 7;
LL qpow(LL a, LL b) {
LL res = 1;
for (; b; a = a * a % MOD, b >>= 1) if (b & 1) res = res * a % MOD;
return res;
}
inline LL inv(LL a) { return qpow(a, MOD - 2); }
LL n, k, ifac[maxk], pre[maxk], suf[maxk];
int main() {
scanf("%lld%lld", &n, &k);
pre[0] = suf[k + 3] = ifac[0] = 1; LL fac = 1;
for (int i = 1; i <= k + 2; i++) pre[i] = pre[i - 1] * (n - i) % MOD, fac = fac * i % MOD;
for (int i = k + 2; i >= 1; i--) suf[i] = suf[i + 1] * (n - i) % MOD;
ifac[k + 2] = inv(fac);
for (int i = k + 1; i >= 1; i--) ifac[i] = ifac[i + 1] * (i + 1) % MOD;
LL ans = 0, sk = 0;
for (int i = 1; i <= k + 2; i++) {
sk = (sk + qpow(i, k)) % MOD;
LL a = sk * pre[i - 1] % MOD * suf[i + 1] % MOD, b = ifac[i - 1] * ifac[k + 2 - i] % MOD;
if ((k + 2 - i) & 1) ans = (ans - a * b % MOD + MOD) % MOD;
else ans = (ans + a * b % MOD) % MOD;
}
printf("%lld\n", ans);
return 0;
}