拉格朗日插值

这是一个很经典的问题，给定 \(k + 1\) 个点值，如何快速确定这个 \(k\) 次多项式？

不难发现可以使用待定系数法然后使用高斯消元即可做到 \(O(n ^ 3)\)。但是有些时候我们的目的并非一定要计算出这个多项式的系数，而仅仅想知道这个多项式在某个位置的点值，那么有没有什么直接通过这给定的 \(k + 1\) 个点表示（构造）出这个 \(k\) 次多项式某个位置上的点值的方法呢？拉格朗日插值就解决了这样一个问题。

可以注意到给定的点值 \((x_i, y_i)\) 的实质其实是告诉你该多项式在 \(x = x_i\) 时取值为 \(y_i\)。那么继续延续上面提到的那个想法，如果我们对这 \(k + 1\) 个点值做一些乘法操作，你会发现很难构造出这样一个表示方法，那么可以尝试一下能否使用加法来表示出这样一个多项式。换句话说，我们需要一些式子（有已知 \(k + 1\) 个点表示）的和来构造出这样一个多项式，即本段开头所说的实质。

不难发现，最简单的方法就是让某个式子在 \(x = x_i\) 时的取值为 \(y_i\) 其他式子的取值为 \(0\) 是最简单的一种方法。那么一个方向就逐渐浮现出来了，我们构造 \(k + 1\) 个由已知点值构成的式子其中其中任意一个式子满足恰好在某个 \(x_i\) 上取值 \(y_i\) 其余位置上取值均为 \(0\)。稍加思考可以发现，对于第 \(i\) 个式子，我们想让其在 \(x_j(j \ne i)\) 上取值为零，不难发现这样样一个构造（其中 \(n\) 为需要求点值的位置）：

\[\prod_{j \ne i} (n - x_j) \]

那么要在 \(x_i\) 上取值为 \(y_i\) 怎么办呢？肯定需要往前面乘一个 \(y_i\) 保证之前的性质成立，但于此同时你发现 \(n = x_i\) 时会多算 \(\prod_{j \ne i} (x_i - x_j)\)，因此还需要除去这个数。那么我们可以得到第 \(i\) 个式子的形式化表达：

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

根据前面的理论，将这 \(k + 1\) 个多项式加起来即为最终所得的多项式：

\[\sum\limits_{i = 1} ^ {k + 1} y_i \prod_{j \ne i} \frac{n - x_j}{x_i - x_j} \]

可见，拉格朗日插值的复杂度为 \(O(k ^ 2)\) 还是非常优秀的。当然，上述的拉格朗日插值只是最基础的形式。实际上拉格朗日插值有很多优化，遇到具体问题可以具体分析。

最基础的拉格朗日插值代码如下：

