在数值分析中，拉格朗日插值法是以法国18世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。如果对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。上面这样的多项式就称为拉格朗日（插值）多项式。

拉格朗日插值

考虑这样一个多项式： $f(x) = f_0 + f_1 x + f_2 x^2 + \dots + f_n x^n$ 。

假设我们已知 $n+1$ 个互不相同的点 $(x_0,x_1,\dots,x_n)$ ，该多项式在这些点的取值为 $(y_0,y_1,\dots,y_n)$ ，我们可以用矩阵描述这 $n+1$ 个等式：

\left( \begin{array}{ccccc} \tag{1} 1 & x_0 & x_0^2 & \cdots & x_0^n \\ 1 & x_1 & x_1^2 & \cdots & x_1^n \\ 1 & x_2 & x_2^2 & \cdots & x_2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \\ \end{array} \right)\left( \begin{array}{c} f_0 \\ f_1 \\ f_2 \\ \vdots \\ f_n \\ \end{array} \right)= \left( \begin{array}{c} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{array} \right)

可以发现系数矩阵是一个范德蒙特矩阵，而对于范德蒙特行列式有如下结论：

\left \vert \begin{array}{ccccc} \tag{2} 1 & x_0 & x_0^2 & \cdots & x_0^n \\ 1 & x_1 & x_1^2 & \cdots & x_1^n \\ 1 & x_2 & x_2^2 & \cdots & x_2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{array} \right \vert = \prod_{j \lt i} (x_i - x_j) \neq 0

结合线性代数知识可知，给定 $n+1$ 个 $x$ 不同的点，可以唯一确定经过这些点的次数不超过 $n$ 的多项式。

由 $(1)$ 式，我们只需要求出范德蒙特矩阵的逆，再在等式两边左乘就可以得到系数向量 $(f_0,f_1,\dots,f_n)$ ，但这样计算量很大。因此在确定多项式时，一般不选择直接求系数，而是用一组基本多项式来表示该多项式。

引入拉格朗日基本多项式：

\begin{aligned} \tag{3} L_i(x) &= \prod_{j \neq i} \frac{x - x_j}{x_i - x_j} \\ &= \frac{(x - x_0) (x - x_1)\dots(x-x_{i-1})(x-x_{i+1})\dots (x-x_n)}{(x_i - x_0) (x_i - x_1)\dots(x_i-x_{i-1})(x_i-x_{i+1})\dots (x_i-x_n)} \end{aligned}

从而 $f(x)$ 可以表示为：

\begin{aligned} \tag{4} y_0 L_0(x) + y_1 L_1(x) + \dots + y_n L_n(x) = \sum \limits_{i=0}^{n}y_iL_i(x) \end{aligned}

可以证明该表示下 $f(x)$ 是经过 $(x0,y0),(x_1,y_1),\dots,(x_n,y_n)$ 的。我们要求 $f(N)$ ，只需求得各基本多项式在 $N$ 处的值，再按权相加即可。

上述过程代码实现如下，时间复杂度为 $O(n^2)$ 。

拉格朗日插值

// 由于结果可能很大，一般要对答案取模，这里模数为 MOD
// qpow 为快速幂，用于快速求逆元
// 参数依次为插值点数，插值点，待求值点横坐标，返回值为待求值点纵坐标
LL lagrange(int n, LL* x, LL* y, LL k) {
    LL ans = 0;
    for (int i = 1; i <= n; i++) {
        LL num = 1, d = 1;
        for (int j = 1; j <= n; j++)
            if (j != i)
                num = (k - x[j]) % MOD * num % MOD,
                d = (x[i] - x[j]) % MOD * d % MOD;
        ans =
            (ans + y[i] % MOD * num % MOD * qpow(d, MOD - 2, MOD) % MOD) % MOD;
    }
    return (ans + MOD) % MOD;
}

求插值多项式系数

由 $(4)$ 式可以确定系数向量 $(f_0,f_1,\dots,f_n)$ 。

