graph/tutte_polynomial.hpp

View this file on GitHub
Last update: 2026-08-29 09:24:19+09:00
Include: #include "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 1 "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 1 "setfunc/sps_exp.hpp"

#line 1 "other/bit.hpp"

int popcnt(int x) { return __builtin_popcount(x); }
int popcnt(u32 x) { return __builtin_popcount(x); }
int popcnt(ll x) { return __builtin_popcountll(x); }
int popcnt(u64 x) { return __builtin_popcountll(x); }
int popcnt_sgn(int x) { return (__builtin_parity(unsigned(x)) & 1 ? -1 : 1); }
int popcnt_sgn(u32 x) { return (__builtin_parity(x) & 1 ? -1 : 1); }
int popcnt_sgn(ll x) { return (__builtin_parityll(x) & 1 ? -1 : 1); }
int popcnt_sgn(u64 x) { return (__builtin_parityll(x) & 1 ? -1 : 1); }
// (0, 1, 2, 3, 4) -> (-1, 0, 1, 1, 2)
int topbit(int x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
int topbit(u32 x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
int topbit(ll x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
int topbit(u64 x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
// (0, 1, 2, 3, 4) -> (-1, 0, 1, 0, 2)
int lowbit(int x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
int lowbit(u32 x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
int lowbit(ll x) { return (x == 0 ? -1 : __builtin_ctzll(x)); }
int lowbit(u64 x) { return (x == 0 ? -1 : __builtin_ctzll(x)); }

template <typename T>
T kth_bit(int k) {
  assert(0 <= k && k < int(8 * sizeof(T)));
  return T(1) << k;
}
template <typename T>
bool has_kth_bit(T x, int k) {
  assert(0 <= k && k < int(8 * sizeof(T)));
  return x >> k & 1;
}

template <typename UINT>
struct all_bit {
  static_assert(is_unsigned<UINT>::value);
  UINT s;
  all_bit(UINT s) : s(s) {}
  struct iter {
    UINT s;
    int operator*() const { return lowbit(s); }
    void operator++() { s &= s - 1; }
    bool operator!=(nullptr_t) const { return s; }
  };
  iter begin() const { return {s}; }
  nullptr_t end() const { return nullptr; }
};

template <typename UINT>
struct all_subset {
  static_assert(is_unsigned<UINT>::value);
  UINT s;
  all_subset(UINT s) : s(s) {}
  struct iter {
    UINT s, t;
    bool done = false;
    UINT operator*() const { return t; }
    void operator++() {
      done = (t == 0);
      t = (t - 1) & s;
    }
    bool operator!=(nullptr_t) const { return !done; }
  };
  iter begin() const { return {s, s}; }
  nullptr_t end() const { return nullptr; }
};

constexpr u64 full_mask(int n) {
  assert(0 <= n && n <= 64);
  return n == 64 ? -1ULL : (1ULL << n) - 1;
}

u64 bit_reverse(u64 x) {
  x = ((x & 0x5555555555555555ULL) << 1) | ((x >> 1) & 0x5555555555555555ULL);
  x = ((x & 0x3333333333333333ULL) << 2) | ((x >> 2) & 0x3333333333333333ULL);
  x = ((x & 0x0f0f0f0f0f0f0f0fULL) << 4) | ((x >> 4) & 0x0f0f0f0f0f0f0f0fULL);
  x = ((x & 0x00ff00ff00ff00ffULL) << 8) | ((x >> 8) & 0x00ff00ff00ff00ffULL);
  x = ((x & 0x0000ffff0000ffffULL) << 16) | ((x >> 16) & 0x0000ffff0000ffffULL);
  x = (x << 32) | (x >> 32);
  return x;
}
#line 1 "setfunc/ranked_zeta.hpp"

#line 1 "setfunc/bitwise_transform.hpp"

namespace bitwise {

enum class trans_type {
  hadamard,
  superset_zeta,
  superset_mobius,
  subset_zeta,
  subset_mobius,
  ranked_zeta,
  ranked_mobius,
  superset_zeta_or
};

template <typename ARR>
inline void ranked_add(ARR& a, const ARR& b) {
  for (int d = 0; d < int(a.size()); ++d) a[d] += b[d];
}

template <typename ARR>
inline void ranked_sub(ARR& a, const ARR& b) {
  for (int d = 0; d < int(a.size()); ++d) a[d] -= b[d];
}

template <trans_type type, int N, typename T>
inline void bitwise_transform_fixed(T* a) {
  static_assert(N >= 1 && (N & (N - 1)) == 0);
  if constexpr (N == 1) {
    return;
  } else {
    constexpr int H = N / 2;
    bitwise_transform_fixed<type, H>(a);
    bitwise_transform_fixed<type, H>(a + H);
    if constexpr (type == trans_type::hadamard) {
      for (int i = 0; i < H; ++i) {
        auto x = a[i], y = a[H + i];
        a[i] = x + y, a[H + i] = x - y;
      }
    }
    if constexpr (type == trans_type::superset_zeta) {
      for (int i = 0; i < H; ++i) a[i] += a[H + i];
    }
    if constexpr (type == trans_type::superset_mobius) {
      for (int i = 0; i < H; ++i) a[i] -= a[H + i];
    }
    if constexpr (type == trans_type::subset_zeta) {
      for (int i = 0; i < H; ++i) a[H + i] += a[i];
    }
    if constexpr (type == trans_type::subset_mobius) {
      for (int i = 0; i < H; ++i) a[H + i] -= a[i];
    }
    if constexpr (type == trans_type::ranked_zeta) {
      for (int i = 0; i < H; ++i) ranked_add(a[H + i], a[i]);
    }
    if constexpr (type == trans_type::ranked_mobius) {
      for (int i = 0; i < H; ++i) ranked_sub(a[H + i], a[i]);
    }
    if constexpr (type == trans_type::superset_zeta_or) {
      for (int i = 0; i < H; ++i) a[i] |= a[H + i];
    }
  }
}

template <trans_type type, int N, typename T>
inline void bitwise_transform_dispatch(vc<T>& a) {
  if (len(a) == N) {
    return bitwise_transform_fixed<type, N>(a.data());
  }
  if constexpr (N > 1) {
    return bitwise_transform_dispatch<type, N / 2>(a);
  }
}

template <trans_type type, typename T>
inline void bitwise_transform(vc<T>& a) {
  int n = len(a);
  assert(n >= 1);
  assert((n & (n - 1)) == 0);
  assert(n <= (1 << 25));
  bitwise_transform_dispatch<type, 1 << 25>(a);
}
}  // namespace bitwise
#line 3 "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];
  bitwise::bitwise_transform<bitwise::trans_type::ranked_zeta>(Rf);
  return Rf;
}

template <typename T, int LIM>
vc<T> ranked_mobius(vc<array<T, LIM + 1>>& Rf) {
  bitwise::bitwise_transform<bitwise::trans_type::ranked_mobius>(Rf);
  vc<T> f(len(Rf));
  for (int s = 0; s < len(f); ++s) f[s] = Rf[s][popcnt(s)];
  return f;
}
#line 3 "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 3 "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 1 "setfunc/ranked_zeta.hpp"

#line 1 "setfunc/bitwise_transform.hpp"

namespace bitwise {

enum class trans_type {
  hadamard,
  superset_zeta,
  superset_mobius,
  subset_zeta,
  subset_mobius,
  ranked_zeta,
  ranked_mobius,
  superset_zeta_or
};

template <typename ARR>
inline void ranked_add(ARR& a, const ARR& b) {
  for (int d = 0; d < int(a.size()); ++d) a[d] += b[d];
}

template <typename ARR>
inline void ranked_sub(ARR& a, const ARR& b) {
  for (int d = 0; d < int(a.size()); ++d) a[d] -= b[d];
}

template <trans_type type, int N, typename T>
inline void bitwise_transform_fixed(T* a) {
  static_assert(N >= 1 && (N & (N - 1)) == 0);
  if constexpr (N == 1) {
    return;
  } else {
    constexpr int H = N / 2;
    bitwise_transform_fixed<type, H>(a);
    bitwise_transform_fixed<type, H>(a + H);
    if constexpr (type == trans_type::hadamard) {
      for (int i = 0; i < H; ++i) {
        auto x = a[i], y = a[H + i];
        a[i] = x + y, a[H + i] = x - y;
      }
    }
    if constexpr (type == trans_type::superset_zeta) {
      for (int i = 0; i < H; ++i) a[i] += a[H + i];
    }
    if constexpr (type == trans_type::superset_mobius) {
      for (int i = 0; i < H; ++i) a[i] -= a[H + i];
    }
    if constexpr (type == trans_type::subset_zeta) {
      for (int i = 0; i < H; ++i) a[H + i] += a[i];
    }
    if constexpr (type == trans_type::subset_mobius) {
      for (int i = 0; i < H; ++i) a[H + i] -= a[i];
    }
    if constexpr (type == trans_type::ranked_zeta) {
      for (int i = 0; i < H; ++i) ranked_add(a[H + i], a[i]);
    }
    if constexpr (type == trans_type::ranked_mobius) {
      for (int i = 0; i < H; ++i) ranked_sub(a[H + i], a[i]);
    }
    if constexpr (type == trans_type::superset_zeta_or) {
      for (int i = 0; i < H; ++i) a[i] |= a[H + i];
    }
  }
}

template <trans_type type, int N, typename T>
inline void bitwise_transform_dispatch(vc<T>& a) {
  if (len(a) == N) {
    return bitwise_transform_fixed<type, N>(a.data());
  }
  if constexpr (N > 1) {
    return bitwise_transform_dispatch<type, N / 2>(a);
  }
}

template <trans_type type, typename T>
inline void bitwise_transform(vc<T>& a) {
  int n = len(a);
  assert(n >= 1);
  assert((n & (n - 1)) == 0);
  assert(n <= (1 << 25));
  bitwise_transform_dispatch<type, 1 << 25>(a);
}
}  // namespace bitwise
#line 3 "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];
  bitwise::bitwise_transform<bitwise::trans_type::ranked_zeta>(Rf);
  return Rf;
}

template <typename T, int LIM>
vc<T> ranked_mobius(vc<array<T, LIM + 1>>& Rf) {
  bitwise::bitwise_transform<bitwise::trans_type::ranked_mobius>(Rf);
  vc<T> f(len(Rf));
  for (int s = 0; s < len(f); ++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 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;
}