拉格朗日插值

这是一个很经典的问题，给定 (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 e i)) 上取值为零，不难发现这样样一个构造（其中 (n) 为需要求点值的位置）：

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

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

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

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

[sumlimits_{i = 1} ^ {k + 1} y_i prod_{j e 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;
}

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

[sumlimits_{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} = sumlimits_{i = 1} ^ n i ^ k)，则可以发现 (S_{n + 1, k} = S_{n, k} + (n + 1) ^ k)。因为我们要证明的是答案是一个关于 (n) 的 (k + 1) 次多项式，因此我们需要往 (k) 上归纳证明。

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

[egin{aligned} S_{n + 1, k + 1} &= S_{n, k + 1} + (n + 1) ^ {k + 1} \ &= S_{n, k + 1} + sumlimits_{i = 0} ^ {k + 1} dbinom{k + 1}{i} n ^ i \ &= sumlimits_{i = 0} ^ {k + 1} dbinom{k + 1}{i} sumlimits_{j = 1} ^ n j ^ i + 1\ &= sumlimits_{i = 0} ^ {k + 1} dbinom{k + 1}{i} S_{n, i} + 1 \ &= sumlimits_{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 = sumlimits_{i = 0} ^ k dbinom{k + 1}{i} S_{n, i} ]

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

[(n + 1) ^ {k + 1} - 1 = sumlimits_{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} - sumlimits_{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) 上的点值，那么可以得到最终的答案：

[sumlimits_{i = 1} ^ {k + 2} y_i prod_{j e i} frac{n - j}{i - j} ]

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