Submit Info #44829

Problem Lang User Status Time Memory
Polynomial Interpolation cpp mnbvmar AC 2530 ms 48.97 MiB

ケース詳細
Name Status Time Memory
example_00 AC 0 ms 0.48 MiB
example_01 AC 1 ms 0.68 MiB
max_random_00 AC 2530 ms 48.97 MiB
max_random_01 AC 2519 ms 48.97 MiB
random_00 AC 2405 ms 47.29 MiB
random_01 AC 1251 ms 26.53 MiB
random_02 AC 925 ms 19.68 MiB

// Follows https://github.com/cuber2460/acmlib07/blob/master/code/finaly/code/FFT.cpp #include <bits/stdc++.h> using namespace std; #define sim template < class c #define ris return * this #define dor > debug & operator << #define eni(x) sim > typename enable_if<sizeof dud<c>(0) x 1, debug&>::type operator<<(c i) { sim > struct rge { c b, e; }; sim > rge<c> range(c i, c j) { return rge<c>{i, j}; } sim > auto dud(c* x) -> decltype(cerr << *x, 0); sim > char dud(...); struct debug { #ifdef LOCAL ~debug() { cerr << endl; } eni(!=) cerr << boolalpha << i; ris; } eni(==) ris << range(begin(i), end(i)); } sim, class b dor(pair < b, c > d) { ris << "(" << d.first << ", " << d.second << ")"; } sim, class b, class z dor(tuple < b, c, z > d) { ris << "(" << get<0>(d) << ", " << get<1>(d) << ", " << get<2>(d) << ")"; } sim dor(rge<c> d) { *this << "["; for (auto it = d.b; it != d.e; ++it) { *this << ", " + 2 * (it == d.b) << *it; } ris << "]"; } #else sim dor(const c&) { ris; } #endif }; #define imie(...) " [" << #__VA_ARGS__ ": " << (__VA_ARGS__) << "] " using ll = long long; using vi = vector<int>; #define SZ(x) ((int)(x).size()) #define ALL(x) (x).begin(), (x).end() // If [mod] needs to be variable, move it to the struct and initialize // it in the beginning of the constructor. template <int mod> struct FormalSeries { private: int maxn; int root; static inline void add_mod(int &a, int b) { a += b; if (a >= mod) { a -= mod; } } static inline void sub_mod(int &a, int b) { a -= b; if (a < 0) { a += mod; } } static inline int mul_mod(int a, int b) { return (int)((ll)a * b % mod); } static inline int pow_mod(int a, int n) { int r = 1; while (n) { if (n & 1) { r = mul_mod(r, a); } n >>= 1; a = mul_mod(a, a); } return r; } static inline int inv_mod(int a) { return pow_mod(a, mod - 2); } static inline int div_mod(int a, int b) { return mul_mod(a, inv_mod(b)); } vector<int> FindPrimeFactors(int num) { vector<int> ans; for (int p = 2; p * p <= num; ++p) { if (num % p == 0) { ans.push_back(p); while (num % p == 0) { num /= p; } } } if (num > 1) { ans.push_back(num); } return ans; } void FindGenerator() { vector<int> prime_factors = FindPrimeFactors(mod - 1); while (true) { const int g = rand() % (mod - 1) + 1; bool ok = true; for (int p : prime_factors) { if (pow_mod(g, mod / p) == 1) { ok = false; break; } } if (ok) { root = pow_mod(g, mod / maxn); break; } } } void DFT(vi &a) const { const int n = a.size(); for (int i = 1, k = 0; i < n; ++i) { for (int bit = n / 2; (k ^= bit) < bit; bit /= 2);;; if (i < k) { swap(a[i], a[k]); } } for (int len = 1, who = 0; len < n; len *= 2, ++who) { static vi t[32]; vi &om = t[who]; if (om.empty()) { om.resize(len + 1); if (who == 0) { om[0] = 1; om[1] = mod - 1; } else { const int stepg = pow_mod(root, maxn / (2 * len)); for (int i = 0; i < len; ++i) { om[i] = t[who - 1][i / 2]; if (i & 1) { om[i] = mul_mod(om[i], stepg); } } } om.back() = om[0]; } for (int i = 0; i < n; i += 2 * len) { for (int k = 0; k < len; ++k) { const int x = a[i + k]; const int y = mul_mod(a[i + k + len], om[k]); add_mod(a[i + k], y); a[i + k + len] = x; sub_mod(a[i + k + len], y); } } } } // Given n, computes some square root of n mod [mod], or -1 if no such // square root exists. If n != 0 and square root [R] exists, the other root // is given by [mod - R]. // The algorithm uses randomness; the choice of the square root may be // non-deterministic. // The implementation uses the Tonelli-Shanks algorithm. int SqrtMod(int n) const { if (n == 0) { return 0; } if (pow_mod(n, (mod - 1) / 2) != 1) { return -1; } int s; int q; for (s = 0, q = mod - 1; q % 2 == 0; s++, q /= 2); if (s == 1) { return pow_mod(n, (mod + 1) / 4); } int z; do { z = rand() % (mod - 1) + 1; } while(pow_mod(z, (mod - 1) / 2) == 1); int c = pow_mod(z, q); int r = pow_mod(n, (q + 1) / 2); int t = pow_mod(n, q); int m = s; while (true) { if (t == 1) { return r; } int i; for (i = 1; pow_mod(t, 1 << i) != 1; ++i); int b = pow_mod(c, 1 << (m - i - 1)); r = mul_mod(r, b); c = mul_mod(b, b); t = mul_mod(t, c); m = i; } } public: FormalSeries(int max_terms) { // (max_terms = max_degree + 1) maxn = 1; while (maxn < max_terms) { maxn *= 2; } assert((mod - 1) % maxn == 0 && "Not enough powers of 2 in mod-1"); FindGenerator(); } vi Mul(const vi &a, const vi &b) const { if (a.empty() || b.empty()) { return {}; } int n = a.size() + b.size(); const int ans_size = n - 1; assert(n <= maxn); while (n & (n - 1)) { ++n; } static vi X, Y; X = a; Y = b; X.resize(n); Y.resize(n); DFT(X); DFT(Y); for (int i = 0; i < n; ++i) { X[i] = mul_mod(X[i], Y[i]); } reverse(X.begin() + 1, X.end()); const int div = inv_mod(n); for (int &x : X) { x = mul_mod(x, div); } DFT(X); X.resize(ans_size); return X; } vi MakeConst(int x) const { assert(x >= 0); return vi{x % mod}; } void NegSelf(vi &A) const { for (int &x : A) { if (x) { x = mod - x; } } } void TrimSelf(vi &A, int num_terms) const { if (SZ(A) > num_terms) { A.resize(num_terms); } } vi Trim(vi A, int num_terms) const { TrimSelf(A, num_terms); return A; } vi Add(const vi &A, const vi &B) const { vi res = A; if (SZ(res) < SZ(B)) { res.resize(SZ(B)); } for (int i = 0; i < SZ(B); ++i) { add_mod(res[i], B[i]); } return res; } vi Sub(const vi &A, const vi &B) const { vi res = A; if (SZ(res) < SZ(B)) { res.resize(SZ(B)); } for (int i = 0; i < SZ(B); ++i) { sub_mod(res[i], B[i]); } return res; } // Returns the derivative of [Q]. [Q] must be nonempty. // Time complexity: O(n). vi Differentiate(vi Q) const { assert(!Q.empty()); const int n = SZ(Q); for (int i = 1; i < n; ++i) { Q[i - 1] = mul_mod(i, Q[i]); } Q.pop_back(); return Q; } // Returns the integral of [Q], with 0 as the lowest term of the result. // Time complexity: O(n). vi Integrate(const vi &Q) const { if (Q.empty()) { return {0}; } const int n = SZ(Q); vi P(n + 1); // Make inverses of 1..n. P[1] = 1; for (int i = 2; i <= n; ++i) { P[i] = mul_mod(mod - mod / i, P[mod % i]); } // Now, multiply them pointwise with Q. for (int i = 1; i <= n; ++i) { P[i] = mul_mod(P[i], Q[i - 1]); } return P; } // Given Q, |Q|=n, finds P, |P|=n such that P*Q = 1 (mod x^n). // If more than [n] terms are needed, simply extend [Q] by zeroes. // Time complexity: O(n log n). // // By Newton's method (x_{n+1} = x_n - f(x_n) / f'(x_n)): // Given Q. Find P, such that P = 1 / Q. // f(P) := Q - 1/P. Then, P_0 = inv(Q[0]), P_{n+1} = P_n * (2 - Q * P_n) mod x^{2^{n+1}}. // But Q * P_n = 1 + x^{2^n}*T (mod x^{2^{n+1}}). // So compute T, and then P_{n+1} = P_n - x^{2^n} * P_n * T (mod x^{2^{n+1}}). vi Inv(const vi &Q) const { const int n = SZ(Q); assert(n > 0 && Q[0] != 0); assert(n <= maxn / 2); auto P = MakeConst(inv_mod(Q[0])); while (SZ(P) < n) { const int sz = SZ(P) * 2; vi Pcopy = P; P.resize(sz); Pcopy.resize(sz); DFT(Pcopy); auto cur = Trim(Q, sz); cur.resize(sz); DFT(cur); for (int i = 0; i < sz; ++i) { cur[i] = mul_mod(Pcopy[i], cur[i]); } reverse(cur.begin() + 1, cur.end()); DFT(cur); const int inv_sz = inv_mod(sz); for (int i = 0; i < sz / 2; ++i) { cur[i] = mul_mod(cur[i + sz / 2], inv_sz); cur[i + sz / 2] = 0; } // Now, [cur] stores T. DFT(cur); for (int i = 0; i < sz; ++i) { cur[i] = mul_mod(Pcopy[i], cur[i]); } reverse(cur.begin() + 1, cur.end()); DFT(cur); for (int i = 0; i < sz / 2; ++i) { cur[i] = mul_mod(cur[i], inv_sz); } // Now, sz/2 terms of [cur] store P_n * T. for (int i = 0; i < sz / 2; ++i) { sub_mod(P[i + sz / 2], cur[i]); } // Or easier, slower variant of the loop (replace everything in the loop): // const int sz = SZ(P) * 2; // auto cur = Mul(Trim(Q, sz), P); // TrimSelf(cur, sz); // NegSelf(cur); // add_mod(cur[0], 2); // P = Trim(Mul(cur, P), sz); } TrimSelf(P, n); return P; } // Given Q, |Q|=n, finds P, |P|=n such that P = log(Q) (mod x^n). // If more than [n] terms are needed, simply extend [Q] by zeroes. // Asserts that Q[0] = 1 (otherwise, the operation is not in Z_p.) // Time complexity: O(n log n). vi Log(const vi &Q) const { // P = log(Q) -> P' = Q' / Q. Just compute Q' / Q and integrate. const int n = Q.size(); assert(n > 0 && Q[0] == 1); assert(n <= maxn / 2); auto Qdif = Differentiate(Q); return Integrate(Trim(Mul(Qdif, Inv(Q)), n - 1)); } // Given Q, |Q|=n, finds P, |P|=n such that P = exp(Q) (mod x^n). // If more than [n] terms are needed, simply extend [Q] by zeroes. // Asserts that Q[0] = 0 (otherwise, the operation is not in Z_p). // Time complexity: O(n log n). vi Exp(const vi &Q) const { const int n = SZ(Q); assert(n > 0 && Q[0] == 0); assert(n <= maxn / 2); auto P = MakeConst(1); while (SZ(P) < n) { const int sz = SZ(P); P.resize(sz * 2); auto Plog = Log(P); TrimSelf(P, sz); NegSelf(Plog); for (int i = 0; i < SZ(Q) && i < sz * 2; ++i) { add_mod(Plog[i], Q[i]); } add_mod(Plog[0], 1); P = Trim(Mul(P, Plog), sz * 2); } TrimSelf(P, n); return P; } // Given Q, |Q|=n, finds P, |P|=n, such that P = Q^n. // If more terms are needed, simply extend [Q] by zeroes. // This method works for any formal series, even if Q[0] = 0. // Time complexity: O((n log n) + log M). vi Pow(vi Q, int M) const { if (Q.empty()) { return {}; } const int n = SZ(Q); int pref_zero = 0; while (pref_zero < n && !Q[pref_zero]) { ++pref_zero; } if ((ll)pref_zero * M >= n) { return vi(n); } Q.erase(Q.begin(), Q.begin() + pref_zero); const int elem = Q[0]; const int inv_elem = inv_mod(elem); for (int &x : Q) { x = mul_mod(x, inv_elem); } assert(Q[0] == 1); Q = Log(Q); for (int &x : Q) { x = mul_mod(x, M); } Q = Exp(Q); const int pwr_elem = pow_mod(elem, M); for (int &x : Q) { x = mul_mod(x, pwr_elem); } Q.resize(n - pref_zero * M); Q.insert(Q.begin(), pref_zero * M, 0); return Q; } // Given Q, |Q|=n, finds some P, |P|=n, such that P^2 = Q. // If more terms are needed, simply extend [Q] by zeroes. // If no such P exists, returns an empty vector. // This method works for any formal series, even if Q[0] = 0. // Time complexity: O(n log n), but with a rather large constant // (could be faster if we instead used a hand-crafted Newton method for sqrt). vi Sqrt(vi Q) const { if (Q.empty()) { return {}; } const int n = SZ(Q); int pref_zero = 0; while (pref_zero < n && !Q[pref_zero]) { ++pref_zero; } if (pref_zero == n) { return vi(n); } if (pref_zero % 2 == 1) { return {}; } const int s2 = Q[pref_zero]; const int s1 = SqrtMod(s2); if (s1 == -1) { return {}; } const int inv_s2 = inv_mod(s2); Q.erase(Q.begin(), Q.begin() + pref_zero); for (int &x : Q) { x = mul_mod(x, inv_s2); } assert(Q[0] == 1); Q.resize(n - pref_zero / 2); Q = Log(Q); for (int &x : Q) { x = mul_mod(x, (mod + 1) / 2); } Q = Exp(Q); for (int &x : Q) { x = mul_mod(x, s1); } Q.insert(Q.begin(), pref_zero / 2, 0); return Q; } // Given P, Q, |P|=n, |Q|=m, Q nonzero, finds some polynomial X such that // |P - Q * X| < m. // Then, P - Q * X is also P mod Q. // If fails (for Q=0), returns the empty vector. // Time complexity: O(n log n). vi DivFloor(vi P, vi Q) const { while (!P.empty() && P.back() == 0) { P.pop_back(); } while (!Q.empty() && Q.back() == 0) { Q.pop_back(); } if (Q.empty()) { return {}; } const int n = SZ(P); const int m = SZ(Q); if (n < m) { return {0}; } reverse(P.begin(), P.end()); P.resize(n - m + 1); reverse(Q.begin(), Q.end()); Q.resize(n - m + 1); vi D = Trim(Mul(P, Inv(Q)), n - m + 1); reverse(D.begin(), D.end()); return D; } // Given P, Q, |P|=n, |Q|=m, Q nonzero, finds polynomials X, Y such that // Y = P - Q * X, and |Y| < m. // If fails (for Q=0), returns two empty vectors. // Time complexity: O(n log n). pair<vi, vi> DivMod(vi P, const vi &Q) const { auto X = DivFloor(P, Q); if (X.empty()) { return {{}, {}}; } if (X == vi{0}) { return {{0}, P}; } auto T = Mul(X, Q); assert(SZ(T) == SZ(P)); for (int i = 0; i < SZ(Q); ++i) { sub_mod(P[i], T[i]); } P.resize(SZ(Q) - 1); return {X, P}; } // Given polynomial [P], |P|=n, and set [xs], |xs|=m, evaluates P // in each point x \in xs. // Time complexity: O(n log n + m log^2 m), space complexity: O(n + m log m). vi MultipointEval(const vi &P, const vi &xs) { const int n = SZ(P); const int m = SZ(xs); if (m == 0) { return {}; } vector<vi> mod_polys(2 * m); for (int i = 0; i < m; ++i) { mod_polys[m + i] = {(mod - xs[i]) % mod, 1}; } for (int i = m - 1; i > 0; --i) { mod_polys[i] = Mul(mod_polys[2 * i], mod_polys[2 * i + 1]); } mod_polys[1] = DivMod(P, mod_polys[1]).second; for (int i = 2; i < 2 * m; ++i) { mod_polys[i] = DivMod(mod_polys[i / 2], mod_polys[i]).second; } vi answer(m); for (int i = 0; i < m; ++i) { assert(SZ(mod_polys[m + i]) == 1); answer[i] = mod_polys[m + i][0]; } return answer; } // Given sets [xs], [ys], |xs|=|ys|=n, finds a polynomial P of degree // <n such that P(xs[i]) = ys[i] for each i. // All [xs] must be different, but this method does not verify that. // Time complexity: O(n log^2 n), space complexity: O(n log n). vi Interpolation(const vi &xs, const vi &ys) { const int n = SZ(xs); vector<vi> mod_polys(2 * n), coef_polys(2 * n); for (int i = 0; i < n; ++i) { mod_polys[n + i] = {(mod - xs[i]) % mod, 1}; } for (int i = n - 1; i > 0; --i) { mod_polys[i] = Mul(mod_polys[2 * i], mod_polys[2 * i + 1]); } coef_polys[1] = {1}; for (int i = 1; i < n; ++i) { coef_polys[2 * i] = DivMod(Mul(coef_polys[i], mod_polys[2 * i + 1]), mod_polys[2 * i]).second; coef_polys[2 * i + 1] = DivMod(Mul(coef_polys[i], mod_polys[2 * i]), mod_polys[2 * i + 1]).second; } for (int i = n; i < 2 * n; ++i) { assert(SZ(coef_polys[i]) == 1); const int coef = coef_polys[i][0]; coef_polys[i][0] = div_mod(ys[i - n], coef); } for (int i = n - 1; i > 0; --i) { coef_polys[i] = Add( Mul(coef_polys[2 * i], mod_polys[2 * i + 1]), Mul(coef_polys[2 * i + 1], mod_polys[2 * i])); } return coef_polys[1]; } }; int main() { ios_base::sync_with_stdio(false); const int p = 998244353; FormalSeries<p> fft(1 << 20); int n; cin >> n; vi xs(n), ys(n); for (int &x : xs) { cin >> x; } for (int &x : ys) { cin >> x; } auto res = fft.Interpolation(xs, ys); for (int x : res) { cout << x << " "; } cout << "\n"; }