拉格朗日插值入门

引子

拉格朗日插值可用于求一个用 \(n+1\) 个点确定了表达式的 \(n\) 次多项式某一个点的纵坐标。

算法本身

例题

\(\rm{luogu}\) 有模版。

现给定 \(n+1\) 个点，其中第 \(i\) 个点用 \((x_i, y_i)\) 表示。这 \(n+1\) 个点显然确定了一个 \(n\) 次多项式 \(f(x)\) ，现在给定 \(k\) ，求 \(f(k)\) 。答案对 \(998244353\) 取模。

\[\begin{aligned} &n\leq 2\times 10^3\\ &\forall i\in \mathbb{N}^+ \cap [1, n+1], x_i, y_i < 998244353\\ & k < 998244353\\ &\forall i,j\in \mathbb{N}^+ \cap [1, n+1] \land i\ne j, x_i \ne x_j \end{aligned} \]

方案 1：高斯消元

既然有了 \(n+1\) 个点，那么就有了 \(n+1\) 个 \(n+1\) 元方程，可以用高斯消元把所有系数都解出来。时间复杂度是 \(O(n^3)\) 的，显然超时。

方案 2: 拉格朗日插值法

这里先直接写出拉格朗日插值的计算式：

\[f(x)=\sum_{i=1}^{n+1} y_i\prod_{j\ne i} \frac{x-x_i}{x_j-x_i} \]

有两种证明方式，第一种是基于中国剩余定理的，第二种是基于构造的，这里讲第二种，因为比较通用。

设 \(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\) ，因此有：

\[a=\frac{y_i}{\prod_{j\ne i} x_i - x_j} \]

因此有：

\[g_i(x)=y_i \prod_{j\ne i} \frac{x-x_j}{x_i-x_j} \]

构造成功，那么就能得到 \(f(x)\) 了：

\[f(x)=\sum_{i=1}^{n+1}g_i(x)=\sum_{i=1}^{n+1} y_i\prod_{j\ne i} \frac{x-x_j}{x_i-x_j} \]

因此可以用 $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;
}

简单应用

计算自然数幂和

\(\Rightarrow \rm luogu\) 链接

形式化的，就是给定两个非负整数 \(n, k\)，求下面的式子：

\[f(n,k) = \sum_{i=1}^n i^k \]

其中 \(k \leq 10^6, n\leq 10^9\) 。

这里首先要证明 \(f(n,k)\) 是一个自变量为 \(n\) ，次数为 \(k+1\) 的多项式。

考虑归纳证明：

当 \(k=0\) ，显然有 \(f(n,0)=n\) ，为 \(1\) 次多项式。
现在证明当 \(k>0\) 的时候，

因为有：

\[\begin{aligned} (i+1)^{k+1} - i^{k+1} &= \left[\sum_{x=0}^{k+1} \binom{k+1}{x} i^x\right]-i^{k+1} \\ &= \sum_{x=0}^k \binom{k+1}{x} i^x \end{aligned} \]

接下来构造一个和：

\[\begin{aligned} \sum_{i=1}^n (i+1)^{k+1} - i^{k+1} &= \sum_{i=1}^n \sum_{x=0}^k \binom{k+1}{x} i^x\\ &= \sum_{x=0}^k \binom{k+1}{x} \sum_{i=1}^n i^x\\ &= \sum_{x=0}^k \binom{k+1}{x} f(n, x) \end{aligned} \]

这个东西明显提出一个 \(f(n,k)\) ，而上面这个和明显等于 \((n+1)^{k+1}-1\) ，因此有：

\[\begin{aligned} (n+1)^{k+1}-1 &= \binom{k+1}{k} f(n, k) +\sum_{x=0}^{k-1} \binom{k+1}{x} f(n,x)\\ (n+1)^{k+1}-1 &= (k+1)f(n,k) + \sum_{x=0}^{k-1} \binom{k+1}{x} f(n,x)\\ f(n,k) &= \frac{1}{k+1}\left[(n+1)^{k+1}-1 - \sum_{x=0}^{k-1} \binom{k+1}{x} f(n,x)\right] \end{aligned} \]

首先，\((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))\) 。接着优化拉格朗日插值的式子：

\[\begin{aligned} f(n,k) &= \sum_{i=1}^{k+2} f(i,k) \prod_{j \ne i} \frac{n-j}{i-j}\\ &= \sum_{i=1}^{k+2} (-1)^{k+2-i} f(i,k) \frac {\left(\prod_{j=1}^{i-1} n-i\right) \left(\prod_{j=i+1}^{k+2} n-i\right)} {n!(k+2-n)!} \end{aligned} \]

分式中的四块式子显然是可以用阶乘，前缀和和后缀和预处理的，预处理时间复杂度为 \(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;
}