This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub maspypy/library
#include "setfunc/sps_log.hpp"
#include "mod/modint.hpp" #include "setfunc/sps_composition.hpp" // exp の逆手順で計算する template <typename mint, int LIM> vc<mint> sps_log(vc<mint>& dp) { const int N = topbit(len(dp)); assert(len(dp) == (1 << N) && dp[0] == mint(1)); vc<mint> s(1 << N); FOR_R(i, N) { vc<mint> a = {dp.begin() + (1 << i), dp.begin() + (2 << i)}; vc<mint> b = {dp.begin(), dp.begin() + (1 << i)}; auto RA = ranked_zeta<mint, LIM>(a); auto RB = ranked_zeta<mint, LIM>(b); FOR(s, 1 << i) { auto &f = RA[s], &g = RB[s]; // assert(g[0] == mint(1)); FOR(d, i + 1) { FOR(i, d) f[d] -= f[i] * g[d - i]; } } a = ranked_mobius<mint, LIM>(RA); copy(all(a), s.begin() + (1 << i)); } return s; }
#line 2 "mod/modint_common.hpp" struct has_mod_impl { template <class T> static auto check(T &&x) -> decltype(x.get_mod(), std::true_type{}); template <class T> static auto check(...) -> std::false_type; }; template <class T> class has_mod : public decltype(has_mod_impl::check<T>(std::declval<T>())) {}; template <typename mint> mint inv(int n) { static const int mod = mint::get_mod(); static vector<mint> dat = {0, 1}; assert(0 <= n); if (n >= mod) n %= mod; while (len(dat) <= n) { int k = len(dat); int q = (mod + k - 1) / k; dat.eb(dat[k * q - mod] * mint::raw(q)); } return dat[n]; } template <typename mint> mint fact(int n) { static const int mod = mint::get_mod(); assert(0 <= n && n < mod); static vector<mint> dat = {1, 1}; while (len(dat) <= n) dat.eb(dat[len(dat) - 1] * mint::raw(len(dat))); return dat[n]; } template <typename mint> mint fact_inv(int n) { static vector<mint> dat = {1, 1}; if (n < 0) return mint(0); while (len(dat) <= n) dat.eb(dat[len(dat) - 1] * inv<mint>(len(dat))); return dat[n]; } template <class mint, class... Ts> mint fact_invs(Ts... xs) { return (mint(1) * ... * fact_inv<mint>(xs)); } template <typename mint, class Head, class... Tail> mint multinomial(Head &&head, Tail &&... tail) { return fact<mint>(head) * fact_invs<mint>(std::forward<Tail>(tail)...); } template <typename mint> mint C_dense(int n, int k) { static vvc<mint> C; static int H = 0, W = 0; auto calc = [&](int i, int j) -> mint { if (i == 0) return (j == 0 ? mint(1) : mint(0)); return C[i - 1][j] + (j ? C[i - 1][j - 1] : 0); }; if (W <= k) { FOR(i, H) { C[i].resize(k + 1); FOR(j, W, k + 1) { C[i][j] = calc(i, j); } } W = k + 1; } if (H <= n) { C.resize(n + 1); FOR(i, H, n + 1) { C[i].resize(W); FOR(j, W) { C[i][j] = calc(i, j); } } H = n + 1; } return C[n][k]; } template <typename mint, bool large = false, bool dense = false> mint C(ll n, ll k) { assert(n >= 0); if (k < 0 || n < k) return 0; if constexpr (dense) return C_dense<mint>(n, k); if constexpr (!large) return multinomial<mint>(n, k, n - k); k = min(k, n - k); mint x(1); FOR(i, k) x *= mint(n - i); return x * fact_inv<mint>(k); } template <typename mint, bool large = false> mint C_inv(ll n, ll k) { assert(n >= 0); assert(0 <= k && k <= n); if (!large) return fact_inv<mint>(n) * fact<mint>(k) * fact<mint>(n - k); return mint(1) / C<mint, 1>(n, k); } // [x^d](1-x)^{-n} template <typename mint, bool large = false, bool dense = false> mint C_negative(ll n, ll d) { assert(n >= 0); if (d < 0) return mint(0); if (n == 0) { return (d == 0 ? mint(1) : mint(0)); } return C<mint, large, dense>(n + d - 1, d); } #line 3 "mod/modint.hpp" template <int mod> struct modint { static constexpr u32 umod = u32(mod); static_assert(umod < u32(1) << 31); u32 val; static modint raw(u32 v) { modint x; x.val = v; return x; } constexpr modint() : val(0) {} constexpr modint(u32 x) : val(x % umod) {} constexpr modint(u64 x) : val(x % umod) {} constexpr modint(u128 x) : val(x % umod) {} constexpr modint(int x) : val((x %= mod) < 0 ? x + mod : x){}; constexpr modint(ll x) : val((x %= mod) < 0 ? x + mod : x){}; constexpr modint(i128 x) : val((x %= mod) < 0 ? x + mod : x){}; bool operator<(const modint &other) const { return val < other.val; } modint &operator+=(const modint &p) { if ((val += p.val) >= umod) val -= umod; return *this; } modint &operator-=(const modint &p) { if ((val += umod - p.val) >= umod) val -= umod; return *this; } modint &operator*=(const modint &p) { val = u64(val) * p.val % umod; return *this; } modint &operator/=(const modint &p) { *this *= p.inverse(); return *this; } modint operator-() const { return modint::raw(val ? mod - val : u32(0)); } modint operator+(const modint &p) const { return modint(*this) += p; } modint operator-(const modint &p) const { return modint(*this) -= p; } modint operator*(const modint &p) const { return modint(*this) *= p; } modint operator/(const modint &p) const { return modint(*this) /= p; } bool operator==(const modint &p) const { return val == p.val; } bool operator!=(const modint &p) const { return val != p.val; } modint inverse() const { int a = val, b = mod, u = 1, v = 0, t; while (b > 0) { t = a / b; swap(a -= t * b, b), swap(u -= t * v, v); } return modint(u); } modint pow(ll n) const { assert(n >= 0); modint ret(1), mul(val); while (n > 0) { if (n & 1) ret *= mul; mul *= mul; n >>= 1; } return ret; } static constexpr int get_mod() { return mod; } // (n, r), r は 1 の 2^n 乗根 static constexpr pair<int, int> ntt_info() { if (mod == 120586241) return {20, 74066978}; if (mod == 167772161) return {25, 17}; if (mod == 469762049) return {26, 30}; if (mod == 754974721) return {24, 362}; if (mod == 880803841) return {23, 211}; if (mod == 943718401) return {22, 663003469}; if (mod == 998244353) return {23, 31}; if (mod == 1004535809) return {21, 582313106}; if (mod == 1012924417) return {21, 368093570}; return {-1, -1}; } static constexpr bool can_ntt() { return ntt_info().fi != -1; } }; #ifdef FASTIO template <int mod> void rd(modint<mod> &x) { fastio::rd(x.val); x.val %= mod; // assert(0 <= x.val && x.val < mod); } template <int mod> void wt(modint<mod> x) { fastio::wt(x.val); } #endif using modint107 = modint<1000000007>; using modint998 = modint<998244353>; #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 2 "setfunc/sps_composition.hpp" // sum_i f_i/i! s^i, s^i is subset-convolution template <typename mint, int LIM> vc<mint> sps_composition_egf(vc<mint>& f, vc<mint>& s) { const int N = topbit(len(s)); assert(len(s) == (1 << N) && s[0] == mint(0)); if (len(f) > N) f.resize(N + 1); int D = len(f) - 1; using ARR = array<mint, LIM + 1>; vvc<ARR> zs(N); FOR(i, N) { zs[i] = ranked_zeta<mint, LIM>({s.begin() + (1 << i), s.begin() + (2 << i)}); } // dp : (d/dt)^df(s) (d=D,D-1,...) vc<mint> dp(1 << (N - D)); dp[0] = f[D]; FOR_R(d, D) { vc<mint> newdp(1 << (N - d)); newdp[0] = f[d]; vc<ARR> zdp = ranked_zeta<mint, LIM>(dp); FOR(i, N - d) { // zs[1<<i:2<<i], zdp[0:1<<i] vc<ARR> znewdp(1 << i); FOR(k, 1 << i) { FOR(p, i + 1) FOR(q, i - p + 1) { znewdp[k][p + q] += zdp[k][p] * zs[i][k][q]; } } auto x = ranked_mobius<mint, LIM>(znewdp); copy(all(x), newdp.begin() + (1 << i)); } swap(dp, newdp); } return dp; } // sum_i f_i s^i, s^i is subset-convolution template <typename mint, int LIM> vc<mint> sps_composition_poly(vc<mint> f, vc<mint> s) { const int N = topbit(len(s)); assert(len(s) == (1 << N)); if (f.empty()) return vc<mint>(1 << N, mint(0)); // convert to egf problem int D = min<int>(len(f) - 1, N); vc<mint> g(D + 1); mint c = s[0]; s[0] = 0; // (x+c)^i vc<mint> pow(D + 1); pow[0] = 1; FOR(i, len(f)) { FOR(j, D + 1) g[j] += f[i] * pow[j]; FOR_R(j, D + 1) pow[j] = pow[j] * c + (j == 0 ? mint(0) : pow[j - 1]); } // to egf mint factorial = 1; FOR(j, D + 1) g[j] *= factorial, factorial *= mint(j + 1); return sps_composition_egf<mint, LIM>(g, s); } #line 3 "setfunc/sps_log.hpp" // exp の逆手順で計算する template <typename mint, int LIM> vc<mint> sps_log(vc<mint>& dp) { const int N = topbit(len(dp)); assert(len(dp) == (1 << N) && dp[0] == mint(1)); vc<mint> s(1 << N); FOR_R(i, N) { vc<mint> a = {dp.begin() + (1 << i), dp.begin() + (2 << i)}; vc<mint> b = {dp.begin(), dp.begin() + (1 << i)}; auto RA = ranked_zeta<mint, LIM>(a); auto RB = ranked_zeta<mint, LIM>(b); FOR(s, 1 << i) { auto &f = RA[s], &g = RB[s]; // assert(g[0] == mint(1)); FOR(d, i + 1) { FOR(i, d) f[d] -= f[i] * g[d - i]; } } a = ranked_mobius<mint, LIM>(RA); copy(all(a), s.begin() + (1 << i)); } return s; }