This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub maspypy/library
#include "setfunc/transposed_sps_composition.hpp"
#include "setfunc/subset_convolution.hpp" #include "poly/convolution_naive.hpp" // for fixed sps s, consider linear map F:a->b = subset-conv(a,s) // given x, calculate transpose(F)(x) template <typename mint, int LIM> vc<mint> transposed_subset_convolution(vc<mint> s, vc<mint> x) { /* sum_{j}x_jb_j = sum_{i subset j}x_ja_is_{j-i} = sum_{i}y_ia_i. y_i = sum_{j supset i}x_js_{j-i} (rev y)_i = sum_{j subset i}(rev x)_js_{i-j} y = rev(conv(rev x), s) */ reverse(all(x)); x = subset_convolution<mint, LIM>(x, s); reverse(all(x)); return x; } // for fixed sps s s.t. s[0] == 0. // consider linear map F:f->t=f(s) for egf f. // given x, calcuate transpose(F)(x) // equivalent: calculate sum_i x_i(s^k/k!)_i for k=0,1,...,N template <typename mint, int LIM> vc<mint> transposed_sps_composition_egf(vc<mint>& s, vc<mint> x) { const int N = topbit(len(s)); assert(len(s) == (1 << N) && len(x) == (1 << N) && s[0] == mint(0)); vc<mint> y(N + 1); y[0] = x[0]; auto& dp = x; FOR(i, N) { vc<mint> newdp(1 << (N - 1 - i)); FOR(j, N - i) { vc<mint> a = {s.begin() + (1 << j), s.begin() + (2 << j)}; vc<mint> b = {dp.begin() + (1 << j), dp.begin() + (2 << j)}; b = transposed_subset_convolution<mint, LIM>(a, b); FOR(k, len(b)) newdp[k] += b[k]; } swap(dp, newdp); y[1 + i] = dp[0]; } return y; } // for fixed sps s s.t. s[0] == 0. // consider linear map F:f->t=f(s) for polynomial f. // given x, calcuate transpose(F)(x) // equivalent: calculate sum_i x_i(s^k/k!)_i for k=0,1,...,M-1 template <typename mint, int LIM> vc<mint> transposed_sps_composition_poly(vc<mint> s, vc<mint> x, int M) { const int N = topbit(len(s)); assert(len(s) == (1 << N) && len(x) == (1 << N)); mint c = s[0]; s[0] -= c; x = transposed_sps_composition_egf<mint, LIM>(s, x); vc<mint> g(M); mint pow = 1; FOR(i, M) { g[i] = pow * fact_inv<mint>(i), pow *= c; } x = convolution_naive<mint>(x, g); x.resize(M); FOR(i, M) x[i] *= fact<mint>(i); return x; }
#line 2 "setfunc/subset_convolution.hpp" #line 2 "setfunc/ranked_zeta.hpp" template <typename T, int LIM> vc<array<T, LIM + 1>> ranked_zeta(const vc<T>& f) { int n = topbit(len(f)); assert(n <= LIM); assert(len(f) == 1 << n); vc<array<T, LIM + 1>> Rf(1 << n); for (int s = 0; s < (1 << n); ++s) Rf[s][popcnt(s)] = f[s]; for (int i = 0; i < n; ++i) { int w = 1 << i; for (int p = 0; p < (1 << n); p += 2 * w) { for (int s = p; s < p + w; ++s) { int t = s | 1 << i; for (int d = 0; d <= n; ++d) Rf[t][d] += Rf[s][d]; } } } return Rf; } template <typename T, int LIM> vc<T> ranked_mobius(vc<array<T, LIM + 1>>& Rf) { int n = topbit(len(Rf)); assert(len(Rf) == 1 << n); for (int i = 0; i < n; ++i) { int w = 1 << i; for (int p = 0; p < (1 << n); p += 2 * w) { for (int s = p; s < p + w; ++s) { int t = s | 1 << i; for (int d = 0; d <= n; ++d) Rf[t][d] -= Rf[s][d]; } } } vc<T> f(1 << n); for (int s = 0; s < (1 << n); ++s) f[s] = Rf[s][popcnt(s)]; return f; } #line 4 "setfunc/subset_convolution.hpp" template <typename T, int LIM = 20> vc<T> subset_convolution_square(const vc<T>& A) { auto RA = ranked_zeta<T, LIM>(A); int n = topbit(len(RA)); FOR(s, len(RA)) { auto& f = RA[s]; FOR_R(d, n + 1) { T x = 0; FOR(i, d + 1) x += f[i] * f[d - i]; f[d] = x; } } return ranked_mobius<T, LIM>(RA); } template <typename T, int LIM = 20> vc<T> subset_convolution(const vc<T>& A, const vc<T>& B) { if (A == B) return subset_convolution_square(A); auto RA = ranked_zeta<T, LIM>(A); auto RB = ranked_zeta<T, LIM>(B); int n = topbit(len(RA)); FOR(s, len(RA)) { auto &f = RA[s], &g = RB[s]; FOR_R(d, n + 1) { T x = 0; FOR(i, d + 1) x += f[i] * g[d - i]; f[d] = x; } } return ranked_mobius<T, LIM>(RA); } #line 2 "poly/convolution_naive.hpp" template <class T, typename enable_if<!has_mod<T>::value>::type* = nullptr> vc<T> convolution_naive(const vc<T>& a, const vc<T>& b) { int n = int(a.size()), m = int(b.size()); if (n > m) return convolution_naive<T>(b, a); if (n == 0) return {}; vector<T> ans(n + m - 1); FOR(i, n) FOR(j, m) ans[i + j] += a[i] * b[j]; return ans; } template <class T, typename enable_if<has_mod<T>::value>::type* = nullptr> vc<T> convolution_naive(const vc<T>& a, const vc<T>& b) { int n = int(a.size()), m = int(b.size()); if (n > m) return convolution_naive<T>(b, a); if (n == 0) return {}; vc<T> ans(n + m - 1); if (n <= 16 && (T::get_mod() < (1 << 30))) { for (int k = 0; k < n + m - 1; ++k) { int s = max(0, k - m + 1); int t = min(n, k + 1); u64 sm = 0; for (int i = s; i < t; ++i) { sm += u64(a[i].val) * (b[k - i].val); } ans[k] = sm; } } else { for (int k = 0; k < n + m - 1; ++k) { int s = max(0, k - m + 1); int t = min(n, k + 1); u128 sm = 0; for (int i = s; i < t; ++i) { sm += u64(a[i].val) * (b[k - i].val); } ans[k] = T::raw(sm % T::get_mod()); } } return ans; } #line 3 "setfunc/transposed_sps_composition.hpp" // for fixed sps s, consider linear map F:a->b = subset-conv(a,s) // given x, calculate transpose(F)(x) template <typename mint, int LIM> vc<mint> transposed_subset_convolution(vc<mint> s, vc<mint> x) { /* sum_{j}x_jb_j = sum_{i subset j}x_ja_is_{j-i} = sum_{i}y_ia_i. y_i = sum_{j supset i}x_js_{j-i} (rev y)_i = sum_{j subset i}(rev x)_js_{i-j} y = rev(conv(rev x), s) */ reverse(all(x)); x = subset_convolution<mint, LIM>(x, s); reverse(all(x)); return x; } // for fixed sps s s.t. s[0] == 0. // consider linear map F:f->t=f(s) for egf f. // given x, calcuate transpose(F)(x) // equivalent: calculate sum_i x_i(s^k/k!)_i for k=0,1,...,N template <typename mint, int LIM> vc<mint> transposed_sps_composition_egf(vc<mint>& s, vc<mint> x) { const int N = topbit(len(s)); assert(len(s) == (1 << N) && len(x) == (1 << N) && s[0] == mint(0)); vc<mint> y(N + 1); y[0] = x[0]; auto& dp = x; FOR(i, N) { vc<mint> newdp(1 << (N - 1 - i)); FOR(j, N - i) { vc<mint> a = {s.begin() + (1 << j), s.begin() + (2 << j)}; vc<mint> b = {dp.begin() + (1 << j), dp.begin() + (2 << j)}; b = transposed_subset_convolution<mint, LIM>(a, b); FOR(k, len(b)) newdp[k] += b[k]; } swap(dp, newdp); y[1 + i] = dp[0]; } return y; } // for fixed sps s s.t. s[0] == 0. // consider linear map F:f->t=f(s) for polynomial f. // given x, calcuate transpose(F)(x) // equivalent: calculate sum_i x_i(s^k/k!)_i for k=0,1,...,M-1 template <typename mint, int LIM> vc<mint> transposed_sps_composition_poly(vc<mint> s, vc<mint> x, int M) { const int N = topbit(len(s)); assert(len(s) == (1 << N) && len(x) == (1 << N)); mint c = s[0]; s[0] -= c; x = transposed_sps_composition_egf<mint, LIM>(s, x); vc<mint> g(M); mint pow = 1; FOR(i, M) { g[i] = pow * fact_inv<mint>(i), pow *= c; } x = convolution_naive<mint>(x, g); x.resize(M); FOR(i, M) x[i] *= fact<mint>(i); return x; }