拉格朗日插值

给定 $n$ 个点，我们可以确定唯一的最高 degree 为 $(n-1)$ 的多项式。

设多项式为 $f(x)$，第 $i$ 个点的坐标为 $(x_i,y_i)$，那么这个多项式在 $k$ 处的取值为：

$$f(k) = \sum\limits_{i=1}^ny_i \prod\limits_{i\neq j} \frac{k-x_j}{x_i-x_j}$$

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

拉格朗日插值板子 
const int maxn = 2005;

int n;
ll k;
ll x[maxn], y[maxn];

ll qpow(ll a, ll b) {
    ll res = 1;
    while (b) {
        if (b & 1)
            res = res * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return res;
}

ll solve(ll k) {
    ll ans = 0;
    for (int i = 1; i <= n; i++) {
        ll nu = 1, de = 1;
        for (int j = 1; j <= n; j++) {
            if (j == i) continue;
            nu = nu * (k - x[j] + mod) % mod;
            de = de * (x[i] - x[j] + mod) % mod;
        }
        de = qpow(de, mod-2);
        ans = (ans + y[i] * nu % mod * de % mod) % mod;
    }
    return ans;
}

int main() {
    cin >> n >> k;
    for (int i = 1; i <= n; i++) cin >> x[i] >> y[i];
    ll ans = solve(k);
    cout << ans << endl;
}

$x$ 坐标为连续整数的拉格朗日插值

如果 $x$ 坐标是连续的整数，那么我们可以在 $O(n)$ 求出 $f(k)$。

我们假设 $x_i = i$，那么：

$$f(k) = \sum\limits_{i=1}^ny_i \prod\limits_{i\neq j} \frac{k-x_j}{x_i-x_j}$$

可以转化为：

$$f(k) = \sum\limits_{i=1}^ny_i \prod\limits_{i\neq j} \frac{k-j}{i-j}$$

分类讨论 $j < i$ 和 $j > i$ 的情况：

$$f(k) = \sum\limits_{i=1}^ny_i \prod\limits_{j=1}^{i-1} \frac{k-j}{i-j} \prod\limits_{j=i+1}^n \frac{k-j}{i-j}$$

所以可得：

$$f(k) = \sum\limits_{i=1}^ny_i \frac{1}{(i-1)!} (k-1)(k-2)…(k-(i-1)) \frac{1}{-1*-2*…*(-(n-i))} (k-i)(k-(i+1))…(k-n)$$

化简得到：

$$f(k) = \sum\limits_{i=1}^n \frac{y_i}{(-1)^{n-i} *i! * (n-i)!} \prod\limits_{j=1}^{i-1}(k-j) \prod\limits_{j=i+1}^n(k-j)$$

• 注意，后面这部分不能变成 $\frac{\prod\limits_{j=1}^{n}(k-j)}{k-i}$，因为我们无法确定 $k-i \neq 0$。

对于 $\prod\limits_{j=1}^{i-1}(k-j) \prod\limits_{j=i+1}^n(k-j)$，我们只要预处理出来一个前缀积和后缀积即可。

所以对于每一项（每个 $i$），我们都可以 $O(1)$ 时间算出对应值，总复杂度 $O(n)$。

重心拉格朗日插值

用于解决动态加点的问题。

利用重心拉格朗日插值，每加入一个新的点，我们可以在 $O(n)$ 时间内算出新的多项式。

$$f(k) = \sum\limits_{i=1}^ny_i \prod\limits_{i\neq j}^n \frac{k-x_j}{x_i-x_j}$$

我们设 $g(k) = \prod\limits_{i=1}^n(k-x_i)$，则有：

$$f(k) = g(k)\sum\limits_{i=1}^n \frac{1}{k-x_i} \prod\limits_{i\neq j}^n \frac{y_i}{x_i-x_j}$$

设 $t_i = \frac{y_i}{\prod\limits_{j \neq i}^n (x_i-x_j)}$，则有：

$$f(k) = g(k) \sum\limits_{i=1}^n \frac{t_i}{k-x_i}$$

所以每次添加一个新的点 $(x_{n+1}, y_{n+1})$ 时：

  1. 重新计算一下所有的 $t_i = t_i * \frac{1}{x_i - x_{n+1}}, i \in [1,n]$。
  2. 计算 $t_{n+1}$。
  3. 计算 $g(k)$。

• 注意，在求 $f(k)$ 时需要先判断一下是否有 $k=x_i$，有的话直接返回 $y_i$。（这个特判仅需要在重心拉格朗日插值中进行）。

重心拉格朗日插值板子 
int Q;
ll x[maxn], y[maxn], t[maxn];
int n;

ll qpow(ll a, ll b) {
    ll res = 1;

    while (b) {
        if (b & 1)
            res = res * a % mod;

        a = a * a % mod;
        b >>= 1;
    }

    return res;
}

ll solve(ll k) {
    for (int i = 1; i <= n; i++) {
        if (x[i] == k)  // 需要特判是否 k = x[i]
            return y[i];
    }

    ll ans = 0;
    ll g = 1;
    ll sum = 0;

    for (int i = 1; i <= n; i++) {
        g = g * (k - x[i] + mod) % mod;
        sum = (sum + t[i] * qpow((k - x[i] + mod) % mod, mod - 2) % mod) % mod;
    }

    ans = g * sum % mod;
    return ans;
}

int main() {
    cin >> Q;
    while (Q--) {
        int op; cin >> op;

        if (op == 1) {  // 添加 (xx,yy) 的点
            ll xx, yy;
            cin >> xx >> yy;
            n++;
            x[n] = xx, y[n] = yy;
            ll de = 1;
            for (int j = 1; j < n; j++) {
                de = de * (x[n] - x[j] + mod) % mod;  // 更新 t[n]
                t[j] = t[j] * qpow((x[j] - x[n] + mod) % mod, mod - 2) % mod;
            }
            de = qpow(de, mod - 2);
            t[n] = y[n] * de % mod;
        } else {
            ll k; cin >> k;  // 求 f(k)
            cout << solve(k) << endl;
        }
    }
}

常用模型

  1. $\sum\limits_{i=1}^k i^n = 1^n+2^n+…k^n$ 是一个多项式 $f(k)$，其中 $deg(f)=n+1$（意味着需要 $n+2$ 个点进行插值）。

例1 CF622F The Sum of the k-th Powers

题意

给定 $k,n$，求：

$$\sum\limits_{i=1}^k i^n = 1^n+2^n+…k^n$$

其中，$0 \leq k \leq 10^9, 1 \leq n \leq 10^6$。

题解

上面说了这是一个多项式 $f(k)$，其中 $deg(f)=n+1$，意味着需要 $n+2$ 个点进行插值。

又因为 $x_i = i, y_i = \sum\limits_{j=1}^i j^n$，我们直接按照上述 $O(n)$ 的连续整数来插值即可。

代码 
#include <bits/stdc++.h>
using namespace std;
const int maxn = 1e6+5;

ll qpow(ll a, ll b) {
    ll res = 1;
    while (b) {
        if (b & 1)
            res = res * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return res;
}

ll k,n;
ll fac[maxn], inv_fac[maxn], t[maxn];
ll y[maxn], pre[maxn], suf[maxn];
int main() {
    fastio;
    cin >> k >> n;
    n += 2;

    fac[0] = inv_fac[0] = 1;
    pre[0] = suf[n+1] = 1;
    for (ll i = 1; i <= n; i++) {
        fac[i] = fac[i-1] * i % mod;
        y[i] = (y[i-1] + qpow(i, n-2)) % mod;
        pre[i] = pre[i-1] * (k - i + mod) % mod;
    }
    for (ll i = n; i >= 1; i--) suf[i] = suf[i+1] * (k - i + mod) % mod;
    inv_fac[n] = qpow(fac[n], mod-2);
    for (ll i = n-1; i >= 1; i--) {
        inv_fac[i] = inv_fac[i+1] * (i+1) % mod;
    }

    ll ans = 0;
    for (ll i = 1; i <= n; i++) {
        int flag = ((n-i) & 1 ? -1 : 1);
        ll res = y[i] % mod * inv_fac[i-1] % mod * inv_fac[n-i] % mod * pre[i-1] % mod * suf[i+1] % mod;
        res = (res * flag + mod) % mod;
        ans = (ans + res) % mod;
    }
    cout << ans << endl;
}