浅谈拉格朗日插值

问题描述

有一个 \(n\) 次多项式 \(f(x)\)。

已知 \(f(x_0) = y_0,f(x_1)=y_1,\dots,f(x_n) = y_n\)（显然 \(x\) 互不相同），还原原多项式 \(f(x)\)。

解法

我们可以看作 \(n + 1\) 个点 \((x_0, y_0),(x_1, y_1),\dots,(x_n, y_n)\)，求一条过这 \(n + 1\) 个点的 \(n\) 次函数。

存在性与唯一性

不会证所以感性理解。

两点唯一确定一条一次函数，三点唯一确定一条二次函数，感性理解 \(n + 1\) 个点就能唯一确定一条 \(n\) 次函数。

高斯消元法

设 \(f(x) = a_0x^{b_0} + a_1x^{b_1} + \dots +a_nx^{b_n}\)，那么可以列出以下 \(n\) 元 \(n\) 次方程

\[\left\{ \begin{array}{lr} a_0x_0^{b_0} + a_1x_0^{b_1} + \dots +a_nx_0^{b_n} = y_0, &\\ a_0x_1^{b_0} + a_1x_0^{b_1} + \dots +a_nx_1^{b_n} = y_1, &\\ \dots, &\\ a_0x_n^{b_0} + a_1x_n^{b_1} + \dots +a_nx_n^{b_n} = y_n, &\\ \end{array} \right. \]

直接高斯消元硬解即可。

时间复杂度 \(\mathcal O(n^3)\)。

拉格朗日插值法

我们考虑构造这个函数 \(f\)。

我们把 \(f\) 拆成若干个函数的和：

\[f(x) = \sum_{i = 0} ^ {n}f_i(x) \]

那么我们只需要构造出满足以下条件的 \(n + 1\) 个 \(f_i\) 函数就能满足题意：

\[\left\{ \begin{array}{lr} f_i(x_i) = y_i \\ f_i(x_j) = 0 ~~ (i \ne j) \end{array} \right. \]

先考虑满足第二个条件。

\[f_i(x) = \prod_{0\le j \le n,j \ne i}(x - x_j) \]

那么如果 \(x \ne x_i\)，则必有一个 \(j\) 满足 \(x = x_j\)，此时 \(x - x_j = 0\) 则原函数值为 \(0\)。

然而当 \(x\) 取 \(x_i\) 我们并不能确保 \(f_i(x_i) = y_i\)。

我们尝试列出来 \(x = x_i\) 时该函数的取值。

\[\prod_{0\le j \le n,j \ne i}(x_i - x_j) \]

然后我们惊奇（好像并不）地发现这是一个常值！

于是我们可以将原函数除以上述值并乘上一个 \(y_i\)，那么就满足条件了，整理一下得：

\[f_i(x) = y_i\prod_{0\le j \le n,j \ne i}\frac{x - x_j}{x_i - x_j} \]

那么：

\[f(x) = \sum_{i = 0}^{n}y_i\prod_{0\le j \le n,j \ne i}\frac{x - x_j}{x_i - x_j} \]

时间复杂度：\(\mathcal O(n^2)\)。

例题

P4781

即模板题。

没什么区别，按照上文所述模拟即可。

#include <bits/stdc++.h>
#define int long long
#define pii pair<int, int>
#define L(i, a, b) for (register int i = (a); i <= (b); i++)
#define R(i, a, b) for (register int i = (a); i >= (b); i--)
#define FRE(x) freopen(x ".in", "r", stdin), freopen(x ".out", "w", stdout)
#define ALL(x) x.begin(), x.end()
using namespace std;

inline void cmax(int& x, int c) {
    x = max(x, c);
}
inline void cmin(int& x, int c) {
    x = min(x, c);
}

int _test_ = 1;

const int N = 2e3 + 5, mod = 998244353;

int n, k, x[N], y[N];

int qpow(int x, int y) {
    if (y == 0)
        return 1;
    int t = qpow(x, y / 2);
    t = (t * t) % mod;
    if (y & 1)
        return (t * x) % mod;
    return t;
}

void init() {}

void clear() {}

void solve() {
    cin >> n >> k;
    L(i, 1, n) cin >> x[i] >> y[i];
    int ans = 0;
    L(i, 1, n) {
        int tmp = y[i] % mod;
        L(j, 1, n) {
            if (i == j)
                continue;
            tmp = (tmp * (k - x[j])) % mod;
            tmp = (tmp * qpow(x[i] - x[j], mod - 2)) % mod;
        }
        ans = (ans + tmp + mod) % mod;  // 这里要加上 mod，因为我们不能保证 tmp 是正数
    }
    cout << ans;
}

signed main() {
    ios::sync_with_stdio(0);
    cin.tie(0), cout.tie(0);
    // cin >> _test_;
    init();
    while (_test_--) {
        clear();
        solve();
    }
    return 0;
}

CF622F

题意

求 \(\sum_{i = 1}^{n} i^k\)，\(n \le 10^9, k \le 10^6\)。

解法

朴素做法枚举快速幂，时间复杂度 \(\mathcal O(n \log k)\)。

---

感性理解一下，当 \(k\) 为 \(1\) 的时候，有公式：\(\frac{n \times (n + 1)}{2}\) 是一个二次多项式。

\(k\) 为 \(2\) 时，有公式：\(\frac{n \times (n + 1) \times (2n + 1)}{6}\)，是一个三次多项式。

\(k\) 为 \(3\) 时，有公式：\(\left[ \frac{n \times (n + 1)}{2} \right]^2\)，是一个四次多项式。

那么 \(\sum_{i = 1}^{n} i^k\) 就是一个 \(k + 1\) 次多项式。

取 \(k + 2\) 个点做拉格朗日插值法可得其关于 \(n\) 的多项式 \(f(n)\)。

把 \(n\) 带进去即可，时间复杂度 \(\mathcal O(k^2\log k)\)。

---

然而，复杂度仍然是不对的。

观察到我们取的 \(k + 2\) 个点是不确定的，所以我们可以钦定第 \(i\) 个点取 \((i, \sum_{j = 1} ^ {n}j^k)\)，令其表示为 \((i, x_i)\)。

那么我们列出最终的多项式：

\[f(n) = \sum_{i = 1} ^ {k + 2}x_i\prod_{1 \le j \le k + 2, j \ne i}\frac{n - j}{i -j} \]

我们发现后面这一大坨非常可算啊。

只看后面的这一部分：

\[\begin{align*} &\\ \prod_{1 \le j \le k + 2, j \ne i}\frac{n - j}{i -j} &= \frac{\prod_{1 \le j \le k + 2, j \ne i}(n - j)}{\prod_{1 \le j \le k + 2, j \ne i}(i -j)}\\ &= \frac{\prod_{1 \le j \le k + 2, j \ne i}(n - j)}{\prod_{1 \le j \le i - 1}(i - j) \times \prod_{i + 1 \le j \le k + 2}(i - j)}\\ &= \frac{\prod_{1 \le j \le k + 2, j \ne i}(n - j)}{(i - 1)\times(i - 2) \times \dots \times 1 \times (-1) \times (-2) \times \dots \times (i - k - 2)} \\ &= \frac{\prod_{1 \le j \le k + 2, j \ne i}(n - j)}{(i - 1)! \times (k + 2 - i) ! \times (-1)^{k + 2 - i}}\\ &= \frac{\prod_{1 \le j \le i - 1}(n - j) \times \prod_{i + 1 \le j \le k + 2}(n - j)}{(i - 1)! \times (k + 2 - i) ! \times (-1)^{k + 2 - i}} \end{align*} \]

分母可以通过预处理解决，而分子则可以拆成前缀、后缀积。

然后 \(-1\) 的次方真的不用快速幂啊。

那么时间复杂度就来到了 \(\mathcal O(k\log k)\)，足以通过此题。

#include <bits/stdc++.h>
#define int long long
#define pii pair<int, int>
#define L(i, a, b) for (register int i = (a); i <= (b); i++)
#define R(i, a, b) for (register int i = (a); i >= (b); i--)
#define FRE(x) freopen(x ".in", "r", stdin), freopen(x ".out", "w", stdout)
#define ALL(x) x.begin(), x.end()
using namespace std;

inline void cmax(int& x, int c) {
    x = max(x, c);
}
inline void cmin(int& x, int c) {
    x = min(x, c);
}

int _test_ = 1;

const int N = 1e6 + 7, mod = 1e9 + 7;

int n, k, pret[N], suft[N], f1[N], fac[N];

int qpow(int x, int y) {
    if (y == 0)
        return 1;
    int t = qpow(x, y / 2);
    t = (t * t) % mod;
    if (y & 1)
        return (t * x) % mod;
    return t;
}

void init() {}

void clear() {}

void solve() {
    cin >> n >> k;
    pret[0] = 1;  // 前缀积
    suft[k + 3] = 1;  // 后缀积
    fac[0] = 1;  // 阶乘
    f1[0] = 1;  // -1 的次方（其实可以放在后面）
    L(i, 1, k + 2) {
        f1[i] = f1[i - 1] * (-1);
        pret[i] = (pret[i - 1] * (n - i)) % mod;
        fac[i] = (fac[i - 1] * i) % mod;
    }
    R(i, k + 2, 1) {
        suft[i] = (suft[i + 1] * (n - i)) % mod;
    }
// 以上为预处理
    int x = 0, ans = 0;
    L(i, 1, k + 2) {
        x = (x + qpow(i, k)) % mod;
        ans = (ans + ((x * (pret[i - 1] * suft[i + 1] % mod) % mod) % mod *
            qpow((fac[i - 1] * fac[k + 2 - i] % mod * f1[k + 2 - i] % mod) % mod,
                mod - 2)) %
            mod) %
        mod;
    }
    cout << (ans + mod) % mod;  // 要加 mod 不然会 WA on#9
}

signed main() {
    ios::sync_with_stdio(0);
    cin.tie(0), cout.tie(0);
    // cin >> _test_;
    init();
    while (_test_--) {
        clear();
        solve();
    }
    return 0;
}