由于 $y_iL_i(x) = \frac{y_i}{\prod \limits_{j \neq i} {x_i-x_j}} \prod \limits_{j\neq i} {x-x_j} = C_i*\prod \limits_{j \neq i} {x-x_j} =C_i* \frac{\prod \limits_{0 \le j \le n} {x-x_j}}{x-x_i}$ ，令 $p_i(x) = \prod \limits_{j \neq i} {x-x_j},q(x)=\prod \limits_{0 \le j \le n} {x-x_j}$ ，则 $f(x) = \sum \limits_{i=0}^{n}y_i*L_i(x) = \sum \limits_{i=0}^{n}C_i*p_i(x)$ 。模拟多项式乘法可以得到 $q(x)$ 的系数向量，记为 $(q_0,q_1,\dots,q_n)$ ，同时记 $p_i(x)$ 的系数向量为 $(p_{i,0},p_{i,1},\dots,p_{i,n})$ ，则由 $q(x) = (x-x_i)p_i(x)$ 可得 $p_{i,j}=\frac{p_{i,j-1}-q_j}{x_j}$ 。由此递推式可 $O(n)$ 求得 $p_i(x)$ 的各次项系数，从而可以 $O(n^2)$ 求得 $f(x)$ 的各次项系数。

求插值多项式系数代码实现如下：

插值多项式系数

// 参数为插值点数，插值点，返回值为系数向量
vector<LL> lagrange_coef(int n, LL* x, LL* y) {
    vector<LL> coef(n), tot(n + 1);
    tot[1] = 1, tot[0] = -x[1] % MOD;
    for (int i = 2; i <= n; i++) {
        vector<LL> v(n + 1);
        for (int j = 1; j <= i; j++) v[j] = tot[j - 1];
        for (int j = 0; j < i; j++)
            v[j] = (v[j] - x[i] % MOD * tot[j] % MOD) % MOD;
        tot = v;
    }
    for (int i = 1; i <= n; i++) {
        LL mul = 1, invx = qpow(x[i] % MOD, MOD - 2, MOD);
        for (int j = 1; j <= n; j++)
            if (j != i) mul = (x[i] - x[j]) % MOD * mul % MOD;
        mul = y[i] % MOD * qpow(mul, MOD - 2, MOD) % MOD;
        LL v = 0;
        for (int j = 0; j < n; j++) {
            v = (v - tot[j]) * invx % MOD;
            coef[j] = (coef[j] + mul * v % MOD) % MOD;
        }
    }
    return coef;
}

注：AtCoder Beginner Contest 208 F 的官方题解提到了插值多项式系数可以利用范德蒙特矩阵在 $O(nlog^2n)$ 时间内求得，AC Library 也有相应代码，感兴趣的朋友可以自行阅读。

拉格朗日插值优化

横坐标连续时插值

如果题目不限制，我们取插值点时一般取横坐标依次递增 $1$ 的 $n+1$ 个连续的插值点，这时可对 $(3)$ 式进行变形：

\tag{5} L_i(x) = (-1)^{n-i}\frac{pre[i-1]*suf[i+1]}{fac[i-1]*fac[n-i]}

这里 $fac[i]=i!, pre[i] = (x-x_0)(x-x_1)\dots(x-x_i), suf[i] = (x-x_i)(x-x_{i+1})\dots(x-x_n)$ 。 $O(n)$ 预处理前后缀乘积以及阶乘逆元即可 $O(1)$ 计算 $L_i(x)$ ，从而在 $O(n)$ 时间内完成插值。

代码实现如下：

横坐标连续时插值

// 参数为插值点数，插值点，待求值点横坐标，返回值为待求值点纵坐标
LL lagrange_cont(int n, LL* x, LL* y, LL k) {
    LL ans = 0;
    vector<LL> pre(n + 1), suf(n + 2), fac(n), inv(n);
    pre[0] = fac[0] = inv[0] = suf[n + 1] = 1;
    for (int i = 1; i <= n; i++) pre[i] = (k - x[i]) % MOD * pre[i - 1] % MOD;
    for (int i = n; i >= 1; i--) suf[i] = (k - x[i]) % MOD * suf[i + 1] % MOD;
    for (int i = 1; i < n; i++) fac[i] = fac[i - 1] * i % MOD;
    inv[n - 1] = qpow(fac[n - 1], MOD - 2, MOD);
    for (int i = n - 2; ~i; i--) inv[i] = inv[i + 1] * (i + 1) % MOD;
    for (int i = 1; i <= n; i++)
        ans =
            (ans + ((n - i) % 2 ? -1ll : 1ll) * y[i] % MOD * pre[i - 1] % MOD *
                       suf[i + 1] % MOD * inv[i - 1] % MOD * inv[n - i] % MOD) %
            MOD;
    return (ans + MOD) % MOD;
}