#include <bits/stdc++.h>
using namespace std;
#define rep(i, l, r) for (int i = l; i <= r; ++i)
const int N = 2000 + 5;
const int Mod = 998244353;
int n, k, ans, x[N], y[N];
int read() {
    char c; int x = 0, f = 1;
    c = getchar();
    while (c > '9' || c < '0') { if(c == '-') f = -1; c = getchar();}
    while (c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
int Inc(int a, int b) { return (a += b) >= Mod ? a - Mod : a;}
int Dec(int a, int b) { return (a -= b) < 0 ? a + Mod : a;}
int Mul(int a, int b) { return 1ll * a * b % Mod;}
int fpow(int a, int b) { int ans = 1; for (; b; a = Mul(a, a), b >>= 1) if(b & 1) ans = Mul(ans, a); return ans;}
int main() {
    n = read(), k = read();
    rep(i, 1, n) x[i] = read(), y[i] = read();
    rep(i, 1, n) {
        int tmp = y[i], d = 1;
        rep(j, 1, n) if(i != j) tmp = Mul(tmp, Dec(k, x[j])), d = Mul(d, Dec(x[i], x[j]));
        ans = Inc(ans, Mul(tmp, fpow(d, Mod - 2)));
    }
    printf("%d", ans);
    return 0;
}

下面来看一个拉格朗日插值最经典的应用，求：

\[\sum\limits_{i = 1} ^ n i ^ k (n \le 10 ^ 9, k \le 10 ^ 3) \]

可以发现，\(k = 1\) 时答案为 \(\dfrac{n(n + 1)}{2}\)；\(k = 2\) 时答案为 \(\dfrac{n(n + 1)(2n + 1)}{6}\)；\(k = 3\) 时答案为 \(\dfrac{n ^ 2 (n + 1) ^ 2}{4}\)。那么特别地，我们可以发现答案应该是一个关于 \(n\) 的 \(k + 1\) 次多项式。那么是不是呢，下面我们来证明这个猜想。

直接证明是不好证明的，但这种找出的规律一般可以使用数学归纳法来证明。可以令 \(S_{n, k} = \sum\limits_{i = 1} ^ n i ^ k\)，则可以发现 \(S_{n + 1, k} = S_{n, k} + (n + 1) ^ k\)。因为我们要证明的是答案是一个关于 \(n\) 的 \(k + 1\) 次多项式，因此我们需要往 \(k\) 上归纳证明。

继续由上面的递推式可以得到如下推导：

\[\begin{aligned} S_{n + 1, k + 1} &= S_{n, k + 1} + (n + 1) ^ {k + 1} \\ &= S_{n, k + 1} + \sum\limits_{i = 0} ^ {k + 1} \dbinom{k + 1}{i} n ^ i \\ &= \sum\limits_{i = 0} ^ {k + 1} \dbinom{k + 1}{i} \sum\limits_{j = 1} ^ n j ^ i + 1\\ &= \sum\limits_{i = 0} ^ {k + 1} \dbinom{k + 1}{i} S_{n, i} + 1 \\ &= \sum\limits_{i = 0} ^ k \dbinom{k + 1}{i} S_{n, i} + S_{n, k + 1} + 1 \\ \end{aligned} \]

将第一条式子和最后一条式子左右同时减去 \(S_{n, k + 1}\) 有：

\[(n + 1) ^ {k + 1} - 1 = \sum\limits_{i = 0} ^ k \dbinom{k + 1}{i} S_{n, i} \]

然后将右边取出 \(S_{n, k}\)，有：

\[(n + 1) ^ {k + 1} - 1 = \sum\limits_{i = 0} ^ {k - 1} \dbinom{k + 1}{i} S_{n, i} + (k + 1)S_{n, k} \]

移项可得：

\[S_{n, k} = \frac{1}{k + 1} ((n + 1) ^ {k + 1} - \sum\limits_{i = 0} ^ {k - 1} S_{n, i} - 1) \]

那么如果 \(S_{n, i} (i < k)\)，满足 \(S_{n, i}\) 是关于 \(n\) 的 \(i + 1\) 次多项式，那么就可以归纳证得 \(S_{n, k}\) 为 \(k + 1\) 次多项式。

那么直接使用拉格朗日插值即可做到 \(O(k ^ 2)\) 的复杂度。但事实上因为这里的点值是你自取的，那么最简单地我们取 \(1, 2, \cdots, k + 2\) 上的点值，那么可以得到最终的答案：

\[\sum\limits_{i = 1} ^ {k + 2} y_i \prod_{j \ne i} \frac{n - j}{i - j} \]

观察后不难得到：后面的连乘部分的分母部分实际上是 \((i - 1)! (n - i)! \times (-1) ^ {n - i}\)，直接维护阶乘即可；而对于分母，你会发现对于所有的 \(i\) 不同的唯一之处是缺少了 \(n - i\) 这个部分，于是我们可以预处理 \(n - i\) 的前缀积后缀积即可。这样我们就在 \(O(k)\) 的时间复杂度内解决了本题。