Skip to the content.

:heavy_check_mark: Power Projection Of Set Power Series
(set/power-projection-of-set-power-series.hpp)

$a$ を集合冪級数とする. $\sum_{s}w_s[x^s]a^i$ を $i=0,\dots,M-1$ について列挙する.

$O(N^22^N+NM)$ 時間.

導入

$a$ を fix したとき,合成 $f(a)$ の計算は $\sum_{i}f_i[x^s]a^i$ の $s$ についての列挙.

この転置問題は $\sum_{s}w_s[x^s]a^i$ の $i$ についての列挙.

多項式 $f$ との合成は以下の要領で計算していた.

$g(x)=f(x+c)$ のとき $g_i=\sum_{j\geq i}\binom{j}{i}c^{j-i}f_j$ である.この転置は難しくない.

また $g$ を EGF にするには $x^k$ の係数を $k!$ 倍すればよいので,特に転置は同じものである.

$[x^0]a=0$ のときの EGF との合成の転置を考える. subset convolution $a\star b$ は $b$ を固定すれば $a$ について線形写像. その転置アルゴリズムは $a\star b=\operatorname{rev}(\operatorname{rev}(a)\star b)$ である.

Depends on

Verified with

Code

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

namespace SetPowerSeries {
template <class mint, int sz = 21>
vector<mint> power_projection(const vector<mint>& a, const vector<mint>& w, int m) {
  static SubsetConvolution<mint, sz> sc;
  int l = __builtin_ctz(a.size());
  assert(a.size() == (1 << l));
  mint c = a[0];
  vector<vector<mint>> g(l + 1);
  g[0] = w;
  for (int i = 1; i <= l; i++) g[i] = vector<mint>(1 << (l - i), 0);
  for (int k = l - 1; k >= 0; 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 h = vector<mint>(g[i].begin() + (1 << k), g[i].begin() + (2 << k));
      reverse(h.begin(), h.end());
      auto q = sc.lift(h);
      sc.ranked_zeta(q);
      sc.ranked_mul(q, p);
      sc.ranked_mobius(q);
      h = sc.unlift(q);
      reverse(h.begin(), h.end());
      for (int j = 0; j < h.size(); j++) g[i + 1][j] += h[j];
    }
  }
  vector<mint> f1(l + 1, 1);
  for (int i = 1; i <= l; i++) f1[i] = f1[i - 1] * i;
  for (int i = 0; i <= l; i++) f1[i] *= g[i][0];
  vector<mint> f(m, 0), binom(l + 1, 0);
  binom[0] = 1;
  mint r0 = 1;
  for (int i = 0; i < m; i++) {
    if (i > l) r0 *= c;
    mint r = r0;
    for (int j = min(i, l); j >= 0; j--, r *= c) f[i] += f1[j] * binom[j] * r;
    for (int j = l; j > 0; j--) binom[j] += binom[j - 1];
  }
  return f;
}
};  // namespace SetPowerSeries
/**
 * @brief Power Projection Of Set Power Series
 * @docs docs/set/power-projection-of-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/power-projection-of-set-power-series.hpp"

namespace SetPowerSeries {
template <class mint, int sz = 21>
vector<mint> power_projection(const vector<mint>& a, const vector<mint>& w, int m) {
  static SubsetConvolution<mint, sz> sc;
  int l = __builtin_ctz(a.size());
  assert(a.size() == (1 << l));
  mint c = a[0];
  vector<vector<mint>> g(l + 1);
  g[0] = w;
  for (int i = 1; i <= l; i++) g[i] = vector<mint>(1 << (l - i), 0);
  for (int k = l - 1; k >= 0; 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 h = vector<mint>(g[i].begin() + (1 << k), g[i].begin() + (2 << k));
      reverse(h.begin(), h.end());
      auto q = sc.lift(h);
      sc.ranked_zeta(q);
      sc.ranked_mul(q, p);
      sc.ranked_mobius(q);
      h = sc.unlift(q);
      reverse(h.begin(), h.end());
      for (int j = 0; j < h.size(); j++) g[i + 1][j] += h[j];
    }
  }
  vector<mint> f1(l + 1, 1);
  for (int i = 1; i <= l; i++) f1[i] = f1[i - 1] * i;
  for (int i = 0; i <= l; i++) f1[i] *= g[i][0];
  vector<mint> f(m, 0), binom(l + 1, 0);
  binom[0] = 1;
  mint r0 = 1;
  for (int i = 0; i < m; i++) {
    if (i > l) r0 *= c;
    mint r = r0;
    for (int j = min(i, l); j >= 0; j--, r *= c) f[i] += f1[j] * binom[j] * r;
    for (int j = l; j > 0; j--) binom[j] += binom[j - 1];
  }
  return f;
}
};  // namespace SetPowerSeries
/**
 * @brief Power Projection Of Set Power Series
 * @docs docs/set/power-projection-of-set-power-series.md
 */
Back to top page