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

太过爱你忘了你带给我的痛 2021-10-23 18:24 283阅读 0赞

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