setfunc/sps_composition.hpp

• View this file on GitHub
• Last update: 2026-09-13 17:20:33+09:00
• Include: #include "setfunc/sps_composition.hpp"

Depends on

• setfunc/bitwise_transform.hpp
• setfunc/ranked_zeta.hpp

Required by

• graph/count/count_biconnected_subgraph.hpp
• graph/count/count_bridgeless_subgraph.hpp
• graph/count/count_connected_subgraph.hpp
• graph/tutte_polynomial.hpp
• setfunc/sps_log.hpp

Verified with

• test/2_library_checker/set_power_series/composition_sps.test.cpp

Code

#include "setfunc/ranked_zeta.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 (f.empty()) return vc<mint>(1 << N, 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 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 (f.empty()) return vc<mint>(1 << N, 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);
}

Back to top page