library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub maspypy/library

:warning: graph/tutte_polynomial.hpp

Depends on

Code

#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/sps_exp.hpp"

#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 4 "setfunc/sps_exp.hpp"

// sum_i 1/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 3 "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;
}
Back to top page