Note/Solution -「洛谷 P5158」「模板」多项式快速插值

$\mathcal{Description}$

  Link.

  给定 $n$ 个点 $(x_i,y_i)$，求一个不超过 $n-1$ 次的多项式 $f(x)$，使得 $f(x_i)\equiv y_i\pmod{998244353}$。

  $n\le10^5$。

$\mathcal{Solution}$

  摆出 Lagrange 插值的式子：

$$f(z)=\sum_{i=1}^ny_i\prod_{j\neq i}\frac{z-x_j}{x_i-x_j}.$$

现在的问题是分母上的 $\prod_{j\neq i}(x_i-x_j)$ 不好求。若令

$$g(x)=\prod_{i=1}^n(x-x_i),$$

则

$$f(z)=\sum_{i=1}^ny_i\left(\lim_{x\rightarrow x_i}\frac{g(x)}{x-x_i}\right)\prod_{i\neq j}(z-x_j).$$

中间的 $\lim$ 可以直接洛出来啊，也可以构造 $\lim_{x\rightarrow x_i}\frac{g(x)}{x-x_i}=\lim_{x\rightarrow x_i}\frac{g(x)-g(x_i)}{x-x_i}$，整理得到

$$f(z)=\sum_{i=1}^n\frac{y_i}{g'(x_i)}\prod_{i\neq j}(z-x_j).$$

先分治求出 $g$，然后多点求值求得 $g'(x_i)$，再分治求出 $f$ 即可。注意求 $g$ 的过程量 $\prod_{i=l}^r(z-x_i)$ 翻转系数就得到多点求值要用的 $\prod_{i=l}^r(1-x_iz)$，可以节约一点常数。最终复杂度 $\mathcal O(n\log^2 n)$。

$\mathcal{Code}$

/*+Rainybunny+*/

#include <bits/stdc++.h>

#define rep(i, l, r) for (int i = l, rep##i = r; i <= rep##i; ++i)
#define per(i, r, l) for (int i = r, per##i = l; i >= per##i; --i)

typedef std::vector<int> Poly;

const int MAXN = 1 << 18, MOD = 998244353;
int n, x[MAXN + 5], y[MAXN + 5];
Poly X[MAXN << 2];

inline int mul(const int u, const int v) { return 1ll * u * v % MOD; }
inline int sub(int u, const int v) { return (u -= v) < 0 ? u + MOD : u; }
inline int add(int u, const int v) { return (u += v) < MOD ? u : u - MOD; }
inline int mpow(int u, int v) {
    int ret = 1;
    for (; v; u = mul(u, u), v >>= 1) ret = mul(ret, v & 1 ? u : 1);
    return ret;
}

namespace PolyOper {

const int G = 3;
int omega[19][MAXN + 5];

inline void init() {
    rep (i, 1, 18) {
        int* wi = omega[i];
        wi[0] = 1, wi[1] = mpow(G, MOD - 1 >> i);
        rep (j, 2, (1 << i) - 1) wi[j] = mul(wi[j - 1], wi[1]);
    }
}

inline void ntt(Poly& u, const int tp) {
    static int rev[MAXN + 5]; int n = u.size();
    rep (i, 0, n - 1) rev[i] = rev[i >> 1] >> 1 | (i & 1) * n >> 1;
    rep (i, 0, n - 1) if (i < rev[i]) std::swap(u[i], u[rev[i]]);
    for (int i = 1, stp = 1; stp < n; ++i, stp <<= 1) {
        int* wi = omega[i];
        for (int j = 0; j < n; j += stp << 1) {
            rep (k, j, j + stp - 1) {
                int ev = u[k], ov = mul(wi[k - j], u[k + stp]);
                u[k] = add(ev, ov), u[k + stp] = sub(ev, ov);
            }
        }
    }
    if (!~tp) {
        int inv = mpow(n, MOD - 2);
        std::reverse(u.begin() + 1, u.end());
        for (int& a: u) a = mul(a, inv);
    }
}

inline Poly padd(Poly u, Poly v) {
    if (u.size() < v.size()) u.swap(v);
    rep (i, 0, int(v.size()) - 1) u[i] = add(u[i], v[i]);
    return u;
}

inline Poly pmul(Poly u, Poly v) {
    int res = u.size() + v.size() - 1, len = 1;
    while (len < res) len <<= 1;
    u.resize(len), v.resize(len);
    ntt(u, 1), ntt(v, 1);
    rep (i, 0, len - 1) u[i] = mul(u[i], v[i]);
    ntt(u, -1);
    return u.resize(res), u;
}

inline Poly pmulT(Poly u, Poly v) {
    int n = u.size(), m = v.size();
    std::reverse(v.begin(), v.end()), v = pmul(u, v);
    rep (i, 0, n - 1) u[i] = v[i + m - 1];
    return u;
}

inline void pinv(const int n, const Poly& u, Poly& r) {
    if (n == 1) return void(r = { { mpow(u[0], MOD - 2) } });
    static Poly tmp; pinv(n >> 1, u, r);
    tmp.resize(n << 1), r.resize(n << 1);
    rep (i, 0, n - 1) tmp[i] = i < u.size() ? u[i] : 0;
    rep (i, n, (n << 1) - 1) tmp[i] = 0;
    ntt(r, 1), ntt(tmp, 1);
    rep (i, 0, (n << 1) - 1) r[i] = mul(r[i], sub(2, mul(tmp[i], r[i])));
    ntt(r, -1), r.resize(n);
}

} // namespace PolyOper.

inline void init(const int u, const int l, const int r) {
    if (l == r) return void(X[u] = { { 1, sub(0, x[l]) } });
    int mid = l + r >> 1;
    init(u << 1, l, mid), init(u << 1 | 1, mid + 1, r);
    X[u] = PolyOper::pmul(X[u << 1], X[u << 1 | 1]);
}

inline void calcG(const int u, const int l, const int r, Poly F) {
    F.resize(r - l + 1);
    if (l == r) return void(y[l] = mul(y[l], mpow(F[0], MOD - 2)));
    int mid = l + r >> 1;
    calcG(u << 1, l, mid, PolyOper::pmulT(F, X[u << 1 | 1]));
    calcG(u << 1 | 1, mid + 1, r, PolyOper::pmulT(F, X[u << 1]));
}

inline Poly calcF(const int u, const int l, const int r) {
    std::reverse(X[u].begin(), X[u].end());
    if (l == r) return { { y[l] } };
    int mid = l + r >> 1;
    Poly &&p(calcF(u << 1, l, mid)), &&q(calcF(u << 1 | 1, mid + 1, r));
    return PolyOper::padd(PolyOper::pmul(p, X[u << 1 | 1]),
      PolyOper::pmul(q, X[u << 1]));
}

int main() {
    scanf("%d", &n);
    rep (i, 0, n - 1) scanf("%d %d", &x[i], &y[i]);

    PolyOper::init(), init(1, 0, n - 1);
    
    int len = 1; while (len < n) len <<= 1;
    Poly T; PolyOper::pinv(len, X[1], T);
    Poly Q(X[1]); std::reverse(Q.begin(), Q.end());
    rep (i, 0, n - 1) Q[i] = mul(i + 1, Q[i + 1]);
    Q.resize(n);

    calcG(1, 0, n - 1, PolyOper::pmulT(Q, T));
    Poly&& ans = calcF(1, 0, n - 1);
    for (int u: ans) printf("%d ", u);
    putchar('\n');
    return 0;
}