给你 n+1n+1 个点，求穿过它们的 nn 次多项式。于是你想到了待定系数高斯消元，不过这 O(n3)O(n^3) 的太烂了。于是我们引入拉格朗日插值法，通过构造法整出多项式。

拉格朗日插值法

记这 n+1n+1 个点为

(x0,y0),(x1,y1),(x2,y2),…,(xn,yn)(x_0,y_0),(x_1,y_1),(x_2,y_2),\dots,(x_n,y_n)

定义拉格朗日基本多项式 ℓi(x)\ell_i(x)

ℓj(x)=∏i≠jx−xixj−xi\ell_j(x)=\prod_{i\ne j}\frac{x-x_i}{x_j-x_i}

很容易发现一个性质：

∀i≠j:ℓj(xi)=0\forall i\ne j :\ell_j(x_i)=0

ℓj(xj)=1\ell_j(x_j) =1

于是我们得到拉格朗日插值多项式为：

L(x)=∑j=0nyjℓj(x)=∑j=0nyj∏i≠jx−xixj−xi\begin{aligned} L(x)&=\sum_{j=0}^{n} y_j\ell_j(x)\\ &=\sum_{j=0}^n y_j\prod_{i\ne j}\frac{x-x_i}{x_j-x_i} \end{aligned}

这显然满足 L(xi)=yiL(x_i)=y_i 。

如果要求出每一项系数，需要多项式乘法，复杂度 O(n2log⁡n)O(n^2\log n) ，若只是求 L(k)L(k) 复杂度 O(n2)O(n^2) 。若希望多点求值，见 重心拉格朗日插值法 。

实现

【模板】拉格朗日插值

求 L(k)L(k)

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
// 拉格朗日插值
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 2005;
const int mod = 998244353;
int n, k;
ll x[N], y[N];
ll fpow(ll a, ll b, ll p) {
    ll r = 1;
    for ( ; b; b >>= 1, a = (a * a) % p) if (b & 1) r = (r * a) % p;
    return r;
}
ll lagrange(int n, ll k, ll x[], ll y[]) {
    ll r = 0;
    for (int j = 1; j <= n; ++j) {
        ll P1 = 1, P2 = 1;
        for (int i = 1; i <= n; ++i) if (i != j) {
            P1 = (P1 * (k - x[i] + mod)) % mod;
            P2 = (P2 * (x[j] - x[i] + mod)) % mod;
        }
        r = (r + y[j] * P1 % mod * fpow(P2, mod - 2, mod) % mod) % mod;
    }
    return r;
}
int main() {
    scanf("%d%d", &n, &k);
    for (int i = 1; i <= n; ++i) scanf("%lld%lld", &x[i], &y[i]);
    printf("%lld\n", lagrange(n, k, x, y));
    return 0;
}

当 x 的值连续

当 xix_i 的取值连续时，可以把复杂度优化到 O(n)O(n) 。令 xi=ix_i=i ，那么：

L(k)=∑j=0nyi∏i≠jk−ij−iL(k)=\sum_{j=0}^n y_i\prod_{i\ne j}\frac{k-i}{j-i}

考虑快速求 ∏i≠jk−ij−i\prod_{i\ne j} \frac{k-i}{j-i} 。

对于分子预处理前缀积、后缀积：

pi=∏j=0i−1k−jsi=∏j=i+1nk−jp_i=\prod_{j=0}^{i-1} k-j\quad s_i=\prod_{j=i+1}^n k-j

分母就是阶乘，同样预处理，但需要注意分母的符号， n−in-i 为奇数时分母取负。

L(k)=∑j=0n(−1)n−jyjpjsjj!(n−j)!L(k)=\sum_{j=0}^n (-1)^{n-j}y_j \frac{p_{j}s_{j}}{j!(n-j)!}

预处理阶乘逆元，时间复杂度的瓶颈就不会在求逆元上。

实现

CF622F，求 ∑i=1nik(n≤109,k≤106)\sum_{i=1}^n i^k\quad(n\le 10^9, k\le 10^6) ，我们熟知 ∑i=1nik=f(n)\sum_{i=1}^n i^k=f(n) 其中 f(n)f(n) 是一个 k+1k+1 次的多项式，那么我们只要暴力算出 n=0,1,2,3,…,k+1n=0,1,2,3,\dots,k+1 的值，再用拉格朗日插值算 f(n)f(n) 即可。

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
// K 次方和
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1000010;
const int mod = 1e9 + 7;
int n, k;
ll fpow(ll a, ll b, ll p = mod) {
    ll r = 1;
    for ( ; b; b >>= 1, a = (a * a) % p) if (b & 1) r = (r * a) % p;
    return r;
}
ll a[N], s[N], p[N], fac[N], inv[N], ifac[N];
void init(int n) {
    fac[0] = fac[1] = inv[0] = inv[1] = ifac[0] = ifac[1] = 1;
    for (int i = 2; i <= n; ++i) fac[i] = fac[i - 1] * i % mod;
    for (int i = 2; i <= n; ++i) inv[i] = (mod - mod / i) * inv[mod % i] % mod;
    for (int i = 2; i <= n; ++i) ifac[i] = ifac[i - 1] * inv[i] % mod;
}
ll lagrange(int n, ll k, ll y[]) {
    init(n);
    // x_i = i (i in [0, n])
    p[0] = 1;
    for (int i = 0; i < n; ++i) p[i + 1] = p[i] * (k - i + mod) % mod;
    s[n] = 1;
    for (int i = n; i >= 1; --i) s[i - 1] = s[i] * (k - i + mod) % mod;
    ll r = 0;
    for (int j = 0; j <= n; ++j) {
        ll tmp = p[j] * s[j] % mod * ifac[j] % mod * ifac[n - j] % mod;
        r = (r + (ll)(((n - j) & 1) ? (mod - 1) : 1) * y[j] % mod * tmp % mod) % mod;
    }
    return r;
}
int main() {
    cin >> n >> k;
    for (int i = 1; i <= k + 1; ++i) a[i] = (a[i - 1] + fpow(i, k)) % mod;
    cout << lagrange(k + 1, n, a) << endl;
    return 0;
}

重心拉格朗日插值法

若每次多加一个点重新求插值，每次重求 O(n2)O(n^2) 的复杂度太烂了。于是我们改进一下：

令：

ℓ(x)=∏i=0n(x−xi)\ell(x)=\prod_{i=0}^n (x-x_i)

那么拉格朗日基本多项式重新写为：

ℓj(x)=ℓ(x)x−xj1∏i≠j(xj−xi)\ell_j(x)=\frac{\ell(x)}{x-x_j} \frac{1}{\prod_{i\ne j} (x_j-x_i)}

定义重心权：

wj=1∏i≠j(xj−xi)w_j=\frac{1}{\prod_{i\ne j}(x_j-x_i)}

则：

L(x)=ℓ(x)∑j=0nwjx−xjyjL(x)=\ell(x)\sum_{j=0}^n \frac{w_j}{x-x_j}y_j

一开始算出所有 wjw_j 复杂度 O(n2)O(n^2) ，之后插一个点只要修改 wjw_j 即可，复杂度 O(n)O(n) 。计算一个 L(k)L(k) 是 O(n)O(n) （注意瓶颈在求逆，对 nn 个数求逆可以做到 O(n+log⁡p)O(n+\log p) 见：数论基础 - 线性求任意-n-个数逆元 。