This documentation is automatically generated by online-judge-tools/verification-helper
#include "setfunc/sps_composition.hpp"
#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 (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 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 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);
}