这是一个很经典的问题,给定 (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)) 的时间复杂度内解决了本题。

GO!