重心拉格朗日插值

重心拉格朗日插值适用于插值点动态变化时，单次插值时间复杂度为 $O(n)$ 或 $O(nlogn)$ ，取决于能否预处理乘法逆元。

由 $(4)$ 式，有：

\tag{6} f(x)= \sum \limits_{i=0}^{n} y_i \prod \limits_{j \neq i} \frac{x-x_j}{x_i-x_j}

令 $g = \prod \limits_{0 \le j \le n} (x-x_j)$ ，有：

\tag{7}f(x)= g\sum \limits_{i=0}^{n} \prod \limits_{j \neq i} \frac{y_i}{(x-x_i)(x_i-x_j)}

再令 $t_i = \prod \limits_{j \neq i} \frac{y_i}{x_i-x_j}$ ，则：

\tag{8} f(x)= g\sum \limits_{i=0}^{n} \frac{t_i}{x-x_i}

对每个点 $i$ ，我们记录它的 $t_i$ 和 $x_i$ 。当加入新的点时，我们只需更新 $t_i$ 。

重心拉格朗日插值

struct CentroidLagrange {
    vector<LL> t, x;
    int size = 0;
    // 插入点 (x,y)
    void insert(LL _x, LL _y) {
        LL tmp = 1;
        for (int i = 0; i < size; i++) {
            tmp = (_x - x[i]) % MOD * tmp % MOD;
            t[i] = t[i] * qpow((x[i] - _x) % MOD, MOD - 2, MOD) % MOD;
        }
        tmp = qpow(tmp, MOD - 2, MOD);
        size++;
        x.push_back(_x);
        t.push_back(_y % MOD * tmp % MOD);
    }
    // 计算 x=k 处的值
    LL calc(LL k) {
        LL ans = 0;
        vector<LL> pre(size), suf(size);
        for (int i = 0; i < size; i++)
            pre[i] = (k - x[i]) % MOD * (i ? pre[i - 1] : 1) % MOD;
        for (int i = size - 1; ~i; i--)
            suf[i] = (k - x[i]) % MOD * (i < size - 1 ? suf[i + 1] : 1) % MOD;
        for (int i = 0; i < size; i++)
            ans = (ans + (i ? pre[i - 1] : 1ll) *
                             (i < size - 1 ? suf[i + 1] : 1ll) % MOD * t[i] %
                             MOD) %
                  MOD;
        return (ans + MOD) % MOD;
    }
};

多项式多点插值与快速插值

待补充…

例题

A 洛谷拉格朗日插值

题目链接

套拉格朗日插值板子即可。

B LOJ 拉格朗日插值

题目链接

套重心拉格朗日插值板子即可。

C The Sum of the k-th Powers

题目链接
$\sum_{i=1}^{n} i^k$ 是 $N$ 的 $k+1$ 次多项式，因此我们取 $k+2$ 个连续的插值点然后 $O(k)$ 插值即可。

证明的话要用到一个定理： $f(x)$ 是 $n$ 次多项式 ⇔ $f(x)$ 的 $n$ 阶差分为常数。注意到 $\sum_{i=1}^{n} i^k$ 的一阶差分是 $k$ 次多项式，因此 $\sum_{i=1}^{n} i^k$ 是 $k+1$ 次多项式。

The Sum of the k-th Powers

#include <bits/stdc++.h>

#define debug(x) cerr << #x << " = " << x << endl

using namespace std;
typedef long long LL;

const int MAXN = 1E6 + 5;
const int MOD = 1e9 + 7;
LL n, m, T;
LL x[MAXN], y[MAXN];

