library/setpowerseries/Compose.hpp
Depends on
Verified with
Code
// https://maspypy.com/%e9%9b%86%e5%90%88%e3%81%b9%e3%81%8d%e7%b4%9a%e6%95%b0%e9%96%a2%e9%80%a3-2-%e5%a4%9a%e9%a0%85%e5%bc%8f%e3%81%a8%e3%81%ae%e5%90%88%e6%88%90
#include "library/setpowerseries/Base.hpp"
#include <vector>
template <typename SPS, typename T = typename SPS::value_type>
SPS SPS_exp_comp(const std::vector<T> &f, const SPS &a) {
int N = a.size();
int n = bitwise::log2(N);
std::vector<SPS> ret(n + 1);
for (int d = n; d >= 0; d--) {
ret[d].resize(1 << (n - d));
ret[d][0] = f[d];
for (int m = 0; m < n - d; m++) {
// ret[d] の [2^m, 2^{m+1}] を求める
SPS pre(ret[d + 1].begin(), ret[d + 1].begin() + (1 << m));
SPS a2(a.begin() + (1 << m), a.begin() + (1 << (m + 1)));
pre *= a2;
std::ranges::copy(pre, ret[d].begin() + (1 << m));
}
}
return ret[0];
}
// sum_k f_k a^k
template <typename SPS, typename T = typename SPS::value_type>
SPS SPS_composition(const std::vector<T> &f, SPS a) {
int N = a.size();
int n = bitwise::log2(N);
std::vector<T> c_pow(f.size() + 1, 1);
for (int i = 1; i <= f.size(); i++)
c_pow[i] = c_pow[i - 1] * a[0];
std::vector<T> f2(n + 1, 0);
for (int i = 0; i < f.size(); i++) {
T now = f[i];
for (int k = 0; k <= std::min(i, n); k++) {
f2[k] += now * c_pow[i - k];
now *= i - k;
}
}
a[0] = 0;
return SPS_exp_comp(f2, a);
}
template <typename SPS, typename T = typename SPS::value_type> SPS exp(SPS a) {
int N = a.size();
int n = bitwise::log2(N);
SPS ret(N, 1);
for (int d = n - 1; d >= 0; d--) {
const int M = 1 << (n - d - 1);
SPS pre(ret.begin(), ret.begin() + M);
SPS a2(a.begin() + M, a.begin() + 2 * M);
pre *= a2;
std::ranges::copy(pre, ret.begin() + M);
}
return ret;
}
template <typename SPS, typename T = typename SPS::value_type> SPS log(SPS a) {
int n = bitwise::log2(a.size());
assert((1 << n) == a.size() and a[0] == 1);
std::vector<T> lg(n + 1, 0);
for (int k = 1; k <= n; k++)
lg[k] = T(-1) / k;
return SPS_composition(lg, SPS(1, 1) - a);
}
#line 1 "library/setpowerseries/Compose.hpp"
// https://maspypy.com/%e9%9b%86%e5%90%88%e3%81%b9%e3%81%8d%e7%b4%9a%e6%95%b0%e9%96%a2%e9%80%a3-2-%e5%a4%9a%e9%a0%85%e5%bc%8f%e3%81%a8%e3%81%ae%e5%90%88%e6%88%90
#line 2 "library/bitwise/Util.hpp"
namespace bitwise{
static int log2(int N){
int n=__builtin_ffs(N)-1;
assert((1<<n)==N);
return n;
}
static bool in(int S,int a){ return (S>>a)&1; }
}
#line 3 "library/bitwise/Ranked.hpp"
#define REP_(i, n) for (int i = 0; i < (n); i++)
#define RREP_(i, n) for (int i = (n) - 1; i >= 0; i--)
class BitwiseRanked {
static int popcount(int S) { return __builtin_popcount(S); }
public:
template <typename T>
static std::vector<std::vector<T>> zeta(const std::vector<T> &A) {
const int n = bitwise::log2(A.size());
std::vector<std::vector<T>> RA(1 << n, std::vector<T>(n + 1, 0));
REP_(S, 1 << n) RA[S][popcount(S)] = A[S];
REP_(i, n)
REP_(S, 1 << n)
if (!bitwise::in(S, i))
REP_(d, n + 1) RA[S | (1 << i)][d] += RA[S][d];
return RA;
}
template <typename T>
static std::vector<T> mobius(std::vector<std::vector<T>> RA) {
const int n = bitwise::log2(RA.size());
REP_(i, n)
REP_(S, 1 << n)
if (!bitwise::in(S, i))
REP_(d, n + 1) RA[S | (1 << i)][d] -= RA[S][d];
std::vector<T> A(1 << n);
REP_(S, 1 << n) A[S] = RA[S][popcount(S)];
return A;
}
template <typename T>
static std::vector<T> convolution(const std::vector<T> &A,
const std::vector<T> &B) {
const int n = bitwise::log2(A.size());
auto RA = zeta(A);
auto RB = zeta(B);
REP_(S, 1 << n) {
auto &ra = RA[S], rb = RB[S];
RREP_(d, n + 1) {
ra[d] *= rb[0];
REP_(i, d) ra[d] += ra[i] * rb[d - i];
}
}
return mobius(RA);
}
};
#undef REP_
#undef RREP_
#line 1 "library/util/Valarray.hpp"
#include <functional>
#include <ranges>
#include <vector>
template <typename T> struct Valarray : std::vector<T> {
using std::vector<T>::vector; // コンストラクタ継承
Valarray(const std::vector<T> &v) : std::vector<T>(v.begin(), v.end()) {}
private:
template <typename Op>
Valarray &apply_inplace(const Valarray &other, Op op) {
if (this->size() < other.size())
this->resize(other.size(), T(0));
for (auto [a, b] : std::views::zip(*this, other))
a = op(a, b);
return *this;
}
public:
Valarray &operator+=(const Valarray &other) {
return apply_inplace(other, std::plus<>());
}
Valarray &operator-=(const Valarray &other) {
return apply_inplace(other, std::minus<>());
}
Valarray &operator*=(const Valarray &other) {
return apply_inplace(other, std::multiplies<>());
}
Valarray &operator/=(const Valarray &other) {
return apply_inplace(other, std::divides<>());
}
friend Valarray operator+(Valarray a, const Valarray &b) { return a += b; }
friend Valarray operator-(Valarray a, const Valarray &b) { return a -= b; }
friend Valarray operator*(Valarray a, const Valarray &b) { return a *= b; }
friend Valarray operator/(Valarray a, const Valarray &b) { return a /= b; }
Valarray operator-() const {
Valarray g = *this;
for (T &a : g)
a = -a;
return g;
}
};
#line 4 "library/setpowerseries/Base.hpp"
template <typename T> struct SetPowerSeries : Valarray<T> {
using SPS = SetPowerSeries;
using Valarray<T>::Valarray;
using Valarray<T>::size;
using Valarray<T>::at;
using value_type = T;
SetPowerSeries(const std::vector<T> &f) : Valarray<T>(f) {}
SPS operator-() const { return SPS(Valarray<T>::operator-()); }
SPS operator*(const SPS &b) const {
return SPS(BitwiseRanked::convolution<T>(*this, b));
}
SPS &operator*=(const SPS &b) { return (*this) = (*this) * b; }
};
#line 4 "library/setpowerseries/Compose.hpp"
template <typename SPS, typename T = typename SPS::value_type>
SPS SPS_exp_comp(const std::vector<T> &f, const SPS &a) {
int N = a.size();
int n = bitwise::log2(N);
std::vector<SPS> ret(n + 1);
for (int d = n; d >= 0; d--) {
ret[d].resize(1 << (n - d));
ret[d][0] = f[d];
for (int m = 0; m < n - d; m++) {
// ret[d] の [2^m, 2^{m+1}] を求める
SPS pre(ret[d + 1].begin(), ret[d + 1].begin() + (1 << m));
SPS a2(a.begin() + (1 << m), a.begin() + (1 << (m + 1)));
pre *= a2;
std::ranges::copy(pre, ret[d].begin() + (1 << m));
}
}
return ret[0];
}
// sum_k f_k a^k
template <typename SPS, typename T = typename SPS::value_type>
SPS SPS_composition(const std::vector<T> &f, SPS a) {
int N = a.size();
int n = bitwise::log2(N);
std::vector<T> c_pow(f.size() + 1, 1);
for (int i = 1; i <= f.size(); i++)
c_pow[i] = c_pow[i - 1] * a[0];
std::vector<T> f2(n + 1, 0);
for (int i = 0; i < f.size(); i++) {
T now = f[i];
for (int k = 0; k <= std::min(i, n); k++) {
f2[k] += now * c_pow[i - k];
now *= i - k;
}
}
a[0] = 0;
return SPS_exp_comp(f2, a);
}
template <typename SPS, typename T = typename SPS::value_type> SPS exp(SPS a) {
int N = a.size();
int n = bitwise::log2(N);
SPS ret(N, 1);
for (int d = n - 1; d >= 0; d--) {
const int M = 1 << (n - d - 1);
SPS pre(ret.begin(), ret.begin() + M);
SPS a2(a.begin() + M, a.begin() + 2 * M);
pre *= a2;
std::ranges::copy(pre, ret.begin() + M);
}
return ret;
}
template <typename SPS, typename T = typename SPS::value_type> SPS log(SPS a) {
int n = bitwise::log2(a.size());
assert((1 << n) == a.size() and a[0] == 1);
std::vector<T> lg(n + 1, 0);
for (int k = 1; k <= n; k++)
lg[k] = T(-1) / k;
return SPS_composition(lg, SPS(1, 1) - a);
}
Back to top page
build/pch/stdc++.hpp