set_power_series.hpp¶
Compose an exponential generating function with a set power series. A coefficient indexed by a mask represents the square-free monomial whose variables are that mask. Splitting masks by their largest variable and by popcount turns every multiplication into ranked subset convolution. The subset zeta transform performs all disjoint decompositions at once, while the rank selects exactly the square-free terms. egf[k] is the multiplier of series^k / k!, and the constant term of series must be zero.
Verified by exp_of_set_power_series, log_of_set_power_series, polynomial_composite_set_power_series, power_projection_of_set_power_series.
\[
\displaystyle F=\prod_i (1+x_i)
\]
Implementation¶
#ifndef NOYA_SET_POWER_SERIES_HPP
#define NOYA_SET_POWER_SERIES_HPP 1
/// @complexity Time: O(n^2 2^n) for EGF composition, exponential, and
/// logarithm; O(n^2 2^n + n m) for m-term power projection.
/// Space: O(n 2^n).
#include "noya/convolution.hpp"
#include <algorithm>
#include <cassert>
#include <vector>
namespace noya {
namespace set_power_series_internal {
template <class T, class Operation>
void subset_transform(std::vector<T> &values, Operation operation) {
int size = int(values.size());
for (int width = 1; width < size; width <<= 1) {
for (int block = 0; block < size; block += width << 1) {
for (int offset = 0; offset < width; offset++) {
operation(values[block + offset],
values[block + width + offset]);
}
}
}
}
template <class Mint>
std::vector<Mint>
transposed_subset_product(const std::vector<Mint> &series,
std::vector<Mint> weights) {
std::reverse(weights.begin(), weights.end());
weights = subset_convolution(weights, series);
std::reverse(weights.begin(), weights.end());
return weights;
}
template <class Mint>
std::vector<Mint>
transposed_egf_composition(const std::vector<Mint> &series,
const std::vector<Mint> &weights) {
int size = int(series.size());
assert(size > 0 && (size & (size - 1)) == 0);
assert(weights.size() == series.size());
assert(series[0] == Mint{});
int variables = size == 1 ? 0 : 31 - __builtin_clz(size);
std::vector<Mint> result(variables + 1);
result[0] = weights[0];
std::vector<Mint> state = weights;
for (int degree = 0; degree < variables; degree++) {
std::vector<Mint> next(1 << (variables - degree - 1));
for (int bit = 0; bit < variables - degree; bit++) {
int length = 1 << bit;
std::vector<Mint> block_series(series.begin() + length,
series.begin() + 2 * length);
std::vector<Mint> block_weights(state.begin() + length,
state.begin() + 2 * length);
block_weights =
transposed_subset_product(block_series, std::move(block_weights));
for (int mask = 0; mask < length; mask++) {
next[mask] += block_weights[mask];
}
}
state = std::move(next);
result[degree + 1] = state[0];
}
return result;
}
} // namespace set_power_series_internal
/// @brief Compose an exponential generating function with a set power series.
/// A coefficient indexed by a mask represents the square-free monomial whose
/// variables are that mask. Splitting masks by their largest variable and by
/// popcount turns every multiplication into ranked subset convolution. The
/// subset zeta transform performs all disjoint decompositions at once, while
/// the rank selects exactly the square-free terms. `egf[k]` is the multiplier
/// of `series^k / k!`, and the constant term of `series` must be zero.
template <class Mint>
std::vector<Mint>
set_power_series_egf_composition(const std::vector<Mint> &egf,
const std::vector<Mint> &series) {
assert(!series.empty() && (series.size() & (series.size() - 1)) == 0);
int variables = series.size() == 1
? 0
: 31 - __builtin_clz(unsigned(series.size()));
assert(int(egf.size()) == variables + 1);
assert(series[0] == Mint{});
std::vector<std::vector<std::vector<Mint>>> states(variables + 1);
for (int level = 0; level <= variables; level++) {
states[level].assign(level + 1,
std::vector<Mint>(std::size_t(1) << level));
states[level][0][0] = egf[variables - level];
}
for (int bit = 0; bit < variables; bit++) {
int half = 1 << bit;
std::vector<std::vector<Mint>> ranked_series(
bit + 1, std::vector<Mint>(half));
for (int mask = 0; mask < half; mask++) {
ranked_series[__builtin_popcount(unsigned(mask))][mask] =
series[half + mask];
}
for (auto &rank : ranked_series) {
set_power_series_internal::subset_transform(
rank, [](Mint &lower, Mint &upper) { upper += lower; });
}
for (int level = bit + 1; level <= variables; level++) {
for (int rank = 0; rank <= level; rank++) {
std::copy_n(states[level][rank].begin(), half,
states[level][rank].begin() + half);
}
}
for (int level = bit; level < variables; level++) {
for (int series_rank = 0; series_rank <= bit; series_rank++) {
for (int state_rank = 0;
series_rank + state_rank <= level; state_rank++) {
for (int mask = 0; mask < half; mask++) {
states[level + 1][series_rank + state_rank + 1][half + mask] +=
ranked_series[series_rank][mask] *
states[level][state_rank][mask];
}
}
}
}
}
for (auto &rank : states[variables]) {
set_power_series_internal::subset_transform(
rank, [](Mint &lower, Mint &upper) { upper -= lower; });
}
std::vector<Mint> result(series.size());
for (int mask = 0; mask < int(series.size()); mask++) {
result[mask] =
states[variables][__builtin_popcount(unsigned(mask))][mask];
}
return result;
}
/// @brief Substitute a set power series into an ordinary polynomial.
/// Taylor expansion at the constant coefficient `c=series[0]` gives the EGF
/// multipliers `f^(k)(c)`. Only the first n derivatives matter because every
/// square-free monomial has degree at most n; the ranked composition routine
/// then evaluates the nilpotent remainder `series-c`.
template <class Mint>
std::vector<Mint>
set_power_series_polynomial_composition(std::vector<Mint> polynomial,
const std::vector<Mint> &series) {
assert(!series.empty() && (series.size() & (series.size() - 1)) == 0);
int variables = series.size() == 1
? 0
: 31 - __builtin_clz(unsigned(series.size()));
std::vector<Mint> derivatives(variables + 1);
Mint constant = series[0];
for (int order = 0; order <= variables && !polynomial.empty(); order++) {
Mint value{};
for (int index = int(polynomial.size()) - 1; index >= 0; index--) {
value = value * constant + polynomial[index];
}
derivatives[order] = value;
for (int index = 1; index < int(polynomial.size()); index++) {
polynomial[index - 1] = polynomial[index] * Mint(index);
}
polynomial.pop_back();
}
std::vector<Mint> nilpotent = series;
nilpotent[0] = Mint{};
return set_power_series_egf_composition(derivatives, nilpotent);
}
/// @brief Return the exponential of a set power series with zero constant.
/// In the square-free quotient every term above degree n vanishes, so exp is
/// the finite EGF with every multiplier equal to one.
template <class Mint>
std::vector<Mint>
set_power_series_exponential(const std::vector<Mint> &series) {
assert(!series.empty() && (series.size() & (series.size() - 1)) == 0);
int variables = series.size() == 1
? 0
: 31 - __builtin_clz(unsigned(series.size()));
assert(series[0] == Mint{});
return set_power_series_egf_composition(
std::vector<Mint>(variables + 1, Mint(1)), series);
}
/// @brief Return the logarithm of a set power series with constant one.
/// Writing the input as `1+u`, the nilpotent series u satisfies
/// log(1+u)=sum (-1)^(k-1)u^k/k. Therefore its kth EGF multiplier is
/// `(-1)^(k-1)(k-1)!`, and ranked EGF composition evaluates all terms.
template <class Mint>
std::vector<Mint>
set_power_series_logarithm(const std::vector<Mint> &series) {
assert(!series.empty() && (series.size() & (series.size() - 1)) == 0);
int variables = series.size() == 1
? 0
: 31 - __builtin_clz(unsigned(series.size()));
assert(series[0] == Mint(1));
std::vector<Mint> multipliers(variables + 1);
Mint factorial = Mint(1);
for (int degree = 1; degree <= variables; degree++) {
multipliers[degree] = (degree & 1) ? factorial : -factorial;
factorial *= Mint(degree);
}
std::vector<Mint> nilpotent = series;
nilpotent[0] -= Mint(1);
return set_power_series_egf_composition(multipliers, nilpotent);
}
/// @brief Project consecutive powers of a set power series onto weights.
/// Reversing masks changes the transpose of subset convolution into an
/// ordinary subset convolution. Repeated transposed ranked composition gives
/// the projections of `u^k/k!` for `u=series-series[0]`; convolving those n+1
/// values with `c^j/j!` and multiplying by m! restores `(c+u)^m`.
template <class Mint>
std::vector<Mint>
set_power_series_power_projection(const std::vector<Mint> &series,
const std::vector<Mint> &weights,
int count) {
assert(count >= 0);
assert(!series.empty() && (series.size() & (series.size() - 1)) == 0);
assert(weights.size() == series.size());
if (count == 0) {
return {};
}
Mint constant = series[0];
std::vector<Mint> nilpotent = series;
nilpotent[0] = Mint{};
std::vector<Mint> projected =
set_power_series_internal::transposed_egf_composition(nilpotent,
weights);
std::vector<Mint> factorial(count, Mint(1));
std::vector<Mint> inverse_factorial(count, Mint(1));
for (int index = 1; index < count; index++) {
factorial[index] = factorial[index - 1] * Mint(index);
}
inverse_factorial.back() = Mint(1) / factorial.back();
for (int index = count - 1; index > 0; index--) {
inverse_factorial[index - 1] = inverse_factorial[index] * Mint(index);
}
std::vector<Mint> constant_egf(count);
Mint power = Mint(1);
for (int index = 0; index < count; index++) {
constant_egf[index] = power * inverse_factorial[index];
power *= constant;
}
std::vector<Mint> result(count);
for (int degree = 0; degree < int(projected.size()); degree++) {
for (int total = degree; total < count; total++) {
result[total] += projected[degree] * constant_egf[total - degree];
}
}
for (int index = 0; index < count; index++) {
result[index] *= factorial[index];
}
return result;
}
} // namespace noya
#endif // NOYA_SET_POWER_SERIES_HPP
#include <algorithm>
#include <array>
#include <cassert>
#include <numeric>
#include <type_traits>
#include <utility>
#include <vector>
/// @complexity Time: O(n^2 2^n) for EGF composition, exponential, and
/// logarithm; O(n^2 2^n + n m) for m-term power projection.
/// Space: O(n 2^n).
/// @complexity Time: O(n log n) for bitwise/divisor transforms; O(n log^2 n) for subset convolution.
/// Space: O(n log n) for subset layers, otherwise O(n).
#ifdef _MSC_VER
#include <intrin.h>
#endif
#if __cplusplus >= 202002L
#include <bit>
#endif
namespace atcoder {
namespace internal {
#if __cplusplus >= 202002L
using std::bit_ceil;
#else
// @return same with std::bit::bit_ceil
unsigned int bit_ceil(unsigned int n) {
unsigned int x = 1;
while (x < (unsigned int)(n)) x *= 2;
return x;
}
#endif
// @param n `1 <= n`
// @return same with std::bit::countr_zero
int countr_zero(unsigned int n) {
#ifdef _MSC_VER
unsigned long index;
_BitScanForward(&index, n);
return index;
#else
return __builtin_ctz(n);
#endif
}
// @param n `1 <= n`
// @return same with std::bit::countr_zero
constexpr int countr_zero_constexpr(unsigned int n) {
int x = 0;
while (!(n & (1 << x))) x++;
return x;
}
} // namespace internal
} // namespace atcoder
#ifdef _MSC_VER
#include <intrin.h>
#endif
#ifdef _MSC_VER
#include <intrin.h>
#endif
namespace atcoder {
namespace internal {
// @param m `1 <= m`
// @return x mod m
constexpr long long safe_mod(long long x, long long m) {
x %= m;
if (x < 0) x += m;
return x;
}
// Fast modular multiplication by barrett reduction
// Reference: https://en.wikipedia.org/wiki/Barrett_reduction
// NOTE: reconsider after Ice Lake
struct barrett {
unsigned int _m;
unsigned long long im;
// @param m `1 <= m`
explicit barrett(unsigned int m) : _m(m), im((unsigned long long)(-1) / m + 1) {}
// @return m
unsigned int umod() const { return _m; }
// @param a `0 <= a < m`
// @param b `0 <= b < m`
// @return `a * b % m`
unsigned int mul(unsigned int a, unsigned int b) const {
// [1] m = 1
// a = b = im = 0, so okay
// [2] m >= 2
// im = ceil(2^64 / m)
// -> im * m = 2^64 + r (0 <= r < m)
// let z = a*b = c*m + d (0 <= c, d < m)
// a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im
// c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2
// ((ab * im) >> 64) == c or c + 1
unsigned long long z = a;
z *= b;
#ifdef _MSC_VER
unsigned long long x;
_umul128(z, im, &x);
#else
unsigned long long x =
(unsigned long long)(((unsigned __int128)(z)*im) >> 64);
#endif
unsigned long long y = x * _m;
return (unsigned int)(z - y + (z < y ? _m : 0));
}
};
// @param n `0 <= n`
// @param m `1 <= m`
// @return `(x ** n) % m`
constexpr long long pow_mod_constexpr(long long x, long long n, int m) {
if (m == 1) return 0;
unsigned int _m = (unsigned int)(m);
unsigned long long r = 1;
unsigned long long y = safe_mod(x, m);
while (n) {
if (n & 1) r = (r * y) % _m;
y = (y * y) % _m;
n >>= 1;
}
return r;
}
// Reference:
// M. Forisek and J. Jancina,
// Fast Primality Testing for Integers That Fit into a Machine Word
// @param n `0 <= n`
constexpr bool is_prime_constexpr(int n) {
if (n <= 1) return false;
if (n == 2 || n == 7 || n == 61) return true;
if (n % 2 == 0) return false;
long long d = n - 1;
while (d % 2 == 0) d /= 2;
constexpr long long bases[3] = {2, 7, 61};
for (long long a : bases) {
long long t = d;
long long y = pow_mod_constexpr(a, t, n);
while (t != n - 1 && y != 1 && y != n - 1) {
y = y * y % n;
t <<= 1;
}
if (y != n - 1 && t % 2 == 0) {
return false;
}
}
return true;
}
template <int n> constexpr bool is_prime = is_prime_constexpr(n);
// @param b `1 <= b`
// @return pair(g, x) s.t. g = gcd(a, b), xa = g (mod b), 0 <= x < b/g
constexpr std::pair<long long, long long> inv_gcd(long long a, long long b) {
a = safe_mod(a, b);
if (a == 0) return {b, 0};
// Contracts:
// [1] s - m0 * a = 0 (mod b)
// [2] t - m1 * a = 0 (mod b)
// [3] s * |m1| + t * |m0| <= b
long long s = b, t = a;
long long m0 = 0, m1 = 1;
while (t) {
long long u = s / t;
s -= t * u;
m0 -= m1 * u; // |m1 * u| <= |m1| * s <= b
// [3]:
// (s - t * u) * |m1| + t * |m0 - m1 * u|
// <= s * |m1| - t * u * |m1| + t * (|m0| + |m1| * u)
// = s * |m1| + t * |m0| <= b
auto tmp = s;
s = t;
t = tmp;
tmp = m0;
m0 = m1;
m1 = tmp;
}
// by [3]: |m0| <= b/g
// by g != b: |m0| < b/g
if (m0 < 0) m0 += b / s;
return {s, m0};
}
// Compile time primitive root
// @param m must be prime
// @return primitive root (and minimum in now)
constexpr int primitive_root_constexpr(int m) {
if (m == 2) return 1;
if (m == 167772161) return 3;
if (m == 469762049) return 3;
if (m == 754974721) return 11;
if (m == 998244353) return 3;
int divs[20] = {};
divs[0] = 2;
int cnt = 1;
int x = (m - 1) / 2;
while (x % 2 == 0) x /= 2;
for (int i = 3; (long long)(i)*i <= x; i += 2) {
if (x % i == 0) {
divs[cnt++] = i;
while (x % i == 0) {
x /= i;
}
}
}
if (x > 1) {
divs[cnt++] = x;
}
for (int g = 2;; g++) {
bool ok = true;
for (int i = 0; i < cnt; i++) {
if (pow_mod_constexpr(g, (m - 1) / divs[i], m) == 1) {
ok = false;
break;
}
}
if (ok) return g;
}
}
template <int m> constexpr int primitive_root = primitive_root_constexpr(m);
// @param n `n < 2^32`
// @param m `1 <= m < 2^32`
// @return sum_{i=0}^{n-1} floor((ai + b) / m) (mod 2^64)
unsigned long long floor_sum_unsigned(unsigned long long n,
unsigned long long m,
unsigned long long a,
unsigned long long b) {
unsigned long long ans = 0;
while (true) {
if (a >= m) {
ans += n * (n - 1) / 2 * (a / m);
a %= m;
}
if (b >= m) {
ans += n * (b / m);
b %= m;
}
unsigned long long y_max = a * n + b;
if (y_max < m) break;
// y_max < m * (n + 1)
// floor(y_max / m) <= n
n = (unsigned long long)(y_max / m);
b = (unsigned long long)(y_max % m);
std::swap(m, a);
}
return ans;
}
} // namespace internal
} // namespace atcoder
namespace atcoder {
namespace internal {
#ifndef _MSC_VER
template <class T>
using is_signed_int128 =
typename std::conditional<std::is_same<T, __int128_t>::value ||
std::is_same<T, __int128>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_unsigned_int128 =
typename std::conditional<std::is_same<T, __uint128_t>::value ||
std::is_same<T, unsigned __int128>::value,
std::true_type,
std::false_type>::type;
template <class T>
using make_unsigned_int128 =
typename std::conditional<std::is_same<T, __int128_t>::value,
__uint128_t,
unsigned __int128>;
template <class T>
using is_integral = typename std::conditional<std::is_integral<T>::value ||
is_signed_int128<T>::value ||
is_unsigned_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_signed_int = typename std::conditional<(is_integral<T>::value &&
std::is_signed<T>::value) ||
is_signed_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_unsigned_int =
typename std::conditional<(is_integral<T>::value &&
std::is_unsigned<T>::value) ||
is_unsigned_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using to_unsigned = typename std::conditional<
is_signed_int128<T>::value,
make_unsigned_int128<T>,
typename std::conditional<std::is_signed<T>::value,
std::make_unsigned<T>,
std::common_type<T>>::type>::type;
#else
template <class T> using is_integral = typename std::is_integral<T>;
template <class T>
using is_signed_int =
typename std::conditional<is_integral<T>::value && std::is_signed<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_unsigned_int =
typename std::conditional<is_integral<T>::value &&
std::is_unsigned<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using to_unsigned = typename std::conditional<is_signed_int<T>::value,
std::make_unsigned<T>,
std::common_type<T>>::type;
#endif
template <class T>
using is_signed_int_t = std::enable_if_t<is_signed_int<T>::value>;
template <class T>
using is_unsigned_int_t = std::enable_if_t<is_unsigned_int<T>::value>;
template <class T> using to_unsigned_t = typename to_unsigned<T>::type;
} // namespace internal
} // namespace atcoder
namespace atcoder {
namespace internal {
struct modint_base {};
struct static_modint_base : modint_base {};
template <class T> using is_modint = std::is_base_of<modint_base, T>;
template <class T> using is_modint_t = std::enable_if_t<is_modint<T>::value>;
} // namespace internal
template <int m, std::enable_if_t<(1 <= m)>* = nullptr>
struct static_modint : internal::static_modint_base {
using mint = static_modint;
public:
static constexpr int mod() { return m; }
static mint raw(int v) {
mint x;
x._v = v;
return x;
}
static_modint() : _v(0) {}
template <class T, internal::is_signed_int_t<T>* = nullptr>
static_modint(T v) {
long long x = (long long)(v % (long long)(umod()));
if (x < 0) x += umod();
_v = (unsigned int)(x);
}
template <class T, internal::is_unsigned_int_t<T>* = nullptr>
static_modint(T v) {
_v = (unsigned int)(v % umod());
}
int val() const { return _v; }
mint& operator++() {
_v++;
if (_v == umod()) _v = 0;
return *this;
}
mint& operator--() {
if (_v == 0) _v = umod();
_v--;
return *this;
}
mint operator++(int) {
mint result = *this;
++*this;
return result;
}
mint operator--(int) {
mint result = *this;
--*this;
return result;
}
mint& operator+=(const mint& rhs) {
_v += rhs._v;
if (_v >= umod()) _v -= umod();
return *this;
}
mint& operator-=(const mint& rhs) {
_v -= rhs._v;
if (_v >= umod()) _v += umod();
return *this;
}
mint& operator*=(const mint& rhs) {
unsigned long long z = _v;
z *= rhs._v;
_v = (unsigned int)(z % umod());
return *this;
}
mint& operator/=(const mint& rhs) { return *this = *this * rhs.inv(); }
mint operator+() const { return *this; }
mint operator-() const { return mint() - *this; }
mint pow(long long n) const {
assert(0 <= n);
mint x = *this, r = 1;
while (n) {
if (n & 1) r *= x;
x *= x;
n >>= 1;
}
return r;
}
mint inv() const {
if (prime) {
assert(_v);
return pow(umod() - 2);
} else {
auto eg = internal::inv_gcd(_v, m);
assert(eg.first == 1);
return eg.second;
}
}
friend mint operator+(const mint& lhs, const mint& rhs) {
return mint(lhs) += rhs;
}
friend mint operator-(const mint& lhs, const mint& rhs) {
return mint(lhs) -= rhs;
}
friend mint operator*(const mint& lhs, const mint& rhs) {
return mint(lhs) *= rhs;
}
friend mint operator/(const mint& lhs, const mint& rhs) {
return mint(lhs) /= rhs;
}
friend bool operator==(const mint& lhs, const mint& rhs) {
return lhs._v == rhs._v;
}
friend bool operator!=(const mint& lhs, const mint& rhs) {
return lhs._v != rhs._v;
}
private:
unsigned int _v;
static constexpr unsigned int umod() { return m; }
static constexpr bool prime = internal::is_prime<m>;
};
template <int id> struct dynamic_modint : internal::modint_base {
using mint = dynamic_modint;
public:
static int mod() { return (int)(bt.umod()); }
static void set_mod(int m) {
assert(1 <= m);
bt = internal::barrett(m);
}
static mint raw(int v) {
mint x;
x._v = v;
return x;
}
dynamic_modint() : _v(0) {}
template <class T, internal::is_signed_int_t<T>* = nullptr>
dynamic_modint(T v) {
long long x = (long long)(v % (long long)(mod()));
if (x < 0) x += mod();
_v = (unsigned int)(x);
}
template <class T, internal::is_unsigned_int_t<T>* = nullptr>
dynamic_modint(T v) {
_v = (unsigned int)(v % mod());
}
int val() const { return _v; }
mint& operator++() {
_v++;
if (_v == umod()) _v = 0;
return *this;
}
mint& operator--() {
if (_v == 0) _v = umod();
_v--;
return *this;
}
mint operator++(int) {
mint result = *this;
++*this;
return result;
}
mint operator--(int) {
mint result = *this;
--*this;
return result;
}
mint& operator+=(const mint& rhs) {
_v += rhs._v;
if (_v >= umod()) _v -= umod();
return *this;
}
mint& operator-=(const mint& rhs) {
_v += mod() - rhs._v;
if (_v >= umod()) _v -= umod();
return *this;
}
mint& operator*=(const mint& rhs) {
_v = bt.mul(_v, rhs._v);
return *this;
}
mint& operator/=(const mint& rhs) { return *this = *this * rhs.inv(); }
mint operator+() const { return *this; }
mint operator-() const { return mint() - *this; }
mint pow(long long n) const {
assert(0 <= n);
mint x = *this, r = 1;
while (n) {
if (n & 1) r *= x;
x *= x;
n >>= 1;
}
return r;
}
mint inv() const {
auto eg = internal::inv_gcd(_v, mod());
assert(eg.first == 1);
return eg.second;
}
friend mint operator+(const mint& lhs, const mint& rhs) {
return mint(lhs) += rhs;
}
friend mint operator-(const mint& lhs, const mint& rhs) {
return mint(lhs) -= rhs;
}
friend mint operator*(const mint& lhs, const mint& rhs) {
return mint(lhs) *= rhs;
}
friend mint operator/(const mint& lhs, const mint& rhs) {
return mint(lhs) /= rhs;
}
friend bool operator==(const mint& lhs, const mint& rhs) {
return lhs._v == rhs._v;
}
friend bool operator!=(const mint& lhs, const mint& rhs) {
return lhs._v != rhs._v;
}
private:
unsigned int _v;
static internal::barrett bt;
static unsigned int umod() { return bt.umod(); }
};
template <int id> internal::barrett dynamic_modint<id>::bt(998244353);
using modint998244353 = static_modint<998244353>;
using modint1000000007 = static_modint<1000000007>;
using modint = dynamic_modint<-1>;
namespace internal {
template <class T>
using is_static_modint = std::is_base_of<internal::static_modint_base, T>;
template <class T>
using is_static_modint_t = std::enable_if_t<is_static_modint<T>::value>;
template <class> struct is_dynamic_modint : public std::false_type {};
template <int id>
struct is_dynamic_modint<dynamic_modint<id>> : public std::true_type {};
template <class T>
using is_dynamic_modint_t = std::enable_if_t<is_dynamic_modint<T>::value>;
} // namespace internal
} // namespace atcoder
namespace atcoder {
namespace internal {
template <class mint,
int g = internal::primitive_root<mint::mod()>,
internal::is_static_modint_t<mint>* = nullptr>
struct fft_info {
static constexpr int rank2 = countr_zero_constexpr(mint::mod() - 1);
std::array<mint, rank2 + 1> root; // root[i]^(2^i) == 1
std::array<mint, rank2 + 1> iroot; // root[i] * iroot[i] == 1
std::array<mint, std::max(0, rank2 - 2 + 1)> rate2;
std::array<mint, std::max(0, rank2 - 2 + 1)> irate2;
std::array<mint, std::max(0, rank2 - 3 + 1)> rate3;
std::array<mint, std::max(0, rank2 - 3 + 1)> irate3;
fft_info() {
root[rank2] = mint(g).pow((mint::mod() - 1) >> rank2);
iroot[rank2] = root[rank2].inv();
for (int i = rank2 - 1; i >= 0; i--) {
root[i] = root[i + 1] * root[i + 1];
iroot[i] = iroot[i + 1] * iroot[i + 1];
}
{
mint prod = 1, iprod = 1;
for (int i = 0; i <= rank2 - 2; i++) {
rate2[i] = root[i + 2] * prod;
irate2[i] = iroot[i + 2] * iprod;
prod *= iroot[i + 2];
iprod *= root[i + 2];
}
}
{
mint prod = 1, iprod = 1;
for (int i = 0; i <= rank2 - 3; i++) {
rate3[i] = root[i + 3] * prod;
irate3[i] = iroot[i + 3] * iprod;
prod *= iroot[i + 3];
iprod *= root[i + 3];
}
}
}
};
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
void butterfly(std::vector<mint>& a) {
int n = int(a.size());
int h = internal::countr_zero((unsigned int)n);
static const fft_info<mint> info;
int len = 0; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed
while (len < h) {
if (h - len == 1) {
int p = 1 << (h - len - 1);
mint rot = 1;
for (int s = 0; s < (1 << len); s++) {
int offset = s << (h - len);
for (int i = 0; i < p; i++) {
auto l = a[i + offset];
auto r = a[i + offset + p] * rot;
a[i + offset] = l + r;
a[i + offset + p] = l - r;
}
if (s + 1 != (1 << len))
rot *= info.rate2[countr_zero(~(unsigned int)(s))];
}
len++;
} else {
// 4-base
int p = 1 << (h - len - 2);
mint rot = 1, imag = info.root[2];
for (int s = 0; s < (1 << len); s++) {
mint rot2 = rot * rot;
mint rot3 = rot2 * rot;
int offset = s << (h - len);
for (int i = 0; i < p; i++) {
auto mod2 = 1ULL * mint::mod() * mint::mod();
auto a0 = 1ULL * a[i + offset].val();
auto a1 = 1ULL * a[i + offset + p].val() * rot.val();
auto a2 = 1ULL * a[i + offset + 2 * p].val() * rot2.val();
auto a3 = 1ULL * a[i + offset + 3 * p].val() * rot3.val();
auto a1na3imag =
1ULL * mint(a1 + mod2 - a3).val() * imag.val();
auto na2 = mod2 - a2;
a[i + offset] = a0 + a2 + a1 + a3;
a[i + offset + 1 * p] = a0 + a2 + (2 * mod2 - (a1 + a3));
a[i + offset + 2 * p] = a0 + na2 + a1na3imag;
a[i + offset + 3 * p] = a0 + na2 + (mod2 - a1na3imag);
}
if (s + 1 != (1 << len))
rot *= info.rate3[countr_zero(~(unsigned int)(s))];
}
len += 2;
}
}
}
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
void butterfly_inv(std::vector<mint>& a) {
int n = int(a.size());
int h = internal::countr_zero((unsigned int)n);
static const fft_info<mint> info;
int len = h; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed
while (len) {
if (len == 1) {
int p = 1 << (h - len);
mint irot = 1;
for (int s = 0; s < (1 << (len - 1)); s++) {
int offset = s << (h - len + 1);
for (int i = 0; i < p; i++) {
auto l = a[i + offset];
auto r = a[i + offset + p];
a[i + offset] = l + r;
a[i + offset + p] =
(unsigned long long)((unsigned int)(l.val() - r.val()) + mint::mod()) *
irot.val();
;
}
if (s + 1 != (1 << (len - 1)))
irot *= info.irate2[countr_zero(~(unsigned int)(s))];
}
len--;
} else {
// 4-base
int p = 1 << (h - len);
mint irot = 1, iimag = info.iroot[2];
for (int s = 0; s < (1 << (len - 2)); s++) {
mint irot2 = irot * irot;
mint irot3 = irot2 * irot;
int offset = s << (h - len + 2);
for (int i = 0; i < p; i++) {
auto a0 = 1ULL * a[i + offset + 0 * p].val();
auto a1 = 1ULL * a[i + offset + 1 * p].val();
auto a2 = 1ULL * a[i + offset + 2 * p].val();
auto a3 = 1ULL * a[i + offset + 3 * p].val();
auto a2na3iimag =
1ULL *
mint((mint::mod() + a2 - a3) * iimag.val()).val();
a[i + offset] = a0 + a1 + a2 + a3;
a[i + offset + 1 * p] =
(a0 + (mint::mod() - a1) + a2na3iimag) * irot.val();
a[i + offset + 2 * p] =
(a0 + a1 + (mint::mod() - a2) + (mint::mod() - a3)) *
irot2.val();
a[i + offset + 3 * p] =
(a0 + (mint::mod() - a1) + (mint::mod() - a2na3iimag)) *
irot3.val();
}
if (s + 1 != (1 << (len - 2)))
irot *= info.irate3[countr_zero(~(unsigned int)(s))];
}
len -= 2;
}
}
}
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
std::vector<mint> convolution_naive(const std::vector<mint>& a,
const std::vector<mint>& b) {
int n = int(a.size()), m = int(b.size());
std::vector<mint> ans(n + m - 1);
if (n < m) {
for (int j = 0; j < m; j++) {
for (int i = 0; i < n; i++) {
ans[i + j] += a[i] * b[j];
}
}
} else {
for (int i = 0; i < n; i++) {
for (int j = 0; j < m; j++) {
ans[i + j] += a[i] * b[j];
}
}
}
return ans;
}
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
std::vector<mint> convolution_fft(std::vector<mint> a, std::vector<mint> b) {
int n = int(a.size()), m = int(b.size());
int z = (int)internal::bit_ceil((unsigned int)(n + m - 1));
a.resize(z);
internal::butterfly(a);
b.resize(z);
internal::butterfly(b);
for (int i = 0; i < z; i++) {
a[i] *= b[i];
}
internal::butterfly_inv(a);
a.resize(n + m - 1);
mint iz = mint(z).inv();
for (int i = 0; i < n + m - 1; i++) a[i] *= iz;
return a;
}
} // namespace internal
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
std::vector<mint> convolution(std::vector<mint>&& a, std::vector<mint>&& b) {
int n = int(a.size()), m = int(b.size());
if (!n || !m) return {};
int z = (int)internal::bit_ceil((unsigned int)(n + m - 1));
assert((mint::mod() - 1) % z == 0);
if (std::min(n, m) <= 60) return convolution_naive(std::move(a), std::move(b));
return internal::convolution_fft(std::move(a), std::move(b));
}
template <class mint, internal::is_static_modint_t<mint>* = nullptr>
std::vector<mint> convolution(const std::vector<mint>& a,
const std::vector<mint>& b) {
int n = int(a.size()), m = int(b.size());
if (!n || !m) return {};
int z = (int)internal::bit_ceil((unsigned int)(n + m - 1));
assert((mint::mod() - 1) % z == 0);
if (std::min(n, m) <= 60) return convolution_naive(a, b);
return internal::convolution_fft(a, b);
}
template <unsigned int mod = 998244353,
class T,
std::enable_if_t<internal::is_integral<T>::value>* = nullptr>
std::vector<T> convolution(const std::vector<T>& a, const std::vector<T>& b) {
int n = int(a.size()), m = int(b.size());
if (!n || !m) return {};
using mint = static_modint<mod>;
int z = (int)internal::bit_ceil((unsigned int)(n + m - 1));
assert((mint::mod() - 1) % z == 0);
std::vector<mint> a2(n), b2(m);
for (int i = 0; i < n; i++) {
a2[i] = mint(a[i]);
}
for (int i = 0; i < m; i++) {
b2[i] = mint(b[i]);
}
auto c2 = convolution(std::move(a2), std::move(b2));
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
c[i] = c2[i].val();
}
return c;
}
std::vector<long long> convolution_ll(const std::vector<long long>& a,
const std::vector<long long>& b) {
int n = int(a.size()), m = int(b.size());
if (!n || !m) return {};
static constexpr unsigned long long MOD1 = 754974721; // 2^24
static constexpr unsigned long long MOD2 = 167772161; // 2^25
static constexpr unsigned long long MOD3 = 469762049; // 2^26
static constexpr unsigned long long M2M3 = MOD2 * MOD3;
static constexpr unsigned long long M1M3 = MOD1 * MOD3;
static constexpr unsigned long long M1M2 = MOD1 * MOD2;
static constexpr unsigned long long M1M2M3 = MOD1 * MOD2 * MOD3;
static constexpr unsigned long long i1 =
internal::inv_gcd(MOD2 * MOD3, MOD1).second;
static constexpr unsigned long long i2 =
internal::inv_gcd(MOD1 * MOD3, MOD2).second;
static constexpr unsigned long long i3 =
internal::inv_gcd(MOD1 * MOD2, MOD3).second;
static constexpr int MAX_AB_BIT = 24;
static_assert(MOD1 % (1ull << MAX_AB_BIT) == 1, "MOD1 isn't enough to support an array length of 2^24.");
static_assert(MOD2 % (1ull << MAX_AB_BIT) == 1, "MOD2 isn't enough to support an array length of 2^24.");
static_assert(MOD3 % (1ull << MAX_AB_BIT) == 1, "MOD3 isn't enough to support an array length of 2^24.");
assert(n + m - 1 <= (1 << MAX_AB_BIT));
auto c1 = convolution<MOD1>(a, b);
auto c2 = convolution<MOD2>(a, b);
auto c3 = convolution<MOD3>(a, b);
std::vector<long long> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++) {
unsigned long long x = 0;
x += (c1[i] * i1) % MOD1 * M2M3;
x += (c2[i] * i2) % MOD2 * M1M3;
x += (c3[i] * i3) % MOD3 * M1M2;
// B = 2^63, -B <= x, r(real value) < B
// (x, x - M, x - 2M, or x - 3M) = r (mod 2B)
// r = c1[i] (mod MOD1)
// focus on MOD1
// r = x, x - M', x - 2M', x - 3M' (M' = M % 2^64) (mod 2B)
// r = x,
// x - M' + (0 or 2B),
// x - 2M' + (0, 2B or 4B),
// x - 3M' + (0, 2B, 4B or 6B) (without mod!)
// (r - x) = 0, (0)
// - M' + (0 or 2B), (1)
// -2M' + (0 or 2B or 4B), (2)
// -3M' + (0 or 2B or 4B or 6B) (3) (mod MOD1)
// we checked that
// ((1) mod MOD1) mod 5 = 2
// ((2) mod MOD1) mod 5 = 3
// ((3) mod MOD1) mod 5 = 4
long long diff =
c1[i] - internal::safe_mod((long long)(x), (long long)(MOD1));
if (diff < 0) diff += MOD1;
static constexpr unsigned long long offset[5] = {
0, 0, M1M2M3, 2 * M1M2M3, 3 * M1M2M3};
x -= offset[diff % 5];
c[i] = x;
}
return c;
}
} // namespace atcoder
namespace noya {
inline int highest_bit(unsigned x) {
return x == 0 ? -1 : 31 - __builtin_clz(x);
}
/// @brief Compute bitwise-AND convolution of two vectors of length 2^e.
template <class T>
std::vector<T> bitwise_and_convolution(std::vector<T> a, std::vector<T> b) {
int N = int(a.size());
assert(int(b.size()) == N);
int e = highest_bit(N);
assert(N == (1 << e));
auto convolution = [&](std::vector<T> &f) -> void {
for (int j = 0; j < e; j++)
for (int i = 0; i < N; i++)
if (i >> j & 1)
f[i ^ 1 << j] += f[i];
};
auto envolution = [&](std::vector<T> &f) -> void {
for (int j = 0; j < e; j++)
for (int i = 0; i < N; i++)
if (i >> j & 1)
f[i ^ 1 << j] -= f[i];
};
convolution(a), convolution(b);
std::vector<T> c(N);
for (int i = 0; i < N; i++)
c[i] = a[i] * b[i];
envolution(c);
return c;
}
/// @brief Compute bitwise-OR convolution of two vectors of length 2^e.
template <class T>
std::vector<T> bitwise_or_convolution(std::vector<T> a, std::vector<T> b) {
int N = int(a.size());
assert(int(b.size()) == N);
int e = highest_bit(N);
assert(N == (1 << e));
std::reverse(a.begin(), a.end());
std::reverse(b.begin(), b.end());
auto c = bitwise_and_convolution(a, b);
std::reverse(c.begin(), c.end());
return c;
}
/// @brief Compute bitwise-XOR convolution of two vectors of length 2^e.
template <class T>
std::vector<T> bitwise_xor_convolution(std::vector<T> a, std::vector<T> b) {
int N = int(a.size());
assert(int(b.size()) == N);
int e = highest_bit(N);
assert(N == (1 << e));
auto convolution = [&](std::vector<T> &f) -> void {
for (int j = 0; j < e; j++)
for (int i = 0; i < N; i++)
if (i >> j & 1) {
T x = f[i], y = f[i ^ (1 << j)];
f[i] = y - x;
f[i ^ (1 << j)] = x + y;
}
};
auto envolution = [&](std::vector<T> &f) -> void {
convolution(f);
T inv = T(N).inv();
for (int i = 0; i < N; i++)
f[i] *= inv;
};
convolution(a), convolution(b);
std::vector<T> c(N);
for (int i = 0; i < N; i++)
c[i] = a[i] * b[i];
envolution(c);
return c;
}
/// @brief Compute c[k] = sum_{gcd(i,j)=k} a[i] b[j]. Index zero is ignored.
template <class T>
std::vector<T> gcd_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
assert(a.size() == b.size());
if (a.empty()) {
return {};
}
int n = int(a.size()) - 1;
std::vector<T> transformed_a(n + 1), transformed_b(n + 1);
for (int divisor = 1; divisor <= n; divisor++) {
for (int multiple = divisor; multiple <= n; multiple += divisor) {
transformed_a[divisor] += a[multiple];
transformed_b[divisor] += b[multiple];
}
}
std::vector<T> result(n + 1);
for (int i = 1; i <= n; i++) {
result[i] = transformed_a[i] * transformed_b[i];
}
for (int divisor = n; divisor >= 1; divisor--) {
for (int multiple = divisor * 2; multiple <= n; multiple += divisor) {
result[divisor] -= result[multiple];
}
}
return result;
}
/// @brief Compute c[k] = sum_{lcm(i,j)=k} a[i] b[j]. Index zero is ignored.
template <class T>
std::vector<T> lcm_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
assert(a.size() == b.size());
if (a.empty()) {
return {};
}
int n = int(a.size()) - 1;
std::vector<T> transformed_a(n + 1), transformed_b(n + 1);
for (int divisor = 1; divisor <= n; divisor++) {
for (int multiple = divisor; multiple <= n; multiple += divisor) {
transformed_a[multiple] += a[divisor];
transformed_b[multiple] += b[divisor];
}
}
std::vector<T> result(n + 1);
for (int i = 1; i <= n; i++) {
result[i] = transformed_a[i] * transformed_b[i];
}
for (int divisor = 1; divisor <= n; divisor++) {
for (int multiple = divisor * 2; multiple <= n; multiple += divisor) {
result[multiple] -= result[divisor];
}
}
return result;
}
/// @brief Compute subset convolution of two vectors of length 2^e.
template <class T>
std::vector<T> subset_convolution(std::vector<T> a, std::vector<T> b) {
int N = int(a.size());
assert(int(b.size()) == N);
int e = highest_bit(N);
assert(N == (1 << e));
auto convolution = [&](std::vector<T> &f) -> void {
for (int j = 0; j < e; j++)
for (int i = 0; i < N; i++)
if (~i >> j & 1)
f[i ^ 1 << j] += f[i];
};
auto envolution = [&](std::vector<T> &f) -> void {
for (int j = 0; j < e; j++)
for (int i = 0; i < N; i++)
if (~i >> j & 1)
f[i ^ 1 << j] -= f[i];
};
std::vector<std::vector<T>> f(e + 1, std::vector<T>(N));
std::vector<std::vector<T>> g(e + 1, std::vector<T>(N));
std::vector<T> ans(N);
for (int i = 0; i < N; i++) {
f[__builtin_popcount(i)][i] = a[i];
g[__builtin_popcount(i)][i] = b[i];
}
for (int i = 0; i < e; i++) {
convolution(f[i]);
convolution(g[i]);
}
std::vector<T> tmp(N);
for (int i = 0; i <= e; i++) {
for (int j = 0; j < N; j++)
tmp[j] = 0;
for (int j = 0; j <= i; j++)
for (int k = 0; k < N; k++)
tmp[k] += f[j][k] * g[i - j][k];
envolution(tmp);
for (int k = 0; k < N; k++)
if (__builtin_popcount(k) == i)
ans[k] = tmp[k];
}
return ans;
}
} // namespace noya
namespace noya {
namespace set_power_series_internal {
template <class T, class Operation>
void subset_transform(std::vector<T> &values, Operation operation) {
int size = int(values.size());
for (int width = 1; width < size; width <<= 1) {
for (int block = 0; block < size; block += width << 1) {
for (int offset = 0; offset < width; offset++) {
operation(values[block + offset],
values[block + width + offset]);
}
}
}
}
template <class Mint>
std::vector<Mint>
transposed_subset_product(const std::vector<Mint> &series,
std::vector<Mint> weights) {
std::reverse(weights.begin(), weights.end());
weights = subset_convolution(weights, series);
std::reverse(weights.begin(), weights.end());
return weights;
}
template <class Mint>
std::vector<Mint>
transposed_egf_composition(const std::vector<Mint> &series,
const std::vector<Mint> &weights) {
int size = int(series.size());
assert(size > 0 && (size & (size - 1)) == 0);
assert(weights.size() == series.size());
assert(series[0] == Mint{});
int variables = size == 1 ? 0 : 31 - __builtin_clz(size);
std::vector<Mint> result(variables + 1);
result[0] = weights[0];
std::vector<Mint> state = weights;
for (int degree = 0; degree < variables; degree++) {
std::vector<Mint> next(1 << (variables - degree - 1));
for (int bit = 0; bit < variables - degree; bit++) {
int length = 1 << bit;
std::vector<Mint> block_series(series.begin() + length,
series.begin() + 2 * length);
std::vector<Mint> block_weights(state.begin() + length,
state.begin() + 2 * length);
block_weights =
transposed_subset_product(block_series, std::move(block_weights));
for (int mask = 0; mask < length; mask++) {
next[mask] += block_weights[mask];
}
}
state = std::move(next);
result[degree + 1] = state[0];
}
return result;
}
} // namespace set_power_series_internal
/// @brief Compose an exponential generating function with a set power series.
/// A coefficient indexed by a mask represents the square-free monomial whose
/// variables are that mask. Splitting masks by their largest variable and by
/// popcount turns every multiplication into ranked subset convolution. The
/// subset zeta transform performs all disjoint decompositions at once, while
/// the rank selects exactly the square-free terms. `egf[k]` is the multiplier
/// of `series^k / k!`, and the constant term of `series` must be zero.
template <class Mint>
std::vector<Mint>
set_power_series_egf_composition(const std::vector<Mint> &egf,
const std::vector<Mint> &series) {
assert(!series.empty() && (series.size() & (series.size() - 1)) == 0);
int variables = series.size() == 1
? 0
: 31 - __builtin_clz(unsigned(series.size()));
assert(int(egf.size()) == variables + 1);
assert(series[0] == Mint{});
std::vector<std::vector<std::vector<Mint>>> states(variables + 1);
for (int level = 0; level <= variables; level++) {
states[level].assign(level + 1,
std::vector<Mint>(std::size_t(1) << level));
states[level][0][0] = egf[variables - level];
}
for (int bit = 0; bit < variables; bit++) {
int half = 1 << bit;
std::vector<std::vector<Mint>> ranked_series(
bit + 1, std::vector<Mint>(half));
for (int mask = 0; mask < half; mask++) {
ranked_series[__builtin_popcount(unsigned(mask))][mask] =
series[half + mask];
}
for (auto &rank : ranked_series) {
set_power_series_internal::subset_transform(
rank, [](Mint &lower, Mint &upper) { upper += lower; });
}
for (int level = bit + 1; level <= variables; level++) {
for (int rank = 0; rank <= level; rank++) {
std::copy_n(states[level][rank].begin(), half,
states[level][rank].begin() + half);
}
}
for (int level = bit; level < variables; level++) {
for (int series_rank = 0; series_rank <= bit; series_rank++) {
for (int state_rank = 0;
series_rank + state_rank <= level; state_rank++) {
for (int mask = 0; mask < half; mask++) {
states[level + 1][series_rank + state_rank + 1][half + mask] +=
ranked_series[series_rank][mask] *
states[level][state_rank][mask];
}
}
}
}
}
for (auto &rank : states[variables]) {
set_power_series_internal::subset_transform(
rank, [](Mint &lower, Mint &upper) { upper -= lower; });
}
std::vector<Mint> result(series.size());
for (int mask = 0; mask < int(series.size()); mask++) {
result[mask] =
states[variables][__builtin_popcount(unsigned(mask))][mask];
}
return result;
}
/// @brief Substitute a set power series into an ordinary polynomial.
/// Taylor expansion at the constant coefficient `c=series[0]` gives the EGF
/// multipliers `f^(k)(c)`. Only the first n derivatives matter because every
/// square-free monomial has degree at most n; the ranked composition routine
/// then evaluates the nilpotent remainder `series-c`.
template <class Mint>
std::vector<Mint>
set_power_series_polynomial_composition(std::vector<Mint> polynomial,
const std::vector<Mint> &series) {
assert(!series.empty() && (series.size() & (series.size() - 1)) == 0);
int variables = series.size() == 1
? 0
: 31 - __builtin_clz(unsigned(series.size()));
std::vector<Mint> derivatives(variables + 1);
Mint constant = series[0];
for (int order = 0; order <= variables && !polynomial.empty(); order++) {
Mint value{};
for (int index = int(polynomial.size()) - 1; index >= 0; index--) {
value = value * constant + polynomial[index];
}
derivatives[order] = value;
for (int index = 1; index < int(polynomial.size()); index++) {
polynomial[index - 1] = polynomial[index] * Mint(index);
}
polynomial.pop_back();
}
std::vector<Mint> nilpotent = series;
nilpotent[0] = Mint{};
return set_power_series_egf_composition(derivatives, nilpotent);
}
/// @brief Return the exponential of a set power series with zero constant.
/// In the square-free quotient every term above degree n vanishes, so exp is
/// the finite EGF with every multiplier equal to one.
template <class Mint>
std::vector<Mint>
set_power_series_exponential(const std::vector<Mint> &series) {
assert(!series.empty() && (series.size() & (series.size() - 1)) == 0);
int variables = series.size() == 1
? 0
: 31 - __builtin_clz(unsigned(series.size()));
assert(series[0] == Mint{});
return set_power_series_egf_composition(
std::vector<Mint>(variables + 1, Mint(1)), series);
}
/// @brief Return the logarithm of a set power series with constant one.
/// Writing the input as `1+u`, the nilpotent series u satisfies
/// log(1+u)=sum (-1)^(k-1)u^k/k. Therefore its kth EGF multiplier is
/// `(-1)^(k-1)(k-1)!`, and ranked EGF composition evaluates all terms.
template <class Mint>
std::vector<Mint>
set_power_series_logarithm(const std::vector<Mint> &series) {
assert(!series.empty() && (series.size() & (series.size() - 1)) == 0);
int variables = series.size() == 1
? 0
: 31 - __builtin_clz(unsigned(series.size()));
assert(series[0] == Mint(1));
std::vector<Mint> multipliers(variables + 1);
Mint factorial = Mint(1);
for (int degree = 1; degree <= variables; degree++) {
multipliers[degree] = (degree & 1) ? factorial : -factorial;
factorial *= Mint(degree);
}
std::vector<Mint> nilpotent = series;
nilpotent[0] -= Mint(1);
return set_power_series_egf_composition(multipliers, nilpotent);
}
/// @brief Project consecutive powers of a set power series onto weights.
/// Reversing masks changes the transpose of subset convolution into an
/// ordinary subset convolution. Repeated transposed ranked composition gives
/// the projections of `u^k/k!` for `u=series-series[0]`; convolving those n+1
/// values with `c^j/j!` and multiplying by m! restores `(c+u)^m`.
template <class Mint>
std::vector<Mint>
set_power_series_power_projection(const std::vector<Mint> &series,
const std::vector<Mint> &weights,
int count) {
assert(count >= 0);
assert(!series.empty() && (series.size() & (series.size() - 1)) == 0);
assert(weights.size() == series.size());
if (count == 0) {
return {};
}
Mint constant = series[0];
std::vector<Mint> nilpotent = series;
nilpotent[0] = Mint{};
std::vector<Mint> projected =
set_power_series_internal::transposed_egf_composition(nilpotent,
weights);
std::vector<Mint> factorial(count, Mint(1));
std::vector<Mint> inverse_factorial(count, Mint(1));
for (int index = 1; index < count; index++) {
factorial[index] = factorial[index - 1] * Mint(index);
}
inverse_factorial.back() = Mint(1) / factorial.back();
for (int index = count - 1; index > 0; index--) {
inverse_factorial[index - 1] = inverse_factorial[index] * Mint(index);
}
std::vector<Mint> constant_egf(count);
Mint power = Mint(1);
for (int index = 0; index < count; index++) {
constant_egf[index] = power * inverse_factorial[index];
power *= constant;
}
std::vector<Mint> result(count);
for (int degree = 0; degree < int(projected.size()); degree++) {
for (int total = degree; total < count; total++) {
result[total] += projected[degree] * constant_egf[total - degree];
}
}
for (int index = 0; index < count; index++) {
result[index] *= factorial[index];
}
return result;
}
} // namespace noya