Polynomial Composite Set Power Series
(set/composite-set-power-series.hpp)

• View this file on GitHub
• Last update: 2025-10-29 02:30:28+09:00
• Include: #include "set/composite-set-power-series.hpp"

多項式 $\sum_{i=0}^{m-1}c_ix^i$，$2^n$ の集合冪級数 $a$ について $\sum_{i=0}^{m-1}c_ia^i$ を $O(n^22^m)$ 時間で計算する．

EGF との合成

$[x^0]a=0$ のとき EGF $f(x)=\sum_{k\geq 0}f_k\frac{x^k}{k!}$ との合成 $f(a)=\sum_{k\geq 0}f_k\frac{a^k}{k!}$ を考える．

• [https://codeforces.com/blog/entry/92183]

$R[x_1,\dots,x_n]/(x_1^2,\dots,x_n^2)$ の元として考えると $[x_n^1]f(a)=[x_n^0]f’(a)[x_n^1]a$ が成り立つ．

従って $[x_n^0]f’(a)$ から $f(a)$ を得られる．

再帰的に計算すると長さ $2^i$ の subset convolution を $n-i$ 回行うことになり，計算量は $O(n^22^n)$ と評価できる．

多項式との合成

$[x^0]a=0$ のとき EGF との合成に帰着される．

$[x^0]a=c\neq 0$ のとき $f(x+c)$ と $a-c$ を合成すればよい． $f(x+c)$ は $n$ 次まで求めれば十分なので $O(nm)$ 時間で計算できる．

別の手法

• [https://atcoder.jp/contests/xmascon20/tasks/xmascon20_h]
• [https://hos-lyric.hatenablog.com/entry/2021/01/14/201231]

ゼータ変換した先の多項式を FPS と合成することでも求められる．

log や pow は定数倍が改善する．

Depends on

• Subset Convolution (set/subset-convolution.hpp)

Verified with

• verify/set/LC_polynomial_composite_set_power_series.test.cpp
• verify/set/UNIT_composite_set_power_series.test.cpp

Code

#pragma once
#include "set/subset-convolution.hpp"

namespace SetPowerSeries {
template <class mint, int sz = 21>
vector<mint> composite_egf(vector<mint> f, vector<mint> a) {
  static SubsetConvolution<mint, sz> sc;
  assert(a[0] == 0);
  if (f.empty()) return vector<mint>(a.size());
  int l = __builtin_ctz(a.size());
  f.resize(l + 1);
  vector<vector<mint>> g(l + 1);
  for (int i = 0; i <= l; i++) g[i] = vector<mint>(1 << (l - i), 0);
  for (int i = 0; i <= l; i++) g[i][0] = f[i];
  for (int k = 0; k < l; k++) {
    auto p = sc.lift(vector<mint>(a.begin() + (1 << k), a.begin() + (2 << k)));
    sc.ranked_zeta(p);
    for (int i = 0; i < l - k; i++) {
      auto q = sc.lift(vector<mint>(g[i + 1].begin(), g[i + 1].begin() + (1 << k)));
      sc.ranked_zeta(q);
      sc.ranked_mul(q, p);
      sc.ranked_mobius(q);
      auto h = sc.unlift(q);
      copy(h.begin(), h.end(), g[i].begin() + (1 << k));
    }
  }
  return g[0];
}
template <class mint, int sz = 21>
vector<mint> composite_polynomial(vector<mint> p, vector<mint> a) {
  if (p.empty()) return vector<mint>(a.size());
  int l = __builtin_ctz(a.size());
  if (a[0] != 0) {
    mint c = a[0];
    a[0] = 0;
    vector<mint> p1(l + 1, 0), binom(l + 1, 0);
    binom[0] = 1;
    for (int i = 0; i < p.size(); i++) {
      mint r = i <= l ? 1 : c.pow(i - l);
      for (int j = min(i, l); j >= 0; j--, r *= c) p1[j] += p[i] * binom[j] * r;
      for (int j = l; j > 0; j--) binom[j] += binom[j - 1];
    }
    swap(p, p1);
  }
  mint r = 1;
  for (int i = 1; i <= l; i++) p[i] *= (r *= i);
  return composite_egf<mint, sz>(p, a);
}

// log(a), [x^0]a=1
// require inverse of 1,...,sz
template <class mint, int sz = 21>
vector<mint> log(vector<mint> a) {
  static SubsetConvolution<mint, sz> sc;
  assert(a[0] == 1);
  int l = __builtin_ctz(a.size());
  if (l == 0) return {0};
  vector<mint> inv(l + 1, 1);
  rep(i, 1, l + 1) inv[i] = mint(i).inv();
  auto p = sc.lift(a);
  sc.ranked_zeta(p);
  for (int k = 0; k < p.size(); k++) {
    auto q = p[k];
    p[k][0] = 0;
    for (int i = 1; i <= l; i++) {
      mint v = i * q[i];
      for (int j = 1; j < i; j++) v -= j * p[k][j] * q[i - j];
      p[k][i] = v * inv[i];
    }
  }
  sc.ranked_mobius(p);
  return sc.unlift(p);
}
// log(a), [x^0]a=1
// not require inverse of 1,...,sz
template <class mint, int sz = 21>
vector<mint> log_arbitrary(vector<mint> a) {
  assert(a[0] == 1);
  int l = __builtin_ctz(a.size());
  if (l == 0) return {0};
  a[0] = 0;
  vector<mint> f(l + 1, 0);
  f[1] = 1;
  for (int i = 2; i <= l; i++) f[i] = f[i - 1] * (1 - i);
  return composite_egf<mint, sz>(f, a);
}
// a^m, [x^0]a=1
// require inverse of 1,...,sz
template <class mint, int sz = 21>
vector<mint> pow(vector<mint> a, mint m) {
  static SubsetConvolution<mint, sz> sc;
  assert(a[0] == 1);
  int l = __builtin_ctz(a.size());
  if (l == 0) return {1};
  vector<mint> inv(l + 1, 1);
  rep(i, 1, l + 1) inv[i] = mint(i).inv();
  auto p = sc.lift(a);
  sc.ranked_zeta(p);
  for (int k = 0; k < p.size(); k++) {
    auto q = p[k];
    p[k][0] = 1;
    for (int i = 1; i <= l; i++) {
      mint v = 0;
      for (int j = 1; j < i; j++) v += (m * j - (i - j)) * p[k][i - j] * q[j];
      v *= inv[i];
      v += m * p[k][0] * q[i];
      p[k][i] = v;
    }
  }
  sc.ranked_mobius(p);
  return sc.unlift(p);
}
// a^m, [x^0]a=1
// not require inverse of 1,...,sz
template <class mint, int sz = 21>
vector<mint> pow_arbitrary(vector<mint> a, mint m) {
  assert(a[0] == 1);
  int l = __builtin_ctz(a.size());
  if (l == 0) return {1};
  a[0] = 0;
  vector<mint> f(l + 1, 0);
  f[0] = 1;
  for (int i = 1; i <= l; i++) {
    f[i] = f[i - 1] * m;
    m -= 1;
  }
  return composite_egf<mint, sz>(f, a);
}
};  // namespace SetPowerSeries
/**
 * @brief Polynomial Composite Set Power Series
 * @docs docs/set/composite-set-power-series.md
 */

#line 2 "set/subset-convolution.hpp"

template <class mint, int n_>
struct SubsetConvolution {
  static constexpr int n = n_;
  using poly = array<mint, n_ + 1>;
  vector<int> pc;
  SubsetConvolution() {
    pc.assign(1 << n, 0);
    for (int i = 1; i < pc.size(); i++) pc[i] = pc[i >> 1] + (i & 1);
  }
  void poly_add(poly& p, const poly& q, int d) {
    for (int i = 0; i < d; i++) p[i] += q[i];
  }
  void poly_sub(poly& p, const poly& q, int d) {
    for (int i = d; i <= n; i++) p[i] -= q[i];
  }
  void poly_mul(poly& p, const poly& q) {
    poly r{};
    for (int i = 0; i <= n; i++)
      for (int j = 0; j <= n - i; j++)
        r[i + j] += p[i] * q[j];
    swap(p, r);
  }
  vector<poly> lift(const vector<mint>& a) {
    int n = a.size();
    assert(n == (n & -n));
    vector<poly> b(n);
    for (int i = 0; i < n; i++) {
      b[i].fill(0);
      b[i][pc[i]] = a[i];
    }
    return b;
  }
  vector<mint> unlift(const vector<poly>& b) {
    int n = b.size();
    assert(n == (n & -n));
    vector<mint> a(n);
    for (int i = 0; i < n; i++) a[i] = b[i][pc[i]];
    return a;
  }
  void ranked_zeta(vector<poly>& a) {
    int n = a.size();
    for (int i = 1; i < n; i <<= 1)
      for (int j = 0; j < n; j += i * 2)
        for (int k = 0; k < i; k++)
          poly_add(a[i + j + k], a[j + k], pc[i + j + k]);
  }
  void ranked_mobius(vector<poly>& a) {
    int n = a.size();
    for (int i = 1; i < n; i <<= 1)
      for (int j = 0; j < n; j += i * 2)
        for (int k = 0; k < i; k++)
          poly_sub(a[i + j + k], a[j + k], pc[i + j + k]);
  }
  void ranked_mul(vector<poly>& a, const vector<poly>& b) {
    for (int i = 0; i < a.size(); i++) poly_mul(a[i], b[i]);
  }
  vector<mint> multiply(const vector<mint>& a, const vector<mint>& b) {
    auto p = lift(a);
    auto q = lift(b);
    ranked_zeta(p);
    ranked_zeta(q);
    ranked_mul(p, q);
    ranked_mobius(p);
    return unlift(p);
  }
};

/**
 * @brief Subset Convolution
 * @docs docs/set/subset-convolution.md
 */
#line 3 "set/composite-set-power-series.hpp"

namespace SetPowerSeries {
template <class mint, int sz = 21>
vector<mint> composite_egf(vector<mint> f, vector<mint> a) {
  static SubsetConvolution<mint, sz> sc;
  assert(a[0] == 0);
  if (f.empty()) return vector<mint>(a.size());
  int l = __builtin_ctz(a.size());
  f.resize(l + 1);
  vector<vector<mint>> g(l + 1);
  for (int i = 0; i <= l; i++) g[i] = vector<mint>(1 << (l - i), 0);
  for (int i = 0; i <= l; i++) g[i][0] = f[i];
  for (int k = 0; k < l; k++) {
    auto p = sc.lift(vector<mint>(a.begin() + (1 << k), a.begin() + (2 << k)));
    sc.ranked_zeta(p);
    for (int i = 0; i < l - k; i++) {
      auto q = sc.lift(vector<mint>(g[i + 1].begin(), g[i + 1].begin() + (1 << k)));
      sc.ranked_zeta(q);
      sc.ranked_mul(q, p);
      sc.ranked_mobius(q);
      auto h = sc.unlift(q);
      copy(h.begin(), h.end(), g[i].begin() + (1 << k));
    }
  }
  return g[0];
}
template <class mint, int sz = 21>
vector<mint> composite_polynomial(vector<mint> p, vector<mint> a) {
  if (p.empty()) return vector<mint>(a.size());
  int l = __builtin_ctz(a.size());
  if (a[0] != 0) {
    mint c = a[0];
    a[0] = 0;
    vector<mint> p1(l + 1, 0), binom(l + 1, 0);
    binom[0] = 1;
    for (int i = 0; i < p.size(); i++) {
      mint r = i <= l ? 1 : c.pow(i - l);
      for (int j = min(i, l); j >= 0; j--, r *= c) p1[j] += p[i] * binom[j] * r;
      for (int j = l; j > 0; j--) binom[j] += binom[j - 1];
    }
    swap(p, p1);
  }
  mint r = 1;
  for (int i = 1; i <= l; i++) p[i] *= (r *= i);
  return composite_egf<mint, sz>(p, a);
}

// log(a), [x^0]a=1
// require inverse of 1,...,sz
template <class mint, int sz = 21>
vector<mint> log(vector<mint> a) {
  static SubsetConvolution<mint, sz> sc;
  assert(a[0] == 1);
  int l = __builtin_ctz(a.size());
  if (l == 0) return {0};
  vector<mint> inv(l + 1, 1);
  rep(i, 1, l + 1) inv[i] = mint(i).inv();
  auto p = sc.lift(a);
  sc.ranked_zeta(p);
  for (int k = 0; k < p.size(); k++) {
    auto q = p[k];
    p[k][0] = 0;
    for (int i = 1; i <= l; i++) {
      mint v = i * q[i];
      for (int j = 1; j < i; j++) v -= j * p[k][j] * q[i - j];
      p[k][i] = v * inv[i];
    }
  }
  sc.ranked_mobius(p);
  return sc.unlift(p);
}
// log(a), [x^0]a=1
// not require inverse of 1,...,sz
template <class mint, int sz = 21>
vector<mint> log_arbitrary(vector<mint> a) {
  assert(a[0] == 1);
  int l = __builtin_ctz(a.size());
  if (l == 0) return {0};
  a[0] = 0;
  vector<mint> f(l + 1, 0);
  f[1] = 1;
  for (int i = 2; i <= l; i++) f[i] = f[i - 1] * (1 - i);
  return composite_egf<mint, sz>(f, a);
}
// a^m, [x^0]a=1
// require inverse of 1,...,sz
template <class mint, int sz = 21>
vector<mint> pow(vector<mint> a, mint m) {
  static SubsetConvolution<mint, sz> sc;
  assert(a[0] == 1);
  int l = __builtin_ctz(a.size());
  if (l == 0) return {1};
  vector<mint> inv(l + 1, 1);
  rep(i, 1, l + 1) inv[i] = mint(i).inv();
  auto p = sc.lift(a);
  sc.ranked_zeta(p);
  for (int k = 0; k < p.size(); k++) {
    auto q = p[k];
    p[k][0] = 1;
    for (int i = 1; i <= l; i++) {
      mint v = 0;
      for (int j = 1; j < i; j++) v += (m * j - (i - j)) * p[k][i - j] * q[j];
      v *= inv[i];
      v += m * p[k][0] * q[i];
      p[k][i] = v;
    }
  }
  sc.ranked_mobius(p);
  return sc.unlift(p);
}
// a^m, [x^0]a=1
// not require inverse of 1,...,sz
template <class mint, int sz = 21>
vector<mint> pow_arbitrary(vector<mint> a, mint m) {
  assert(a[0] == 1);
  int l = __builtin_ctz(a.size());
  if (l == 0) return {1};
  a[0] = 0;
  vector<mint> f(l + 1, 0);
  f[0] = 1;
  for (int i = 1; i <= l; i++) {
    f[i] = f[i - 1] * m;
    m -= 1;
  }
  return composite_egf<mint, sz>(f, a);
}
};  // namespace SetPowerSeries
/**
 * @brief Polynomial Composite Set Power Series
 * @docs docs/set/composite-set-power-series.md
 */

Back to top page