[洛谷P5158]【模板】多项式快速插值

最新推荐文章于 2021-05-02 16:12:30 发布

weixin_30610755

最新推荐文章于 2021-05-02 16:12:30 发布

阅读量117

点赞数

原文链接：http://www.cnblogs.com/Memory-of-winter/p/11290246.html

版权

题目大意：有$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;
}

转载于:https://www.cnblogs.com/Memory-of-winter/p/11290246.html