本节点所需要的前置研发项目:多项式基础

引入

问题

先给定 $n$ 个点以及一个数 $k$,有一个函数经过了这 $n$ 个点,请求出 $k$ 在函数上对应的函数值。

如何解决这个问题呢?

下面给出一种较为暴力的方法。

我们设这个函数的解析式为 $f(x)=\sum\limits_{i=0}^{n-1}a_ix^i$,将点值代入,直接高斯消元解出系数,然后把 $k$ 带入式子中去即可。时间复杂度为 $O(n^3)$。

然而有一种神奇的方法可以在 $O(n^2)$ 的时间内解决此问题。这种方法就是拉格朗日插值法

公式

数值分析中,拉格朗日插值法是以法国18世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示各结果之间某种内在联系或规律,而不少函数都只能通过繁复实验和多次观测来了解。而,如果对实践中的某个物理量进行观测,在若干个不同的地方得到相应的观测值,拉格朗日插值法可以找到一个多项式,其恰好在各个观测的点取到观测到的值。上面这样的多项式就称为拉格朗日(插值)多项式(Lagrange polynomial)。

——摘自 维基百科-拉格朗日插值法

首先放上公式。
$$
f(k)=\sum\limits_{i=1}^ny_i\prod_{i\ne j}\dfrac{k-x_j}{x_i-x_j}
$$
其中的 $x_i,y_j$ 都是给出的点值。

为什么这个式子就可以保证 $k$ 在 $n$ 个点确定的函数上呢?我们带入一个 $x_p$ 试试。($p$ 是给出的点值)

这里有一个 Sigma,我们分类讨论。

如果 $i=p$,那么 Prod 后面的分式分子分母都相同,取值为 1,这一部分的值就是 $y_p$

如果 $i\ne p$,那么总有一个 $j$ 满足 $j=p$ 使得 Prod 后面的分式分子为 0,这一部分取值为 0.

综上所述,如果带入 $x_p$,得到的函数值为 $y_p$,这就保证了 $k$ 在 $n$ 个点确定的函数上。

代码实现

下面给出 C++ 实现的代码。代码中的函数值对 998244353 取模,所以使用了乘法逆元。代码中的qpow(fac, P - 2) 在取模 P 的意义下等价于 $fac^{-1}$(费马小定理)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
// Author: SocialZxy
#include<bits/stdc++.h>
#define int long long
using namespace std;

int get() {
int x = 0, f = 1; char c = getchar();
while(!isdigit(c)) { if(c == '-') f = -1; c = getchar(); }
while(isdigit(c)) { x = x * 10 + c - '0'; c = getchar(); }
return x * f;
}

const int N = 2e3 + 5, P = 998244353;
int n, k, x[N], y[N];

int qpow(int x, int y) { // 快速幂
int res = 1;
while(y) {
if(y & 1) res = (res * x) % P;
x = x * x % P;
y >>= 1;
}
return res;
}

signed main() {
n = get(), k = get();
for(int i = 1; i <= n; i++) x[i] = get(), y[i] = get(); // 读入点值
int ans = 0;
for(int i = 1; i <= n; i++) { // 拉格朗日插值
int res = 1, fac = 1;
for(int j = 1; j <= n; j++) {
if(j == i) continue;
res = res * (k - x[j]) % P;
fac = fac * (x[i] - x[j]) % P;
}
res = y[i] * res % P * qpow(fac, P - 2) % P;
ans = (ans + res) % P;
}
cout << (ans + P) % P;
return 0;
}