LL qpow(LL a, LL b, LL MOD) {
    LL ans = 1, base = a;
    while (b) {
        if (b & 1) ans = ans * base % MOD;
        b >>= 1, base = base * base % MOD;
    }
    return ans;
}

// 参数为插值点数，插值点，待求值点横坐标，返回值为待求值点纵坐标
LL lagrange_cont(int n, LL* x, LL* y, LL k) {
    LL ans = 0;
    vector<LL> pre(n + 1), suf(n + 2), fac(n), inv(n);
    pre[0] = fac[0] = inv[0] = suf[n + 1] = 1;
    for (int i = 1; i <= n; i++) pre[i] = (k - x[i]) % MOD * pre[i - 1] % MOD;
    for (int i = n; i >= 1; i--) suf[i] = (k - x[i]) % MOD * suf[i + 1] % MOD;
    for (int i = 1; i < n; i++) fac[i] = fac[i - 1] * i % MOD;
    inv[n - 1] = qpow(fac[n - 1], MOD - 2, MOD);
    for (int i = n - 2; ~i; i--) inv[i] = inv[i + 1] * (i + 1) % MOD;
    for (int i = 1; i <= n; i++)
        ans =
            (ans + ((n - i) % 2 ? -1ll : 1ll) * y[i] % MOD * pre[i - 1] % MOD *
                       suf[i + 1] % MOD * inv[i - 1] % MOD * inv[n - i] % MOD) %
            MOD;
    return (ans + MOD) % MOD;
}

signed main() {
    ios::sync_with_stdio(false);
    cin >> n >> m;
    for (int i = 1; i <= m + 2; i++)
        x[i] = i, y[i] = (y[i - 1] + qpow(i, m, MOD)) % MOD;
    cout << lagrange_cont(m + 2, x, y, n);
    return 0;
}

D 教科书般的亵渎

题目链接
假设 $x$ 个怪物的血量为 $1,2,\dots,x$ ，那么使用一张亵渎可以将这 $x$ 个怪全部击杀，获得分数为 $\sum_{i=1}^{x} i^k$ 。其中 $k$ 为需要使用的亵渎张数，易知 $k = m+1$ ，故 $\sum_{i=1}^{x} i^k$ 是 $m+2$ 次多项式，取 $m+3$ 个连续的点插值即可。

下面考虑如果某些血量不存在该如何处理：

设血量为 $a_i$ 的怪不存在（规定 $a_0=0$ 不存在），那么当血量比它少的怪全部被击杀后，其余怪的剩余血量为 $1,2,\dots,n-a_i$ ，由此可计算再次使用一张亵渎的得分。当然，大于 $a_i$ 的某些血量可能不存在，从小到大记为 $a_{i+1},a_{i+2},\dots,a_m$ ，我们需要把相应的得分扣除。由于此次亵渎会停在 $a_{i+1}$ ，接下来只需要从 $a_{i+1}$ 开始重复上述过程。

时间复杂度： $O(m^2)$ 。

教科书般的亵渎

#include <bits/stdc++.h>

#define debug(x) cerr << #x << " = " << x << endl

using namespace std;
typedef long long LL;

const int MAXN = 55;
const int MOD = 1e9 + 7;
LL n, m, T;
LL x[MAXN], y[MAXN], a[MAXN];

LL qpow(LL a, LL b, LL MOD) {
    LL ans = 1, base = a;
    while (b) {
        if (b & 1) ans = ans * base % MOD;
        b >>= 1, base = base * base % MOD;
    }
    return ans;
}

// 参数为插值点数，插值点，待求值点横坐标，返回值为待求值点纵坐标
LL lagrange_cont(int n, LL* x, LL* y, LL k) {
    LL ans = 0;
    vector<LL> pre(n + 1), suf(n + 2), fac(n), inv(n);
    pre[0] = fac[0] = inv[0] = suf[n + 1] = 1;
    for (int i = 1; i <= n; i++) pre[i] = (k - x[i]) % MOD * pre[i - 1] % MOD;
    for (int i = n; i >= 1; i--) suf[i] = (k - x[i]) % MOD * suf[i + 1] % MOD;
    for (int i = 1; i < n; i++) fac[i] = fac[i - 1] * i % MOD;
    inv[n - 1] = qpow(fac[n - 1], MOD - 2, MOD);
    for (int i = n - 2; ~i; i--) inv[i] = inv[i + 1] * (i + 1) % MOD;
    for (int i = 1; i <= n; i++)
        ans =
            (ans + ((n - i) % 2 ? -1ll : 1ll) * y[i] % MOD * pre[i - 1] % MOD *
                       suf[i + 1] % MOD * inv[i - 1] % MOD * inv[n - i] % MOD) %
            MOD;
    return (ans + MOD) % MOD;
}

