This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub suisen-cp/cp-library-cpp
#include "library/convolution/array_subset_convolution.hpp"
#ifndef SUISEN_ARRAY_SUBSET_CONVOLUTION #define SUISEN_ARRAY_SUBSET_CONVOLUTION #include "library/polynomial/array_fps_naive.hpp" #include "library/transform/subset.hpp" namespace suisen::array_ranked_subset_transform { template <typename T, std::size_t N> using polynomial_t = ArrayFPSNaive<T, N>; namespace internal { template <typename T, std::size_t N> std::vector<polynomial_t<T, N>> ranked(const std::vector<T>& a) { const int n = a.size(); assert((-n & n) == n); std::vector fs(n, polynomial_t<T, N>{}); for (int i = 0; i < n; ++i) fs[i][__builtin_popcount(i)] = a[i]; return fs; } template <typename T, std::size_t N> std::vector<T> deranked(const std::vector<polynomial_t<T, N>>& polys) { const int n = polys.size(); assert((-n & n) == n); std::vector<T> a(n); for (int i = 0; i < n; ++i) a[i] = polys[i][__builtin_popcount(i)]; return a; } } // namespace suisen::array_ranked_subset_transform::internal template <typename T, std::size_t N> std::vector<polynomial_t<T, N>> ranked_zeta(const std::vector<T>& a) { auto ranked = internal::ranked<T, N>(a); subset_transform::zeta(ranked); return ranked; } template <typename T, std::size_t N> std::vector<T> deranked_mobius(std::vector<polynomial_t<T, N>>& ranked) { subset_transform::mobius(ranked); return internal::deranked<T, N>(ranked); } } // namespace suisen::array_ranked_subset_transform #endif // SUISEN_ARRAY_SUBSET_CONVOLUTION
#line 1 "library/convolution/array_subset_convolution.hpp" #line 1 "library/polynomial/array_fps_naive.hpp" #include <cassert> #include <cmath> #include <limits> #include <type_traits> #include <array> #line 1 "library/type_traits/type_traits.hpp" #line 5 "library/type_traits/type_traits.hpp" #include <iostream> #line 7 "library/type_traits/type_traits.hpp" namespace suisen { template <typename ...Constraints> using constraints_t = std::enable_if_t<std::conjunction_v<Constraints...>, std::nullptr_t>; template <typename T, typename = std::nullptr_t> struct bitnum { static constexpr int value = 0; }; template <typename T> struct bitnum<T, constraints_t<std::is_integral<T>>> { static constexpr int value = std::numeric_limits<std::make_unsigned_t<T>>::digits; }; template <typename T> static constexpr int bitnum_v = bitnum<T>::value; template <typename T, size_t n> struct is_nbit { static constexpr bool value = bitnum_v<T> == n; }; template <typename T, size_t n> static constexpr bool is_nbit_v = is_nbit<T, n>::value; template <typename T, typename = std::nullptr_t> struct safely_multipliable { using type = T; }; template <typename T> struct safely_multipliable<T, constraints_t<std::is_signed<T>, is_nbit<T, 32>>> { using type = long long; }; template <typename T> struct safely_multipliable<T, constraints_t<std::is_signed<T>, is_nbit<T, 64>>> { using type = __int128_t; }; template <typename T> struct safely_multipliable<T, constraints_t<std::is_unsigned<T>, is_nbit<T, 32>>> { using type = unsigned long long; }; template <typename T> struct safely_multipliable<T, constraints_t<std::is_unsigned<T>, is_nbit<T, 64>>> { using type = __uint128_t; }; template <typename T> using safely_multipliable_t = typename safely_multipliable<T>::type; template <typename T, typename = void> struct rec_value_type { using type = T; }; template <typename T> struct rec_value_type<T, std::void_t<typename T::value_type>> { using type = typename rec_value_type<typename T::value_type>::type; }; template <typename T> using rec_value_type_t = typename rec_value_type<T>::type; template <typename T> class is_iterable { template <typename T_> static auto test(T_ e) -> decltype(e.begin(), e.end(), std::true_type{}); static std::false_type test(...); public: static constexpr bool value = decltype(test(std::declval<T>()))::value; }; template <typename T> static constexpr bool is_iterable_v = is_iterable<T>::value; template <typename T> class is_writable { template <typename T_> static auto test(T_ e) -> decltype(std::declval<std::ostream&>() << e, std::true_type{}); static std::false_type test(...); public: static constexpr bool value = decltype(test(std::declval<T>()))::value; }; template <typename T> static constexpr bool is_writable_v = is_writable<T>::value; template <typename T> class is_readable { template <typename T_> static auto test(T_ e) -> decltype(std::declval<std::istream&>() >> e, std::true_type{}); static std::false_type test(...); public: static constexpr bool value = decltype(test(std::declval<T>()))::value; }; template <typename T> static constexpr bool is_readable_v = is_readable<T>::value; } // namespace suisen #line 11 "library/polynomial/array_fps_naive.hpp" #line 1 "library/math/modint_extension.hpp" #line 5 "library/math/modint_extension.hpp" #include <optional> /** * refernce: https://37zigen.com/tonelli-shanks-algorithm/ * calculates x s.t. x^2 = a mod p in O((log p)^2). */ template <typename mint> std::optional<mint> safe_sqrt(mint a) { static int p = mint::mod(); if (a == 0) return std::make_optional(0); if (p == 2) return std::make_optional(a); if (a.pow((p - 1) / 2) != 1) return std::nullopt; mint b = 1; while (b.pow((p - 1) / 2) == 1) ++b; static int tlz = __builtin_ctz(p - 1), q = (p - 1) >> tlz; mint x = a.pow((q + 1) / 2); b = b.pow(q); for (int shift = 2; x * x != a; ++shift) { mint e = a.inv() * x * x; if (e.pow(1 << (tlz - shift)) != 1) x *= b; b *= b; } return std::make_optional(x); } /** * calculates x s.t. x^2 = a mod p in O((log p)^2). * if not exists, raises runtime error. */ template <typename mint> auto sqrt(mint a) -> decltype(mint::mod(), mint()) { return *safe_sqrt(a); } template <typename mint> auto log(mint a) -> decltype(mint::mod(), mint()) { assert(a == 1); return 0; } template <typename mint> auto exp(mint a) -> decltype(mint::mod(), mint()) { assert(a == 0); return 1; } template <typename mint, typename T> auto pow(mint a, T b) -> decltype(mint::mod(), mint()) { return a.pow(b); } template <typename mint> auto inv(mint a) -> decltype(mint::mod(), mint()) { return a.inv(); } #line 1 "library/math/inv_mods.hpp" #include <vector> namespace suisen { template <typename mint> class inv_mods { public: inv_mods() = default; inv_mods(int n) { ensure(n); } const mint& operator[](int i) const { ensure(i); return invs[i]; } static void ensure(int n) { int sz = invs.size(); if (sz < 2) invs = { 0, 1 }, sz = 2; if (sz < n + 1) { invs.resize(n + 1); for (int i = sz; i <= n; ++i) invs[i] = mint(mod - mod / i) * invs[mod % i]; } } private: static std::vector<mint> invs; static constexpr int mod = mint::mod(); }; template <typename mint> std::vector<mint> inv_mods<mint>::invs{}; template <typename mint> std::vector<mint> get_invs(const std::vector<mint>& vs) { const int n = vs.size(); mint p = 1; for (auto& e : vs) { p *= e; assert(e != 0); } mint ip = p.inv(); std::vector<mint> rp(n + 1); rp[n] = 1; for (int i = n - 1; i >= 0; --i) { rp[i] = rp[i + 1] * vs[i]; } std::vector<mint> res(n); for (int i = 0; i < n; ++i) { res[i] = ip * rp[i + 1]; ip *= vs[i]; } return res; } } #line 14 "library/polynomial/array_fps_naive.hpp" namespace suisen { template <typename T, std::size_t N> struct ArrayFPSNaive : std::array<T, N> { static constexpr int SIZE = N; static constexpr int DEG = SIZE - 1; using value_type = T; using element_type = rec_value_type_t<T>; ArrayFPSNaive() { this->fill(value_type{ 0 }); } ArrayFPSNaive(const std::initializer_list<value_type> l) : ArrayFPSNaive() { std::copy(l.begin(), l.end(), this->begin()); } ArrayFPSNaive operator+() const { return ArrayFPSNaive(*this); } ArrayFPSNaive operator-() const { ArrayFPSNaive f(*this); for (auto& e : f) e = -e; return f; } ArrayFPSNaive& operator++() { return ++(*this)[0], * this; } ArrayFPSNaive& operator--() { return --(*this)[0], * this; } ArrayFPSNaive& operator+=(const value_type x) { return (*this)[0] += x, *this; } ArrayFPSNaive& operator-=(const value_type x) { return (*this)[0] -= x, *this; } ArrayFPSNaive& operator+=(const ArrayFPSNaive& g) { for (int i = 0; i < SIZE; ++i) (*this)[i] += g[i]; return *this; } ArrayFPSNaive& operator-=(const ArrayFPSNaive& g) { for (int i = 0; i < SIZE; ++i) (*this)[i] -= g[i]; return *this; } ArrayFPSNaive& operator*=(const ArrayFPSNaive& g) { return *this = *this * g; } ArrayFPSNaive& operator*=(const value_type x) { for (auto& e : *this) e *= x; return *this; } ArrayFPSNaive& operator/=(const ArrayFPSNaive& g) { return *this = *this / g; } ArrayFPSNaive& operator%=(const ArrayFPSNaive& g) { return *this = *this % g; } ArrayFPSNaive& operator<<=(int shamt) { shamt = std::min(shamt, SIZE); for (int i = SIZE - 1; i >= shamt; --i) std::swap((*this)[i], (*this)[i - shamt]); std::fill(this->begin(), this->begin() + shamt, value_type{ 0 }); return *this; } ArrayFPSNaive& operator>>=(int shamt) { shamt = std::min(shamt, SIZE); for (int i = 0; i < SIZE - shamt; ++i) std::swap((*this)[i], (*this)[i + shamt]); std::fill(this->begin() + (SIZE - shamt), this->end(), value_type{ 0 }); return *this; } friend ArrayFPSNaive operator+(ArrayFPSNaive f, const ArrayFPSNaive& g) { f += g; return f; } friend ArrayFPSNaive operator+(ArrayFPSNaive f, const value_type& x) { f += x; return f; } friend ArrayFPSNaive operator-(ArrayFPSNaive f, const ArrayFPSNaive& g) { f -= g; return f; } friend ArrayFPSNaive operator-(ArrayFPSNaive f, const value_type& x) { f -= x; return f; } friend ArrayFPSNaive operator*(const ArrayFPSNaive& f, const ArrayFPSNaive& g) { ArrayFPSNaive h; for (int i = 0; i < SIZE; ++i) for (int j = 0; i + j < SIZE; ++j) h[i + j] += f[i] * g[j]; return h; } friend ArrayFPSNaive operator*(ArrayFPSNaive f, const value_type& x) { f *= x; return f; } friend ArrayFPSNaive operator/(ArrayFPSNaive f, ArrayFPSNaive g) { return std::move(div_mod(std::move(f), std::move(g)).first); } friend ArrayFPSNaive operator%(ArrayFPSNaive f, ArrayFPSNaive g) { return std::move(div_mod(std::move(f), std::move(g)).second); } friend ArrayFPSNaive operator*(const value_type x, ArrayFPSNaive f) { f *= x; return f; } friend ArrayFPSNaive operator<<(ArrayFPSNaive f, const int shamt) { f <<= shamt; return f; } friend ArrayFPSNaive operator>>(ArrayFPSNaive f, const int shamt) { f >>= shamt; return f; } friend std::pair<ArrayFPSNaive, ArrayFPSNaive> div_mod(ArrayFPSNaive f, const ArrayFPSNaive& g) { int fd = DEG, gd = DEG; while (fd >= 0 and f[fd] == value_type{ 0 }) --fd; while (gd >= 0 and g[gd] == value_type{ 0 }) --gd; assert(gd >= 0); if (fd < gd) return { ArrayFPSNaive{}, f }; if (gd == 0) return { f *= g[0].inv(), ArrayFPSNaive{} }; const int k = fd - gd; value_type head_inv = g[gd].inv(); ArrayFPSNaive q; for (int i = k; i >= 0; --i) { value_type div = f[i + gd] * head_inv; q[i] = div; for (int j = 0; j <= gd; ++j) f[i + j] -= div * g[j]; } std::fill(f.begin() + gd, f.end(), value_type{ 0 }); return { std::move(q), std::move(f) }; } ArrayFPSNaive mul(const ArrayFPSNaive& g) const { return (*this) * g; } ArrayFPSNaive diff() const { ArrayFPSNaive g; for (int i = 1; i <= DEG; ++i) g[i - 1] = (*this)[i] * i; g[DEG] = 0; return g; } ArrayFPSNaive intg() const { ArrayFPSNaive g; for (int i = 0; i < DEG; ++i) g[i + 1] = (*this)[i] * invs[i + 1]; return g; } ArrayFPSNaive inv() const { ArrayFPSNaive g; const value_type inv_f0 = ::inv((*this)[0]); g[0] = inv_f0; for (int i = 1; i <= DEG; ++i) { for (int j = 1; j <= i; ++j) g[i] -= g[i - j] * (*this)[j]; g[i] *= inv_f0; } return g; } ArrayFPSNaive exp() const { assert((*this)[0] == value_type{ 0 }); ArrayFPSNaive g; g[0] = value_type{ 1 }; for (int i = 1; i <= DEG; ++i) { for (int j = 1; j <= i; ++j) g[i] += j * g[i - j] * (*this)[j]; g[i] *= invs[i]; } return g; } ArrayFPSNaive log() const { assert((*this)[0] == value_type{ 1 }); ArrayFPSNaive g; g[0] = value_type{ 0 }; for (int i = 1; i <= DEG; ++i) { g[i] = i * (*this)[i]; for (int j = 1; j < i; ++j) g[i] -= (i - j) * g[i - j] * (*this)[j]; g[i] *= invs[i]; } return g; } ArrayFPSNaive pow(const long long k) const { if (k == 0) { ArrayFPSNaive g; g[0] = 1; return g; } int z = 0; while (z < SIZE and (*this)[z] == value_type{ 0 }) ++z; if (z >= DEG / k + 1) return ArrayFPSNaive{}; const int d = DEG - z * k; const int bf = z, bg = z * k; ArrayFPSNaive g; const value_type inv_f0 = ::inv((*this)[bf]); g[bg] = (*this)[bf].pow(k); for (int i = 1; i <= d; ++i) { for (int j = 1; j <= i; ++j) g[bg + i] += (element_type{ k } * j - (i - j)) * g[bg + i - j] * (*this)[bf + j]; g[bg + i] *= inv_f0 * invs[i]; } return g; } ArrayFPSNaive sqrt() const { int dl = 0; while (dl < SIZE and (*this)[dl] == value_type{ 0 }) ++dl; if (dl == SIZE) return ArrayFPSNaive{}; if (dl & 1) assert(false); const int d = DEG - dl / 2; const int bf = dl, bg = bf / 2; ArrayFPSNaive g; g[bg] = ::sqrt((*this)[bf]); value_type inv_2g0 = ::inv(2 * g[bg]); for (int i = 1; i <= d; ++i) { g[bg + i] = (*this)[bf + i]; for (int j = 1; j < i; ++j) g[bg + i] -= g[bg + j] * g[bg + i - j]; g[bg + i] *= inv_2g0; } return g; } value_type eval(value_type x) const { value_type y = 0; for (int i = DEG; i >= 0; --i) y = y * x + (*this)[i]; return y; } private: static inline inv_mods<element_type> invs; }; } // namespace suisen template <typename mint, std::size_t N> auto sqrt(suisen::ArrayFPSNaive<mint, N> a) -> decltype(mint::mod(), suisen::ArrayFPSNaive<mint, N>{}) { return a.sqrt(); } template <typename mint, std::size_t N> auto log(suisen::ArrayFPSNaive<mint, N> a) -> decltype(mint::mod(), suisen::ArrayFPSNaive<mint, N>{}) { return a.log(); } template <typename mint, std::size_t N> auto exp(suisen::ArrayFPSNaive<mint, N> a) -> decltype(mint::mod(), suisen::ArrayFPSNaive<mint, N>{}) { return a.exp(); } template <typename mint, std::size_t N, typename T> auto pow(suisen::ArrayFPSNaive<mint, N> a, const T& b) -> decltype(mint::mod(), suisen::ArrayFPSNaive<mint, N>{}) { return a.pow(b); } template <typename mint, std::size_t N> auto inv(suisen::ArrayFPSNaive<mint, N> a) -> decltype(mint::mod(), suisen::ArrayFPSNaive<mint, N>{}) { return a.inv(); } #line 1 "library/transform/subset.hpp" #line 1 "library/transform/kronecker_power.hpp" #line 5 "library/transform/kronecker_power.hpp" #include <utility> #line 7 "library/transform/kronecker_power.hpp" #line 1 "library/util/default_operator.hpp" namespace suisen { namespace default_operator { template <typename T> auto zero() -> decltype(T { 0 }) { return T { 0 }; } template <typename T> auto one() -> decltype(T { 1 }) { return T { 1 }; } template <typename T> auto add(const T &x, const T &y) -> decltype(x + y) { return x + y; } template <typename T> auto sub(const T &x, const T &y) -> decltype(x - y) { return x - y; } template <typename T> auto mul(const T &x, const T &y) -> decltype(x * y) { return x * y; } template <typename T> auto div(const T &x, const T &y) -> decltype(x / y) { return x / y; } template <typename T> auto mod(const T &x, const T &y) -> decltype(x % y) { return x % y; } template <typename T> auto neg(const T &x) -> decltype(-x) { return -x; } template <typename T> auto inv(const T &x) -> decltype(one<T>() / x) { return one<T>() / x; } } // default_operator namespace default_operator_noref { template <typename T> auto zero() -> decltype(T { 0 }) { return T { 0 }; } template <typename T> auto one() -> decltype(T { 1 }) { return T { 1 }; } template <typename T> auto add(T x, T y) -> decltype(x + y) { return x + y; } template <typename T> auto sub(T x, T y) -> decltype(x - y) { return x - y; } template <typename T> auto mul(T x, T y) -> decltype(x * y) { return x * y; } template <typename T> auto div(T x, T y) -> decltype(x / y) { return x / y; } template <typename T> auto mod(T x, T y) -> decltype(x % y) { return x % y; } template <typename T> auto neg(T x) -> decltype(-x) { return -x; } template <typename T> auto inv(T x) -> decltype(one<T>() / x) { return one<T>() / x; } } // default_operator } // namespace suisen #line 9 "library/transform/kronecker_power.hpp" namespace suisen { namespace kronecker_power_transform { namespace internal { template <typename UnitTransform, typename ReferenceGetter, std::size_t... Seq> void unit_transform(UnitTransform transform, ReferenceGetter ref_getter, std::index_sequence<Seq...>) { transform(ref_getter(Seq)...); } } template <typename T, std::size_t D, auto unit_transform> void kronecker_power_transform(std::vector<T> &x) { const std::size_t n = x.size(); for (std::size_t block = 1; block < n; block *= D) { for (std::size_t l = 0; l < n; l += D * block) { for (std::size_t offset = l; offset < l + block; ++offset) { const auto ref_getter = [&](std::size_t i) -> T& { return x[offset + i * block]; }; internal::unit_transform(unit_transform, ref_getter, std::make_index_sequence<D>()); } } } } template <typename T, typename UnitTransform> void kronecker_power_transform(std::vector<T> &x, const std::size_t D, UnitTransform unit_transform) { const std::size_t n = x.size(); std::vector<T> work(D); for (std::size_t block = 1; block < n; block *= D) { for (std::size_t l = 0; l < n; l += D * block) { for (std::size_t offset = l; offset < l + block; ++offset) { for (std::size_t i = 0; i < D; ++i) work[i] = x[offset + i * block]; unit_transform(work); for (std::size_t i = 0; i < D; ++i) x[offset + i * block] = work[i]; } } } } template <typename T, auto e = default_operator::zero<T>, auto add = default_operator::add<T>, auto mul = default_operator::mul<T>> auto kronecker_power_transform(std::vector<T> &x, const std::vector<std::vector<T>> &A) -> decltype(e(), add(std::declval<T>(), std::declval<T>()), mul(std::declval<T>(), std::declval<T>()), void()) { const std::size_t D = A.size(); assert(D == A[0].size()); auto unit_transform = [&](std::vector<T> &x) { std::vector<T> y(D, e()); for (std::size_t i = 0; i < D; ++i) for (std::size_t j = 0; j < D; ++j) { y[i] = add(y[i], mul(A[i][j], x[j])); } x.swap(y); }; kronecker_power_transform<T>(x, D, unit_transform); } } } // namespace suisen #line 5 "library/transform/subset.hpp" namespace suisen::subset_transform { namespace internal { template <typename T, auto add = default_operator::add<T>> void zeta_unit_transform(T &x0, T &x1) { // 1, 0 x1 = add(x1, x0); // 1, 1 } template <typename T, auto sub = default_operator::sub<T>> void mobius_unit_transform(T &x0, T &x1) { // 1, 0 x1 = sub(x1, x0); // -1, 1 } } // namespace internal using kronecker_power_transform::kronecker_power_transform; template <typename T, auto add = default_operator::add<T>> void zeta(std::vector<T> &a) { kronecker_power_transform<T, 2, internal::zeta_unit_transform<T, add>>(a); } template <typename T, auto sub = default_operator::sub<T>> void mobius(std::vector<T> &a) { kronecker_power_transform<T, 2, internal::mobius_unit_transform<T, sub>>(a); } } // namespace suisen::subset_transform #line 6 "library/convolution/array_subset_convolution.hpp" namespace suisen::array_ranked_subset_transform { template <typename T, std::size_t N> using polynomial_t = ArrayFPSNaive<T, N>; namespace internal { template <typename T, std::size_t N> std::vector<polynomial_t<T, N>> ranked(const std::vector<T>& a) { const int n = a.size(); assert((-n & n) == n); std::vector fs(n, polynomial_t<T, N>{}); for (int i = 0; i < n; ++i) fs[i][__builtin_popcount(i)] = a[i]; return fs; } template <typename T, std::size_t N> std::vector<T> deranked(const std::vector<polynomial_t<T, N>>& polys) { const int n = polys.size(); assert((-n & n) == n); std::vector<T> a(n); for (int i = 0; i < n; ++i) a[i] = polys[i][__builtin_popcount(i)]; return a; } } // namespace suisen::array_ranked_subset_transform::internal template <typename T, std::size_t N> std::vector<polynomial_t<T, N>> ranked_zeta(const std::vector<T>& a) { auto ranked = internal::ranked<T, N>(a); subset_transform::zeta(ranked); return ranked; } template <typename T, std::size_t N> std::vector<T> deranked_mobius(std::vector<polynomial_t<T, N>>& ranked) { subset_transform::mobius(ranked); return internal::deranked<T, N>(ranked); } } // namespace suisen::array_ranked_subset_transform