#line 1 "test/formal_power_series/composition_of_formal_power_series_large.0.test.cpp" // competitive-verifier: PROBLEM https://judge.yosupo.jp/problem/composition_of_formal_power_series_large #line 2 "fps_composition.hpp" #line 2 "binomial.hpp" #include #include template class Binomial { std::vector factorial_, invfactorial_; Binomial() : factorial_{Tp(1)}, invfactorial_{Tp(1)} {} void preprocess(int n) { if (const int nn = factorial_.size(); nn < n) { int k = nn; while (k < n) k *= 2; k = std::min(k, Tp::mod()); factorial_.resize(k); invfactorial_.resize(k); for (int i = nn; i < k; ++i) factorial_[i] = factorial_[i - 1] * i; invfactorial_.back() = factorial_.back().inv(); for (int i = k - 2; i >= nn; --i) invfactorial_[i] = invfactorial_[i + 1] * (i + 1); } } public: static const Binomial &get(int n) { static Binomial bin; bin.preprocess(n); return bin; } Tp binom(int n, int m) const { return n < m ? Tp() : factorial_[n] * invfactorial_[m] * invfactorial_[n - m]; } Tp inv(int n) const { return factorial_[n - 1] * invfactorial_[n]; } Tp factorial(int n) const { return factorial_[n]; } Tp inv_factorial(int n) const { return invfactorial_[n]; } }; #line 2 "fft.hpp" #line 4 "fft.hpp" #include #include #include #line 8 "fft.hpp" template class FftInfo { static Tp least_quadratic_nonresidue() { for (int i = 2;; ++i) if (Tp(i).pow((Tp::mod() - 1) / 2) == -1) return Tp(i); } const int ordlog2_; const Tp zeta_; const Tp invzeta_; const Tp imag_; const Tp invimag_; mutable std::vector root_; mutable std::vector invroot_; FftInfo() : ordlog2_(__builtin_ctzll(Tp::mod() - 1)), zeta_(least_quadratic_nonresidue().pow((Tp::mod() - 1) >> ordlog2_)), invzeta_(zeta_.inv()), imag_(zeta_.pow(1LL << (ordlog2_ - 2))), invimag_(-imag_), root_{Tp(1), imag_}, invroot_{Tp(1), invimag_} {} public: static const FftInfo &get() { static FftInfo info; return info; } Tp imag() const { return imag_; } Tp inv_imag() const { return invimag_; } Tp zeta() const { return zeta_; } Tp inv_zeta() const { return invzeta_; } const std::vector &root(int n) const { // [0, n) assert((n & (n - 1)) == 0); if (const int s = root_.size(); s < n) { root_.resize(n); for (int i = __builtin_ctz(s); (1 << i) < n; ++i) { const int j = 1 << i; root_[j] = zeta_.pow(1LL << (ordlog2_ - i - 2)); for (int k = j + 1; k < j * 2; ++k) root_[k] = root_[k - j] * root_[j]; } } return root_; } const std::vector &inv_root(int n) const { // [0, n) assert((n & (n - 1)) == 0); if (const int s = invroot_.size(); s < n) { invroot_.resize(n); for (int i = __builtin_ctz(s); (1 << i) < n; ++i) { const int j = 1 << i; invroot_[j] = invzeta_.pow(1LL << (ordlog2_ - i - 2)); for (int k = j + 1; k < j * 2; ++k) invroot_[k] = invroot_[k - j] * invroot_[j]; } } return invroot_; } }; inline int fft_len(int n) { --n; n |= n >> 1, n |= n >> 2, n |= n >> 4, n |= n >> 8; return (n | n >> 16) + 1; } namespace detail { template inline void butterfly_n(Iterator a, int n, const std::vector::value_type> &root) { assert(n > 0); assert((n & (n - 1)) == 0); const int bn = __builtin_ctz(n); if (bn & 1) { for (int i = 0; i < n / 2; ++i) { const auto a0 = a[i], a1 = a[i + n / 2]; a[i] = a0 + a1, a[i + n / 2] = a0 - a1; } } for (int i = n >> (bn & 1); i >= 4; i /= 4) { const int i4 = i / 4; for (int k = 0; k < i4; ++k) { const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1]; const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3]; const auto a02p = a0 + a2, a02m = a0 - a2; const auto a13p = a1 + a3, a13m = (a1 - a3) * root[1]; a[k + i4 * 0] = a02p + a13p, a[k + i4 * 1] = a02p - a13p; a[k + i4 * 2] = a02m + a13m, a[k + i4 * 3] = a02m - a13m; } for (int j = i, m = 2; j < n; j += i, m += 2) { const auto r = root[m], r2 = r * r, r3 = r2 * r; for (int k = j; k < j + i4; ++k) { const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1] * r; const auto a2 = a[k + i4 * 2] * r2, a3 = a[k + i4 * 3] * r3; const auto a02p = a0 + a2, a02m = a0 - a2; const auto a13p = a1 + a3, a13m = (a1 - a3) * root[1]; a[k + i4 * 0] = a02p + a13p, a[k + i4 * 1] = a02p - a13p; a[k + i4 * 2] = a02m + a13m, a[k + i4 * 3] = a02m - a13m; } } } } template inline void inv_butterfly_n(Iterator a, int n, const std::vector::value_type> &root) { assert(n > 0); assert((n & (n - 1)) == 0); const int bn = __builtin_ctz(n); for (int i = 4; i <= (n >> (bn & 1)); i *= 4) { const int i4 = i / 4; for (int k = 0; k < i4; ++k) { const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1]; const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3]; const auto a01p = a0 + a1, a01m = a0 - a1; const auto a23p = a2 + a3, a23m = (a2 - a3) * root[1]; a[k + i4 * 0] = a01p + a23p, a[k + i4 * 1] = a01m + a23m; a[k + i4 * 2] = a01p - a23p, a[k + i4 * 3] = a01m - a23m; } for (int j = i, m = 2; j < n; j += i, m += 2) { const auto r = root[m], r2 = r * r, r3 = r2 * r; for (int k = j; k < j + i4; ++k) { const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1]; const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3]; const auto a01p = a0 + a1, a01m = a0 - a1; const auto a23p = a2 + a3, a23m = (a2 - a3) * root[1]; a[k + i4 * 0] = a01p + a23p, a[k + i4 * 1] = (a01m + a23m) * r; a[k + i4 * 2] = (a01p - a23p) * r2, a[k + i4 * 3] = (a01m - a23m) * r3; } } } if (bn & 1) { for (int i = 0; i < n / 2; ++i) { const auto a0 = a[i], a1 = a[i + n / 2]; a[i] = a0 + a1, a[i + n / 2] = a0 - a1; } } } } // namespace detail // FFT_n: A(x) |-> bit-reversed order of [A(1), A(zeta_n), ..., A(zeta_n^(n-1))] template inline void fft_n(Iterator a, int n) { using Tp = typename std::iterator_traits::value_type; detail::butterfly_n(a, n, FftInfo::get().root(n / 2)); } template inline void fft(std::vector &a) { fft_n(a.begin(), a.size()); } // IFFT_n: bit-reversed order of [A(1), A(zeta_n), ..., A(zeta_n^(n-1))] |-> A(x) template inline void inv_fft_n(Iterator a, int n) { using Tp = typename std::iterator_traits::value_type; detail::inv_butterfly_n(a, n, FftInfo::get().inv_root(n / 2)); const Tp iv = Tp::mod() - (Tp::mod() - 1) / n; for (int i = 0; i < n; ++i) a[i] *= iv; } template inline void inv_fft(std::vector &a) { inv_fft_n(a.begin(), a.size()); } // IFFT_n^T: A(x) |-> 1/n FFT_n((x^n A(x^(-1))) mod (x^n - 1)) template inline void transposed_inv_fft_n(Iterator a, int n) { using Tp = typename std::iterator_traits::value_type; const Tp iv = Tp::mod() - (Tp::mod() - 1) / n; for (int i = 0; i < n; ++i) a[i] *= iv; detail::butterfly_n(a, n, FftInfo::get().inv_root(n / 2)); } template inline void transposed_inv_fft(std::vector &a) { transposed_inv_fft_n(a.begin(), a.size()); } // FFT_n^T : FFT_n((x^n A(x^(-1))) mod (x^n - 1)) |-> n A(x) template inline void transposed_fft_n(Iterator a, int n) { using Tp = typename std::iterator_traits::value_type; detail::inv_butterfly_n(a, n, FftInfo::get().root(n / 2)); } template inline void transposed_fft(std::vector &a) { transposed_fft_n(a.begin(), a.size()); } template inline std::vector convolution_fft(std::vector a, std::vector b) { if (a.empty() || b.empty()) return {}; const int n = a.size(); const int m = b.size(); const int len = fft_len(n + m - 1); a.resize(len); b.resize(len); fft(a); fft(b); for (int i = 0; i < len; ++i) a[i] *= b[i]; inv_fft(a); a.resize(n + m - 1); return a; } template inline std::vector square_fft(std::vector a) { if (a.empty()) return {}; const int n = a.size(); const int len = fft_len(n * 2 - 1); a.resize(len); fft(a); for (int i = 0; i < len; ++i) a[i] *= a[i]; inv_fft(a); a.resize(n * 2 - 1); return a; } template inline std::vector convolution_naive(const std::vector &a, const std::vector &b) { if (a.empty() || b.empty()) return {}; const int n = a.size(); const int m = b.size(); std::vector res(n + m - 1); for (int i = 0; i < n; ++i) for (int j = 0; j < m; ++j) res[i + j] += a[i] * b[j]; return res; } template inline std::vector convolution(const std::vector &a, const std::vector &b) { if (std::min(a.size(), b.size()) < 60) return convolution_naive(a, b); if (std::addressof(a) == std::addressof(b)) return square_fft(a); return convolution_fft(a, b); } #line 2 "fps_basic.hpp" #line 2 "semi_relaxed_conv.hpp" #line 5 "semi_relaxed_conv.hpp" #include #include #line 8 "semi_relaxed_conv.hpp" template inline std::enable_if_t &>, std::vector> semi_relaxed_convolution_naive(const std::vector &A, Closure gen, int n) { std::vector B(n), AB(n); for (int i = 0; i < n; ++i) { for (int j = std::max(0, i - (int)A.size() + 1); j < i; ++j) AB[i] += A[i - j] * B[j]; B[i] = gen(i, AB); if (!A.empty()) AB[i] += A[0] * B[i]; } return B; } // returns coefficients generated by closure // closure: gen(index, current_product) template inline std::enable_if_t &>, std::vector> semi_relaxed_convolution(const std::vector &A, Closure gen, int n) { if (A.size() < 60) return semi_relaxed_convolution_naive(A, gen, n); enum { BaseCaseSize = 32 }; static_assert((BaseCaseSize & (BaseCaseSize - 1)) == 0); static const int Block[] = {16, 16, 16, 16, 16}; static const int BlockSize[] = { BaseCaseSize, BaseCaseSize * Block[0], BaseCaseSize * Block[0] * Block[1], BaseCaseSize * Block[0] * Block[1] * Block[2], BaseCaseSize * Block[0] * Block[1] * Block[2] * Block[3], BaseCaseSize * Block[0] * Block[1] * Block[2] * Block[3] * Block[4], }; // returns (which_block, level) auto blockinfo = [](int ind) { int i = ind / BaseCaseSize, lv = 0; while ((i & (Block[lv] - 1)) == 0) i /= Block[lv++]; return std::make_pair(i & (Block[lv] - 1), lv); }; std::vector B(n), AB(n); std::vector>> dftA, dftB; for (int i = 0; i < n; ++i) { const int s = i & (BaseCaseSize - 1); // block contribution if (i >= BaseCaseSize && s == 0) { const auto [j, lv] = blockinfo(i); const int blocksize = BlockSize[lv]; if (blocksize * j == i) { if ((int)dftA.size() == lv) { dftA.emplace_back(); dftB.emplace_back(Block[lv] - 1); } if ((j - 1) * blocksize < (int)A.size()) { dftA[lv] .emplace_back(A.begin() + (j - 1) * blocksize, A.begin() + std::min((j + 1) * blocksize, A.size())) .resize(blocksize * 2); fft(dftA[lv][j - 1]); } } if (!dftA[lv].empty()) { dftB[lv][j - 1].resize(blocksize * 2); std::copy_n(B.begin() + (i - blocksize), blocksize, dftB[lv][j - 1].begin()); std::fill_n(dftB[lv][j - 1].begin() + blocksize, blocksize, Tp(0)); fft(dftB[lv][j - 1]); // middle product std::vector mp(blocksize * 2); for (int k = 0; k < std::min(j, dftA[lv].size()); ++k) for (int l = 0; l < blocksize * 2; ++l) mp[l] += dftA[lv][k][l] * dftB[lv][j - 1 - k][l]; inv_fft(mp); for (int k = 0; k < blocksize && i + k < n; ++k) AB[i + k] += mp[k + blocksize]; } } // basecase contribution for (int j = std::max(i - s, i - (int)A.size() + 1); j < i; ++j) AB[i] += A[i - j] * B[j]; B[i] = gen(i, AB); if (!A.empty()) AB[i] += A[0] * B[i]; } return B; } #line 8 "fps_basic.hpp" template inline int order(const std::vector &a) { for (int i = 0; i < (int)a.size(); ++i) if (a[i] != 0) return i; return -1; } template inline std::vector fps_inv(const std::vector &a, int n) { assert(order(a) == 0); if (n <= 0) return {}; if (std::min(a.size(), n) < 60) return semi_relaxed_convolution( a, [v = a[0].inv()](int n, auto &&c) { return n == 0 ? v : -c[n] * v; }, n); enum { Threshold = 32 }; const int len = fft_len(n); std::vector invA, shopA(len), shopB(len); invA = semi_relaxed_convolution( a, [v = a[0].inv()](int n, auto &&c) { return n == 0 ? v : -c[n] * v; }, Threshold); invA.resize(len); for (int i = Threshold * 2; i <= len; i *= 2) { std::fill(std::copy_n(a.begin(), std::min(a.size(), i), shopA.begin()), shopA.begin() + i, Tp(0)); std::copy_n(invA.begin(), i, shopB.begin()); fft_n(shopA.begin(), i); fft_n(shopB.begin(), i); for (int j = 0; j < i; ++j) shopA[j] *= shopB[j]; inv_fft_n(shopA.begin(), i); std::fill_n(shopA.begin(), i / 2, Tp(0)); fft_n(shopA.begin(), i); for (int j = 0; j < i; ++j) shopA[j] *= shopB[j]; inv_fft_n(shopA.begin(), i); for (int j = i / 2; j < i; ++j) invA[j] = -shopA[j]; } invA.resize(n); return invA; } template inline std::vector fps_div(const std::vector &a, const std::vector &b, int n) { assert(order(b) == 0); if (n <= 0) return {}; return semi_relaxed_convolution( b, [&, v = b[0].inv()](int n, auto &&c) { if (n < (int)a.size()) return (a[n] - c[n]) * v; return -c[n] * v; }, n); } template inline std::vector deriv(const std::vector &a) { const int n = (int)a.size() - 1; if (n <= 0) return {}; std::vector res(n); for (int i = 1; i <= n; ++i) res[i - 1] = a[i] * i; return res; } template inline std::vector integr(const std::vector &a, Tp c = {}) { const int n = a.size() + 1; auto &&bin = Binomial::get(n); std::vector res(n); res[0] = c; for (int i = 1; i < n; ++i) res[i] = a[i - 1] * bin.inv(i); return res; } template inline std::vector fps_log(const std::vector &a, int n) { return integr(fps_div(deriv(a), a, n - 1)); } template inline std::vector fps_exp(const std::vector &a, int n) { if (n <= 0) return {}; assert(a.empty() || a[0] == 0); return semi_relaxed_convolution( deriv(a), [bin = Binomial::get(n)](int n, auto &&c) { return n == 0 ? Tp(1) : c[n - 1] * bin.inv(n); }, n); } template inline std::vector fps_pow(std::vector a, long long e, int n) { if (n <= 0) return {}; if (e == 0) { std::vector res(n); res[0] = 1; return res; } const int o = order(a); if (o < 0 || o > n / e || (o == n / e && n % e == 0)) return std::vector(n); if (o != 0) a.erase(a.begin(), a.begin() + o); const Tp ia0 = a[0].inv(); const Tp a0e = a[0].pow(e); const Tp me = e; for (int i = 0; i < (int)a.size(); ++i) a[i] *= ia0; a = fps_log(a, n - o * e); for (int i = 0; i < (int)a.size(); ++i) a[i] *= me; a = fps_exp(a, n - o * e); for (int i = 0; i < (int)a.size(); ++i) a[i] *= a0e; a.insert(a.begin(), o * e, Tp(0)); return a; } #line 10 "fps_composition.hpp" // returns f(g) mod x^n // see: // [1]: Yasunori Kinoshita, Baitian Li. Power Series Composition in Near-Linear Time. // https://arxiv.org/abs/2404.05177 template inline std::vector composition(const std::vector &f, const std::vector &g, int n) { if (n <= 0) return {}; if (g.empty()) { std::vector res(n); if (!f.empty()) res[0] = f[0]; return res; } // [y^(-1)] (f(y) / (-g(x) + y)) mod x^n in R[x]((y^(-1))) auto rec = [g0 = g[0]](auto &&rec, const std::vector &P, const std::vector &Q, int d, int n) { if (n == 1) { std::vector invQ(d + 1); auto &&bin = Binomial::get(d * 2); Tp gg = 1; for (int i = 0; i <= d; ++i) invQ[d - i] = bin.binom(d + i - 1, d - 1) * gg, gg *= g0; // invQ[i] = [y^(-2d + i)]Q^(-1) // P[0,d-1] * invQ[-2d,-d] => [0,d-1] * [0,d] // take [-d,-1] => take [d,2d-1] auto PinvQ = convolution(P, invQ); PinvQ.erase(PinvQ.begin(), PinvQ.begin() + d); PinvQ.resize(d); return PinvQ; } std::vector dftQ(d * n * 4); for (int i = 0; i < d; ++i) for (int j = 0; j < n; ++j) dftQ[i * (n * 2) + j] = Q[i * n + j]; dftQ[d * n * 2] = 1; fft(dftQ); std::vector V(d * n * 2); for (int i = 0; i < d * n * 4; i += 2) V[i / 2] = dftQ[i] * dftQ[i + 1]; inv_fft(V); V[0] -= 1; for (int i = 1; i < d * 2; ++i) for (int j = 0; j < n / 2; ++j) V[i * (n / 2) + j] = V[i * n + j]; V.resize(d * n); const auto T = rec(rec, P, std::move(V), d * 2, n / 2); std::vector dftT(d * n * 2); for (int i = 0; i < d * 2; ++i) for (int j = 0; j < n / 2; ++j) dftT[i * n + j] = T[i * (n / 2) + j]; fft(dftT); std::vector U(d * n * 4); for (int i = 0; i < d * n * 4; i += 2) { U[i] = dftT[i / 2] * dftQ[i + 1]; U[i + 1] = dftT[i / 2] * dftQ[i]; } inv_fft(U); // [-2d,d-1] => [0,3d-1] // take [-d,-1] => take [d,2d-1] for (int i = 0; i < d; ++i) for (int j = 0; j < n; ++j) U[i * n + j] = U[(i + d) * (n * 2) + j]; U.resize(d * n); return U; }; const int k = fft_len(std::max(n, f.size())); std::vector Q(k); for (int i = 0; i < std::min(k, g.size()); ++i) Q[i] = -g[i]; auto res = rec(rec, f, Q, 1, k); res.resize(n); return res; } // returns [x^k]gf^0, [x^k]gf, ..., [x^k]gf^(n-1) // see: // [1]: noshi91. FPS の合成と逆関数、冪乗の係数列挙 Θ(n (log(n))^2) // https://noshi91.hatenablog.com/entry/2024/03/16/224034 template inline std::vector enum_kth_term_of_power(const std::vector &f, const std::vector &g, int k, int n) { if (k < 0 || n <= 0) return {}; if (f.empty()) { std::vector res(n); if (k < (int)g.size()) res[0] = g[k]; return res; } // [x^k] (g(x) / (-f(x) + y)) in R[x]((y^(-1))) std::vector P(g), Q(k + 1); P.resize(k + 1); for (int i = 0; i < std::min(k + 1, f.size()); ++i) Q[i] = -f[i]; int d = 1; for (; k; d *= 2, k /= 2) { const int len = fft_len((d * 2) * ((k + 1) * 2) - 1); std::vector dftP(len), dftQ(len); for (int i = 0; i < d; ++i) for (int j = 0; j <= k; ++j) { dftP[i * ((k + 1) * 2) + j] = P[i * (k + 1) + j]; dftQ[i * ((k + 1) * 2) + j] = Q[i * (k + 1) + j]; } dftQ[d * (k + 1) * 2] = 1; fft(dftP); fft(dftQ); P.resize(len / 2); Q.resize(len / 2); if (k & 1) { auto &&root = FftInfo::get().inv_root(len / 2); for (int i = 0; i < len; i += 2) { P[i / 2] = (dftP[i] * dftQ[i + 1] - dftP[i + 1] * dftQ[i]).div_by_2() * root[i / 2]; Q[i / 2] = dftQ[i] * dftQ[i + 1]; } } else { for (int i = 0; i < len; i += 2) { P[i / 2] = (dftP[i] * dftQ[i + 1] + dftP[i + 1] * dftQ[i]).div_by_2(); Q[i / 2] = dftQ[i] * dftQ[i + 1]; } } inv_fft(P); inv_fft(Q); if (d * (k + 1) * 4 >= len) Q[(d * (k + 1) * 4) % len] -= 1; for (int i = 1; i < d * 2; ++i) for (int j = 0; j <= k / 2; ++j) { P[i * (k / 2 + 1) + j] = P[i * (k + 1) + j]; Q[i * (k / 2 + 1) + j] = Q[i * (k + 1) + j]; } P.resize(d * 2 * (k / 2 + 1)); Q.resize(d * 2 * (k / 2 + 1)); } std::vector invQ(n + 1); auto &&bin = Binomial::get(d + n); Tp ff = 1; for (int i = 0; i <= n; ++i) invQ[n - i] = bin.binom(d + i - 1, d - 1) * ff, ff *= f[0]; // invQ[i] = [y^(-2d + i)]Q^(-1) // P[0,d-1] * invQ[-(d+n),-d] => [0,d-1] * [0,n] auto PinvQ = convolution(P, invQ); // take [-n,-1] => take [d,d+n-1] PinvQ.erase(PinvQ.begin(), PinvQ.begin() + d); PinvQ.resize(n); // output => [-1,-n] reverse // before I just reverse it and mistaken something. std::reverse(PinvQ.begin(), PinvQ.end()); return PinvQ; } // returns g s.t. f(g) = g(f) = x mod x^n template inline std::vector reversion(std::vector f, int n) { if (n <= 0 || f.size() < 2) return {}; assert(order(f) == 1); const auto if1 = f[1].inv(); if (n == 1) return {Tp(0)}; f.resize(n); Tp ff = 1; for (int i = 1; i < n; ++i) f[i] *= ff *= if1; auto a = enum_kth_term_of_power(f, {Tp(1)}, n - 1, n); auto &&bin = Binomial::get(n); for (int i = 1; i < n; ++i) a[i] *= (n - 1) * bin.inv(i); auto b = fps_pow(std::vector(a.rbegin(), a.rend() - 1), Tp(1 - n).inv().val(), n - 1); for (int i = 0; i < n - 1; ++i) b[i] *= if1; b.insert(b.begin(), Tp(0)); return b; } #line 2 "modint.hpp" #include #line 5 "modint.hpp" // clang-format off template class ModInt { static_assert((Mod >> 31) == 0, "`Mod` must less than 2^(31)"); template static std::enable_if_t, unsigned> safe_mod(Int v) { using D = std::common_type_t; return (v %= (int)Mod) < 0 ? (D)(v + (int)Mod) : (D)v; } struct PrivateConstructor {} static inline private_constructor{}; ModInt(PrivateConstructor, unsigned v) : v_(v) {} unsigned v_; public: static unsigned mod() { return Mod; } static ModInt from_raw(unsigned v) { return ModInt(private_constructor, v); } static ModInt zero() { return from_raw(0); } static ModInt one() { return from_raw(1); } bool is_zero() const { return v_ == 0; } bool is_one() const { return v_ == 1; } ModInt() : v_() {} template, int> = 0> ModInt(Int v) : v_(safe_mod(v)) {} template, int> = 0> ModInt(Int v) : v_(v % Mod) {} unsigned val() const { return v_; } ModInt operator-() const { return from_raw(v_ == 0 ? v_ : Mod - v_); } ModInt pow(long long e) const { if (e < 0) return inv().pow(-e); for (ModInt x(*this), res(from_raw(1));; x *= x) { if (e & 1) res *= x; if ((e >>= 1) == 0) return res; }} ModInt inv() const { int x1 = 1, x3 = 0, a = val(), b = Mod; while (b) { const int q = a / b, x1_old = x1, a_old = a; x1 = x3, x3 = x1_old - x3 * q, a = b, b = a_old - b * q; } return from_raw(x1 < 0 ? x1 + (int)Mod : x1); } template std::enable_if_t div_by_2() const { if (v_ & 1) return from_raw((v_ + Mod) >> 1); return from_raw(v_ >> 1); } ModInt &operator+=(const ModInt &a) { if ((v_ += a.v_) >= Mod) v_ -= Mod; return *this; } ModInt &operator-=(const ModInt &a) { if ((v_ += Mod - a.v_) >= Mod) v_ -= Mod; return *this; } ModInt &operator*=(const ModInt &a) { v_ = (unsigned long long)v_ * a.v_ % Mod; return *this; } ModInt &operator/=(const ModInt &a) { return *this *= a.inv(); } ModInt &operator++() { return *this += one(); } ModInt operator++(int) { ModInt o(*this); *this += one(); return o; } ModInt &operator--() { return *this -= one(); } ModInt operator--(int) { ModInt o(*this); *this -= one(); return o; } friend ModInt operator+(const ModInt &a, const ModInt &b) { return ModInt(a) += b; } friend ModInt operator-(const ModInt &a, const ModInt &b) { return ModInt(a) -= b; } friend ModInt operator*(const ModInt &a, const ModInt &b) { return ModInt(a) *= b; } friend ModInt operator/(const ModInt &a, const ModInt &b) { return ModInt(a) /= b; } friend bool operator==(const ModInt &a, const ModInt &b) { return a.v_ == b.v_; } friend bool operator!=(const ModInt &a, const ModInt &b) { return a.v_ != b.v_; } friend std::istream &operator>>(std::istream &a, ModInt &b) { int v; a >> v; b.v_ = safe_mod(v); return a; } friend std::ostream &operator<<(std::ostream &a, const ModInt &b) { return a << b.val(); } }; // clang-format on #line 7 "test/formal_power_series/composition_of_formal_power_series_large.0.test.cpp" template inline std::vector convolution_naive_trunc(const std::vector &a, const std::vector &b, int k) { if (a.empty() || b.empty()) return std::vector(k); const int n = a.size(); const int m = b.size(); std::vector res(k); for (int i = 0; i < std::min(n, k); ++i) for (int j = 0; j < std::min(m, k) && i + j < k; ++j) res[i + j] += a[i] * b[j]; return res; } template inline std::vector convolution_trunc(const std::vector &a, const std::vector &b, int k) { if (k < 60) return convolution_naive_trunc(a, b, k); auto ab = convolution(std::vector(a.begin(), std::min(a.end(), a.begin() + k)), std::vector(b.begin(), std::min(b.end(), b.begin() + k))); ab.resize(k); return ab; } // returns V(z) s.t. V(U(z)) = uV(z) // We ensure Schroeder function of U exists. template std::vector schroeder_function(const std::vector &U, int n) { // see: TAOCP Vol 2 Chap 4.7. // returns V(z) s.t. V(U(z)) = W(z)V(z) + S(z) + O(z^n) auto brent_traub = [](auto &&brent_traub, const std::vector &U, const std::vector &W, const std::vector &S, int n) -> std::vector { assert((n & (n - 1)) == 0); if (n == 1) { auto S0 = S.empty() ? Tp() : S[0]; auto W0 = W.empty() ? Tp() : W[0]; if (S0 == 0 && W0 == 1) return {Tp(1)}; return {S0 / (1 - W0)}; } auto V = brent_traub(brent_traub, U, W, S, n / 2); auto VU = composition(V, U, n); auto WV = convolution_trunc(W, V, n); std::vector R(n / 2); for (int i = n / 2; i < n; ++i) R[i - n / 2] = WV[i] + (i < (int)S.size() ? S[i] : Tp(0)) - VU[i]; std::vector U_hat = U; U_hat.erase(U_hat.begin()); U_hat = fps_pow(fps_inv(U_hat, n / 2), n / 2, n / 2); auto W_hat = convolution_trunc(W, U_hat, n / 2); auto S_hat = convolution_trunc(R, U_hat, n / 2); auto V_hat = brent_traub(brent_traub, U, W_hat, S_hat, n / 2); V.insert(V.end(), V_hat.begin(), V_hat.end()); return V; }; auto V = brent_traub(brent_traub, U, {U.at(1)}, {}, fft_len(n)); V.resize(n); return V; } int main() { std::ios::sync_with_stdio(false); std::cin.tie(nullptr); using mint = ModInt<998244353>; int n; std::cin >> n; std::vector f(n); for (int i = 0; i < n; ++i) std::cin >> f[i]; const auto g = schroeder_function(f, n); const auto h = reversion(g, n); for (int i = 0; i < n; ++i) { if (i) std::cout << ' '; std::cout << g[i]; } std::cout << '\n'; for (int i = 0; i < n; ++i) { if (i) std::cout << ' '; std::cout << h[i]; } std::cout << '\n'; return 0; }