signed main() {
    ios::sync_with_stdio(false);
    cin >> T;
    while (T--) {
        cin >> n >> m;
        for (int i = 1; i <= m; i++) cin >> a[i];
        a[0] = 0;
        sort(a + 1, a + m + 1);
        LL ans = 0;
        for (int i = 1; i <= m + 3; i++)
            x[i] = i, y[i] = (y[i - 1] + qpow(i, m + 1, MOD)) % MOD;
        for (int i = 0; i <= m; i++) {
            ans = (ans + lagrange_cont(m + 3, x, y, n - a[i])) % MOD;
            for (int j = i + 1; j <= m; j++)
                ans = (ans - qpow(a[j] - a[i], m + 1, MOD)) % MOD;
        }
        cout << (ans + MOD) % MOD << endl;
    }
    return 0;
}

E 排名估算

题目链接
关键是推式子，这里要用到全概率公式和贝叶斯公式。

设 $M$ 表示抽到 $m$ 个人排名比小 C 低， $R_i$ 表示排名为 $i$ ，则：

\begin{aligned} P(R_i|M) &= \frac{P(R_i)P(M|R_i)}{\sum_{j=1}^{n}P(R_j)P(M|R_j)} \\ &= \frac{(n-i)^m}{\sum_{j=1}^{n}(n-j)^m} \end{aligned}

因此答案为：

\sum_{i=1}^{n}i*P(R_i|M)=\frac{\sum_{i=1}^{n}i(n-i)^m}{\sum_{j=1}^{n}(n-j)^m}

换元可证分子是 $m+2$ 次多项式，分母是 $m+1$ 次多项式，使用拉格朗日插值求解即可。

排名估算

#include <bits/stdc++.h>

#define debug(x) cerr << #x << " = " << x << endl

using namespace std;
typedef long long LL;

const int MAXN = 5005;
const int MOD = 998244353;
LL n, m, x[MAXN], y[MAXN];

LL qpow(LL a, LL b, LL MOD) {
    LL ans = 1, base = a;
    while (b) {
        if (b & 1) ans = ans * base % MOD;
        b >>= 1, base = base * base % MOD;
    }
    return ans;
}

LL lagrange_cont(int n, LL* x, LL* y, LL k) {
    LL ans = 0;
    vector<LL> pre(n + 1), suf(n + 2), fac(n), inv(n);
    pre[0] = fac[0] = inv[0] = suf[n + 1] = 1;
    for (int i = 1; i <= n; i++) pre[i] = (k - x[i]) % MOD * pre[i - 1] % MOD;
    for (int i = n; i >= 1; i--) suf[i] = (k - x[i]) % MOD * suf[i + 1] % MOD;
    for (int i = 1; i < n; i++) fac[i] = fac[i - 1] * i % MOD;
    inv[n - 1] = qpow(fac[n - 1], MOD - 2, MOD);
    for (int i = n - 2; ~i; i--) inv[i] = inv[i + 1] * (i + 1) % MOD;
    for (int i = 1; i <= n; i++)
        ans =
            (ans + ((n - i) % 2 ? -1ll : 1ll) * y[i] * pre[i - 1] % MOD *
                       suf[i + 1] % MOD * inv[i - 1] % MOD * inv[n - i] % MOD) %
            MOD;
    return (ans + MOD) % MOD;
}

