[洛谷P5158]【模板】多项式快速插值
题目大意：有$n$个点$(x_i,y_i)$，求一个$n-1$次的多项式满足$f(x_i)\equiv y_i\pmod{998244353}$。$n\leqslant10^5$
题解：先有拉格朗日插值公式：
$$
\sum\limits_iy_i\prod\limits_{j\not=i}\dfrac{x-x_j}{x_i-x_j}
$$
发现，如果这么去做，至少是$O(n^2)$的，并不可以通过。考虑优化这个式子，把它拆开
$$
\begin{align*}
&\sum\limits_{i=1}^ny_i\prod\limits_{j\not=i}\dfrac{x-x_j}{x_i-x_j}\\
=&\sum\limits_{i=1}^n\dfrac{y_i}{\prod_{j\not=i}(x_i-x_j)}\prod\limits_{j\not=i}(x-x_j)
\end{align*}
$$
令$G(x)=\prod\limits_{i=1}^n(x-x_i)$，可知$\prod\limits_{j\not=i}(x_i-x_j)=\dfrac{G(x)}{x-x_i}$，因为分母为$0$，不可以直接求。由于上下在$x_i$处均为$0$，可使用洛必达法则：
$$
\lim\limits_{x\to x_i}\dfrac{G(x)}{x-x_i}=G'(x)
$$
故$\prod\limits_{j\not=i}(x_i-x_j)=G'(x_i)$，所以可以构造$G'(x)$然后多点求值求得$G'(x)$在每个$x_i$处的点值。此时，式子变成了：
$$
\sum\limits_{i=1}^n\dfrac{y_i}{G'(x_i)}\prod\limits_{j\not=i}(x-x_i)
$$
令$F_{L,R}$为$[L,R)$内的答案，$P_{L,R}=\prod\limits_{i=L}^{R-1}(x-x_i)$，则$F_{L,R}=P_{L,M}Q_{M,R}+P_{M,R}Q_{L,M}$
卡点：调试记录断点没有删除，然后调了好久
C++ Code：
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <iostream>
#include <vector>
#define mul(a, b) (static_cast<long long> (a) * (b) % mod)
#define _mul(a, b) (a = static_cast<long long> (a) * (b) % mod)
#define clr(A, l, r) (std::memset(A + (l), 0, (r) - (l) << 2))
#define cpy(A, B, len) (std::memcpy(A, B, (len) << 2))
#define cpyclr(A, B, len) (cpy(A, B, len), clr(A, len, lim))
const int maxn = 1 << 18, mod = 998244353;
typedef std::vector<int> VI;
inline void reduce(int &x) { x += x >> 31 & mod; }
namespace Math {
inline int pw(int base, int p) {
static int res;
for (res = 1; p; p >>= 1, _mul(base, base)) if (p & 1) _mul(res, base);
return res;
}
inline int inv(int x) { return pw(x, mod - 2); }
}
namespace Poly {
#define N maxn
int lim, s, rev[N], Wn[N];
inline void FFTINIT(int n) {
lim = 1, s = 0; while (lim < n) lim <<= 1, ++s;
for (int i = 1; i < lim; ++i) rev[i] = rev[i >> 1] >> 1 | (i & 1) << (s - 1);
const int t = Math::pw(3, (mod - 1) / lim);
Wn[lim >> 1] = 1;
for (int *i = Wn + 1 + (lim >> 1); i != Wn + lim; ++i) *i = mul(*(i - 1), t);
for (int i = lim >> 1; --i; ) Wn[i] = Wn[i << 1];
}
inline void INIT(int n) { lim = 1; while (lim < n) lim <<= 1; }
inline void FFT(int *A, const int op = 1) {
static unsigned long long t[N];
int shift = s - __builtin_ctz(lim);
for (int i = 0; i < lim; ++i) t[rev[i] >> shift] = A[i];
for (int mid = 1; mid < lim; mid <<= 1)
for (int i = 0; i < lim; i += mid << 1)
for (int j = 0; j < mid; ++j) {
const unsigned long long Y = t[i + j + mid] * Wn[j + mid] % mod;
t[i + j + mid] = t[i + j] - Y + mod, t[i + j] += Y;
}
for (int i = 0; i < lim; ++i) A[i] = t[i] % mod;
if (!op) {
const int ilim = Math::inv(lim);
for (int *i = A; i != A + lim; ++i) _mul(*i, ilim);
std::reverse(A + 1, A + lim);
}
}
void INV(int *A, int *B, int n) {
if (n == 1) { *B = Math::inv(*A); return ; }
static int C[N], D[N];
const int len = n + 1 >> 1;
INV(A, B, len), INIT(len * 3);
cpyclr(C, A, n), cpyclr(D, B, len);
FFT(C), FFT(D);
for (int i = 0; i < lim; ++i) D[i] = (2 - mul(C[i], D[i]) + mod) * D[i] % mod;
FFT(D, 0), cpy(B + len, D + len, n - len);
}
void DIF(int *A, int *B, int n) {
B[n - 1] = 0; for (int i = 1; i < n; ++i) B[i - 1] = mul(A[i], i);
}
void MUL(VI &res, VI P1, VI P2) {
static int A[N], B[N];
const int n = P1.size(), m = P2.size();
INIT(n + m - 1);
std::copy(P1.begin(), P1.end(), A), clr(A, n, lim);
std::copy(P2.begin(), P2.end(), B), clr(B, m, lim);
FFT(A), FFT(B);
for (int i = 0; i < lim; ++i) _mul(A[i], B[i]);
FFT(A, 0);
res.assign(A, A + n + m - 1);
}
int pos[N];
VI P[N << 1];
void DC_FFT(int rt, int l, int r) {
if (l == r) { P[rt] = { pos[l], 1 }; return ; }
const int mid = l + r >> 1, L = rt << 1, R = rt << 1 | 1;
DC_FFT(L, l, mid), DC_FFT(R, mid + 1, r);
MUL(P[rt], P[L], P[R]);
}
namespace Evaluation {
int res[N];
VI S[N << 1];
void DIV(int A, int n, int B, int m, int *F) {
const int len = n - m + 1;
static int C[N], D[N];
std::reverse_copy(S[A].begin(), S[A].end(), C);
std::reverse_copy(P[B].begin(), P[B].end(), D);
INV(D, F, len), INIT(len << 1);
clr(C, len, lim), clr(F, len, lim);
FFT(C), FFT(F);
for (int i = 0; i < lim; ++i) _mul(F[i], C[i]);
FFT(F, 0), std::reverse(F, F + len);
}
void __DIVMOD(int res, int A, int n, int B, int m) {
if (n < m) {
S[res].assign(S[A].begin(), S[A].end());
return ;
}
static int C[N], D[N];
DIV(A, n, B, m, C), INIT(n);
std::copy(P[B].begin(), P[B].end(), D);
clr(C, n - m + 1, lim), clr(D, m, lim);
FFT(C), FFT(D);
for (int i = 0; i < lim; ++i) _mul(C[i], D[i]);
FFT(C, 0);
for (int i = 0; i < m - 1; ++i) reduce(C[i] = S[A][i] - C[i]);
S[res].assign(C, C + m - 1);
}
void DIVMOD(int res, int A) {
int n = S[A].size(), m = P[res].size();
__DIVMOD(res, A, n, res, m);
}
void calc(int rt, int l, int r) {
if (l == r) { res[l] = S[rt][0]; return ; }
int mid = l + r >> 1;
DIVMOD(rt << 1, rt), DIVMOD(rt << 1 | 1, rt);
calc(rt << 1, l, mid), calc(rt << 1 | 1, mid + 1, r);
}
void EVAL(int *f, int n, int m, int *__pos, int *__res) {
for (int i = 1; i <= m; ++i) reduce(pos[i] = -__pos[i]);
DC_FFT(1, 1, m);
S[0].assign(f, f + n), DIVMOD(1, 0);
calc(1, 1, m), cpy(__res + 1, res + 1, m);
}
void EVAL_for_INTER(int *f, int n, int m, int *__res) {
S[0].assign(f, f + n), DIVMOD(1, 0);
calc(1, 1, m), cpy(__res + 1, res + 1, m);
}
}
using Evaluation::EVAL_for_INTER;
using Evaluation::EVAL;
namespace Interpolation {
int g[N], Y[N];
VI S[N << 1];
void calc(int rt, int l, int r) {
if (l == r) { S[rt] = { static_cast<int> mul(Y[l], Math::inv(g[l])) }; return ; }
const int mid = l + r >> 1, L = rt << 1, R = rt << 1 | 1;
calc(L, l, mid), calc(R, mid + 1, r);
static VI t;
MUL(S[rt], S[L], P[R]), MUL(t, S[R], P[L]);
for (int i = 0; i <= r - l; ++i) reduce(S[rt][i] += t[i] - mod);
}
void INTER(int *X, int *__Y, int n, int *res) {
static int C[N], D[N];
for (int i = 1; i <= n; ++i)
reduce(pos[i] = -X[i]), Y[i] = __Y[i];
DC_FFT(1, 1, n), std::copy(P[1].begin(), P[1].end(), C);
DIF(C, D, n + 1), EVAL_for_INTER(D, n, n, g);
calc(1, 1, n), std::copy(S[1].begin(), S[1].end(), res);
}
}
using Interpolation::INTER;
#undef N
}
int n, X[maxn], Y[maxn], res[maxn];
int main() {
std::ios::sync_with_stdio(false), std::cin.tie(0), std::cout.tie(0);
std::cin >> n;
Poly::FFTINIT(n + n);
for (int i = 1; i <= n; ++i) std::cin >> X[i] >> Y[i];
Poly::INTER(X, Y, n, res);
for (int i = 0; i < n; ++i) std::cout << res[i] << ' ';
std::cout << '\n';
return 0;
}
