This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub suisen-cp/cp-library-cpp
#include "library/polynomial/array_fps_naive.hpp"
#ifndef SUISEN_ARRAY_FPS_NAIVE #define SUISEN_ARRAY_FPS_NAIVE #include <cassert> #include <cmath> #include <limits> #include <type_traits> #include <array> #include "library/type_traits/type_traits.hpp" #include "library/math/modint_extension.hpp" #include "library/math/inv_mods.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(); } #endif // SUISEN_ARRAY_FPS_NAIVE
#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(); }