signed main() {
    ios::sync_with_stdio(false);
    cin >> n >> m;
    for (int i = 1; i <= m + 2; i++) {
        x[i] = i;
        y[i] = (y[i - 1] + qpow((n - i) % MOD, m, MOD)) % MOD;
    }
    LL den = qpow(lagrange_cont(m + 2, x, y, n), MOD - 2, MOD);
    for (int i = 1; i <= m + 3; i++) {
        x[i] = i;
        y[i] = (y[i - 1] + qpow((n - i) % MOD, m, MOD) * i % MOD) % MOD;
    }
    LL num = lagrange_cont(m + 3, x, y, n);
    cout << num * den % MOD;
    return 0;
}

F Cumulative Sum

题目链接
做法见题解 AtCoder Beginner Contest 208F

G Cowmpany Compensation

题目链接

首先考虑用动态规划求解：设 $dp[i][j]$ 表示以 $i$ 为根，根节点工资不超过 $j$ 的方案数，从而有 $dp[i][j] = dp[i][j-1] + \prod \limits_{k∈Son_i}dp[k][j]$ 。

上述动态规划时间复杂度为 $O(nD)$ ，由于 $D$ 很大无法通过。实际上 $dp[1][D]$ 是 $D$ 的 $n$ 次多项式，因此我们只需要 $O(n^2)$ 求解 $dp[1][1..n+1]$ ，然后使用拉格朗日插值求解 $dp[1][D]$ 。

Cowmpany Compensation

#include <bits/stdc++.h>

#define debug(x) cerr << #x << " = " << x << endl

using namespace std;
typedef long long LL;

const int MAXN = 3005;
const int MOD = 1e9 + 7;
LL n, d;
LL x[MAXN], dp[MAXN][MAXN];
vector<int> G[MAXN];

LL qpow(LL a, LL b, LL MOD) {
    LL ans = 1, base = a;
    while (b) {
        if (b & 1) ans = ans * base % MOD;
        b >>= 1, base = base * base % MOD;
    }
    return ans;
}

LL lagrange_cont(int n, LL* x, LL* y, LL k) {
    LL ans = 0;
    vector<LL> pre(n + 1), suf(n + 2), fac(n), inv(n);
    pre[0] = fac[0] = inv[0] = suf[n + 1] = 1;
    for (int i = 1; i <= n; i++) pre[i] = (k - x[i]) % MOD * pre[i - 1] % MOD;
    for (int i = n; i >= 1; i--) suf[i] = (k - x[i]) % MOD * suf[i + 1] % MOD;
    for (int i = 1; i < n; i++) fac[i] = fac[i - 1] * i % MOD;
    inv[n - 1] = qpow(fac[n - 1], MOD - 2, MOD);
    for (int i = n - 2; ~i; i--) inv[i] = inv[i + 1] * (i + 1) % MOD;
    for (int i = 1; i <= n; i++)
        ans =
            (ans + ((n - i) % 2 ? -1ll : 1ll) * y[i] % MOD * pre[i - 1] % MOD *
                       suf[i + 1] % MOD * inv[i - 1] % MOD * inv[n - i] % MOD) %
            MOD;
    return (ans + MOD) % MOD;
}

void dfs(int now) {
    for (auto& to : G[now]) dfs(to);
    for (int i = 1; i <= n + 1; i++) {
        dp[now][i] = dp[now][i - 1];
        LL tmp = 1;
        for (auto& to : G[now]) tmp = tmp * dp[to][i] % MOD;
        dp[now][i] = (dp[now][i] + tmp) % MOD;
    }
}

signed main() {
    ios::sync_with_stdio(false);
    cin >> n >> d;
    for (int i = 2; i <= n; i++) {
        int p;
        cin >> p;
        G[p].push_back(i);
    }
    dfs(1);
    for (int i = 1; i <= n + 1; i++) x[i] = i;
    cout << lagrange_cont(n + 1, x, dp[1], d);
    return 0;
}

拉格朗日插值

求插值多项式系数

拉格朗日插值优化

横坐标连续时插值

重心拉格朗日插值

多项式多点插值与快速插值

例题

A 洛谷 拉格朗日插值

B LOJ 拉格朗日插值

C The Sum of the k-th Powers

D 教科书般的亵渎

E 排名估算

F Cumulative Sum

G Cowmpany Compensation

A 洛谷拉格朗日插值