This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub maspypy/library
#include "graph/tutte_polynomial.hpp"
#include "ds/unionfind/unionfind.hpp" #include "setfunc/sps_exp.hpp" #include "setfunc/sps_composition.hpp" template <typename mint, int NMAX> mint Tutte_Polynomial_Eval_connected(Graph<int, 0> G, mint X, mint Y) { int N = G.N; X -= 1, Y -= 1; /* V の連結成分分解を考えると, 各部分集合に S に対して, S を span する A の選び方に対する y^{cycle} の sum を F[S] として c[n]F^n/n!, c[n] = X^{n-k(E)} として EGF composition でできる. F[S] の計算 1 点ずつ足していく 集合に辺の個数に応じた重みをつけて exp 重み C(N,1) + C(N,2)Y + C(N,3)YY+... */ vv(mint, bin, N + 1, N + 1); bin[0][0] = 1; FOR(i, N) FOR(j, i + 1) bin[i + 1][j] += bin[i][j], bin[i + 1][j + 1] += bin[i][j]; vc<mint> wt(N + 1); FOR(n, 1, N + 1) { mint pow = 1; FOR(m, 1, n + 1) { wt[n] += bin[n][m] * pow, pow *= Y; } } vc<mint> F(1 << N); FOR(v, N) { u32 nbd = 0; for (auto& e: G[v]) nbd |= 1 << e.to; vc<mint> f(1 << v); FOR(s, 1 << v) { f[s] = F[s] * wt[popcnt(s & nbd)]; } f = sps_exp<mint, NMAX>(f); FOR(s, 1 << v) { F[s | 1 << v] = f[s]; } } if (X == mint(0)) { return F.back(); } // X で割れないときはこうすれば動く. 何もかもが環で動作する. // vc<mint> c(N + 1); // mint pow = 1; // FOR(n, 1, N + 1) { c[n] = pow, pow *= X; } // F = sps_composition_egf<mint, NMAX>(c, F); // return F.back(); FOR(s, 1 << N) F[s] *= X; F = sps_exp<mint, NMAX>(F); return F.back() * X.inverse(); } // QOJ 45 template <typename mint, int NMAX> mint Tutte_Polynomial_Eval(Graph<int, 0> G, mint X, mint Y) { int N = G.N; UnionFind uf(N); for (auto& e: G.edges) uf.merge(e.frm, e.to); vvc<int> vs(N); FOR(v, N) vs[uf[v]].eb(v); mint ANS = 1; for (auto& V: vs) { if (V.empty()) continue; Graph<int, 0> H = G.rearrange(V); ANS *= Tutte_Polynomial_Eval_connected<mint, NMAX>(H, X, Y); } return ANS; }
#line 2 "ds/unionfind/unionfind.hpp" struct UnionFind { int n, n_comp; vc<int> dat; // par or (-size) UnionFind(int n = 0) { build(n); } void build(int m) { n = m, n_comp = m; dat.assign(n, -1); } void reset() { build(n); } int operator[](int x) { while (dat[x] >= 0) { int pp = dat[dat[x]]; if (pp < 0) { return dat[x]; } x = dat[x] = pp; } return x; } ll size(int x) { x = (*this)[x]; return -dat[x]; } bool merge(int x, int y) { x = (*this)[x], y = (*this)[y]; if (x == y) return false; if (-dat[x] < -dat[y]) swap(x, y); dat[x] += dat[y], dat[y] = x, n_comp--; return true; } vc<int> get_all() { vc<int> A(n); FOR(i, n) A[i] = (*this)[i]; return A; } }; #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 "setfunc/sps_exp.hpp" // sum_i f_i/i! s^i, s^i is subset-convolution template <typename mint, int LIM> vc<mint> sps_exp(vc<mint>& s) { const int N = topbit(len(s)); assert(len(s) == (1 << N) && s[0] == mint(0)); vc<mint> dp(1 << N); dp[0] = mint(1); FOR(i, N) { vc<mint> a = {s.begin() + (1 << i), s.begin() + (2 << i)}; vc<mint> b = {dp.begin(), dp.begin() + (1 << i)}; a = subset_convolution<mint, LIM>(a, b); copy(all(a), dp.begin() + (1 << i)); } return dp; } #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 4 "graph/tutte_polynomial.hpp" template <typename mint, int NMAX> mint Tutte_Polynomial_Eval_connected(Graph<int, 0> G, mint X, mint Y) { int N = G.N; X -= 1, Y -= 1; /* V の連結成分分解を考えると, 各部分集合に S に対して, S を span する A の選び方に対する y^{cycle} の sum を F[S] として c[n]F^n/n!, c[n] = X^{n-k(E)} として EGF composition でできる. F[S] の計算 1 点ずつ足していく 集合に辺の個数に応じた重みをつけて exp 重み C(N,1) + C(N,2)Y + C(N,3)YY+... */ vv(mint, bin, N + 1, N + 1); bin[0][0] = 1; FOR(i, N) FOR(j, i + 1) bin[i + 1][j] += bin[i][j], bin[i + 1][j + 1] += bin[i][j]; vc<mint> wt(N + 1); FOR(n, 1, N + 1) { mint pow = 1; FOR(m, 1, n + 1) { wt[n] += bin[n][m] * pow, pow *= Y; } } vc<mint> F(1 << N); FOR(v, N) { u32 nbd = 0; for (auto& e: G[v]) nbd |= 1 << e.to; vc<mint> f(1 << v); FOR(s, 1 << v) { f[s] = F[s] * wt[popcnt(s & nbd)]; } f = sps_exp<mint, NMAX>(f); FOR(s, 1 << v) { F[s | 1 << v] = f[s]; } } if (X == mint(0)) { return F.back(); } // X で割れないときはこうすれば動く. 何もかもが環で動作する. // vc<mint> c(N + 1); // mint pow = 1; // FOR(n, 1, N + 1) { c[n] = pow, pow *= X; } // F = sps_composition_egf<mint, NMAX>(c, F); // return F.back(); FOR(s, 1 << N) F[s] *= X; F = sps_exp<mint, NMAX>(F); return F.back() * X.inverse(); } // QOJ 45 template <typename mint, int NMAX> mint Tutte_Polynomial_Eval(Graph<int, 0> G, mint X, mint Y) { int N = G.N; UnionFind uf(N); for (auto& e: G.edges) uf.merge(e.frm, e.to); vvc<int> vs(N); FOR(v, N) vs[uf[v]].eb(v); mint ANS = 1; for (auto& V: vs) { if (V.empty()) continue; Graph<int, 0> H = G.rearrange(V); ANS *= Tutte_Polynomial_Eval_connected<mint, NMAX>(H, X, Y); } return ANS; }