#line 1 "other/a.cc" // #undef _GLIBCXX_DEBUG // #define NDEBUG #include #line 2 "Library/lib/alias" /** * @file alias * @brief Alias */ #line 13 "Library/lib/alias" #line 1 "Library/lib/bit" #if __cplusplus > 201703L #include #else #ifndef _GLIBCXX_BIT #define _GLIBCXX_BIT 1 #include #include namespace std { template constexpr _Tp __rotl(_Tp __x, int __s) noexcept { constexpr auto _Nd = numeric_limits<_Tp>::digits; const int __r = __s % _Nd; if (__r == 0) return __x; else if (__r > 0) return (__x << __r) | (__x >> ((_Nd - __r) % _Nd)); else return (__x >> -__r) | (__x << ((_Nd + __r) % _Nd)); // rotr(x, -r) } template constexpr _Tp __rotr(_Tp __x, int __s) noexcept { constexpr auto _Nd = numeric_limits<_Tp>::digits; const int __r = __s % _Nd; if (__r == 0) return __x; else if (__r > 0) return (__x >> __r) | (__x << ((_Nd - __r) % _Nd)); else return (__x << -__r) | (__x >> ((_Nd + __r) % _Nd)); // rotl(x, -r) } template constexpr int __countl_zero(_Tp __x) noexcept { constexpr auto _Nd = numeric_limits<_Tp>::digits; if (__x == 0) return _Nd; constexpr auto _Nd_ull = numeric_limits::digits; constexpr auto _Nd_ul = numeric_limits::digits; constexpr auto _Nd_u = numeric_limits::digits; if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_u) { constexpr int __diff = _Nd_u - _Nd; return __builtin_clz(__x) - __diff; } else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ul) { constexpr int __diff = _Nd_ul - _Nd; return __builtin_clzl(__x) - __diff; } else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ull) { constexpr int __diff = _Nd_ull - _Nd; return __builtin_clzll(__x) - __diff; } else // (_Nd > _Nd_ull) { static_assert(_Nd <= (2 * _Nd_ull), "Maximum supported integer size is 128-bit"); unsigned long long __high = __x >> _Nd_ull; if (__high != 0) { constexpr int __diff = (2 * _Nd_ull) - _Nd; return __builtin_clzll(__high) - __diff; } constexpr auto __max_ull = numeric_limits::max(); unsigned long long __low = __x & __max_ull; return (_Nd - _Nd_ull) + __builtin_clzll(__low); } } template constexpr int __countl_one(_Tp __x) noexcept { if (__x == numeric_limits<_Tp>::max()) return numeric_limits<_Tp>::digits; return __countl_zero<_Tp>((_Tp)~__x); } template constexpr int __countr_zero(_Tp __x) noexcept { constexpr auto _Nd = numeric_limits<_Tp>::digits; if (__x == 0) return _Nd; constexpr auto _Nd_ull = numeric_limits::digits; constexpr auto _Nd_ul = numeric_limits::digits; constexpr auto _Nd_u = numeric_limits::digits; if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_u) return __builtin_ctz(__x); else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ul) return __builtin_ctzl(__x); else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ull) return __builtin_ctzll(__x); else // (_Nd > _Nd_ull) { static_assert(_Nd <= (2 * _Nd_ull), "Maximum supported integer size is 128-bit"); constexpr auto __max_ull = numeric_limits::max(); unsigned long long __low = __x & __max_ull; if (__low != 0) return __builtin_ctzll(__low); unsigned long long __high = __x >> _Nd_ull; return __builtin_ctzll(__high) + _Nd_ull; } } template constexpr int __countr_one(_Tp __x) noexcept { if (__x == numeric_limits<_Tp>::max()) return numeric_limits<_Tp>::digits; return __countr_zero((_Tp)~__x); } template constexpr int __popcount(_Tp __x) noexcept { constexpr auto _Nd = numeric_limits<_Tp>::digits; if (__x == 0) return 0; constexpr auto _Nd_ull = numeric_limits::digits; constexpr auto _Nd_ul = numeric_limits::digits; constexpr auto _Nd_u = numeric_limits::digits; if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_u) return __builtin_popcount(__x); else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ul) return __builtin_popcountl(__x); else if _GLIBCXX17_CONSTEXPR (_Nd <= _Nd_ull) return __builtin_popcountll(__x); else // (_Nd > _Nd_ull) { static_assert(_Nd <= (2 * _Nd_ull), "Maximum supported integer size is 128-bit"); constexpr auto __max_ull = numeric_limits::max(); unsigned long long __low = __x & __max_ull; unsigned long long __high = __x >> _Nd_ull; return __builtin_popcountll(__low) + __builtin_popcountll(__high); } } template constexpr bool __has_single_bit(_Tp __x) noexcept { return __popcount(__x) == 1; } template constexpr _Tp __bit_ceil(_Tp __x) noexcept { constexpr auto _Nd = numeric_limits<_Tp>::digits; if (__x == 0 || __x == 1) return 1; auto __shift_exponent = _Nd - __countl_zero((_Tp)(__x - 1u)); #ifdef _GLIBCXX_HAVE_BUILTIN_IS_CONSTANT_EVALUATED if (!__builtin_is_constant_evaluated()) { __glibcxx_assert(__shift_exponent != numeric_limits<_Tp>::digits); } #endif using __promoted_type = decltype(__x << 1); if _GLIBCXX17_CONSTEXPR (!is_same<__promoted_type, _Tp>::value) { const int __extra_exp = sizeof(__promoted_type) / sizeof(_Tp) / 2; __shift_exponent |= (__shift_exponent & _Nd) << __extra_exp; } return (_Tp)1u << __shift_exponent; } template constexpr _Tp __bit_floor(_Tp __x) noexcept { constexpr auto _Nd = numeric_limits<_Tp>::digits; if (__x == 0) return 0; return (_Tp)1u << (_Nd - __countl_zero((_Tp)(__x >> 1))); } template constexpr _Tp __bit_width(_Tp __x) noexcept { constexpr auto _Nd = numeric_limits<_Tp>::digits; return _Nd - __countl_zero(__x); } } // namespace std #endif #endif #line 15 "Library/lib/alias" namespace workspace { constexpr char eol = '\n'; using namespace std; using i32 = int_least32_t; using u32 = uint_least32_t; using i64 = int_least64_t; using u64 = uint_least64_t; #ifdef __SIZEOF_INT128__ using i128 = __int128_t; using u128 = __uint128_t; #else #warning 128bit integer is not available. #endif template > using priority_queue = std::priority_queue, Comp>; template using stack = std::stack>; template constexpr _Tp __bsf(_Tp __x) noexcept { return std::__countr_zero(__x); } template constexpr _Tp __bsr(_Tp __x) noexcept { return std::__bit_width(__x) - 1; } } // namespace workspace #line 6 "other/a.cc" // #include "lib/cxx20" // #include "lib/direct" // #include "lib/opt" #line 2 "Library/src/sys/clock.hpp" /* * @fn clock.hpp * @brief Clock */ #line 9 "Library/src/sys/clock.hpp" namespace workspace { using namespace std::chrono; namespace internal { // The start time of the program. const auto start_time{system_clock::now()}; } // namespace internal /* * @fn elapsed * @return elapsed time of the program */ int64_t elapsed() { const auto end_time{system_clock::now()}; return duration_cast(end_time - internal::start_time).count(); } } // namespace workspace #line 2 "Library/src/sys/ejection.hpp" /** * @file ejection.hpp * @brief Ejection */ #line 9 "Library/src/sys/ejection.hpp" namespace workspace { namespace internal { struct ejection { bool exit = 0; }; } // namespace internal /** * @brief eject from a try block, throw nullptr * @param arg output */ template void eject(Tp const &arg) { std::cout << arg << "\n"; throw internal::ejection{}; } void exit() { throw internal::ejection{true}; } } // namespace workspace #line 2 "Library/src/sys/iteration.hpp" /** * @file iteration.hpp * @brief Case Iteration */ #line 9 "Library/src/sys/iteration.hpp" #line 11 "Library/src/sys/iteration.hpp" namespace workspace { void main(); struct { // 1-indexed unsigned current{0}; unsigned total{1}; void read() { (std::cin >> total).ignore(); } int iterate() { static bool once = false; assert(!once); once = true; while (current++ < total) { try { main(); } catch (internal::ejection const &status) { if (status.exit) break; } } return 0; } } case_info; } // namespace workspace #line 2 "Library/src/utils/cat.hpp" /** * @file cat.hpp * @brief Cat */ #line 9 "Library/src/utils/cat.hpp" namespace workspace { template constexpr C1 &&cat(C1 &&__c1, C2 const &__c2) noexcept { __c1.insert(__c1.end(), std::begin(__c2), std::end(__c2)); return __c1; } } // namespace workspace #line 2 "Library/src/utils/chval.hpp" /* * @file chval.hpp * @brief Change Less/Greater */ #line 9 "Library/src/utils/chval.hpp" namespace workspace { /* * @fn chle * @brief Substitute y for x if comp(y, x) is true. * @param x Reference * @param y Const reference * @param comp Compare function * @return Whether or not x is updated */ template > bool chle(Tp &x, const Tp &y, Comp comp = Comp()) { return comp(y, x) ? x = y, true : false; } /* * @fn chge * @brief Substitute y for x if comp(x, y) is true. * @param x Reference * @param y Const reference * @param comp Compare function * @return Whether or not x is updated */ template > bool chge(Tp &x, const Tp &y, Comp comp = Comp()) { return comp(x, y) ? x = y, true : false; } } // namespace workspace #line 4 "Library/lib/utils" // #include "src/utils/fixed_point.hpp" // #include "src/utils/grid.hpp" // #include "src/utils/hash.hpp" #line 2 "Library/src/utils/io/istream.hpp" /** * @file istream.hpp * @brief Input Stream */ #include #line 13 "Library/src/utils/io/istream.hpp" #line 2 "Library/src/utils/sfinae.hpp" /** * @file sfinae.hpp * @brief SFINAE */ #line 10 "Library/src/utils/sfinae.hpp" #include #ifdef __SIZEOF_INT128__ #define __INT128_DEFINED__ 1 #else #define __INT128_DEFINED__ 0 #endif namespace std { #if __INT128_DEFINED__ template <> struct make_signed<__uint128_t> { using type = __int128_t; }; template <> struct make_signed<__int128_t> { using type = __int128_t; }; template <> struct make_unsigned<__uint128_t> { using type = __uint128_t; }; template <> struct make_unsigned<__int128_t> { using type = __uint128_t; }; #endif } // namespace std namespace workspace { template class trait> using enable_if_trait_type = typename std::enable_if::value>::type; template using element_type = typename std::decay()))>::type; template struct has_begin : std::false_type {}; template struct has_begin()), nullptr)> : std::true_type {}; template struct mapped_of { using type = element_type; }; template struct mapped_of::first_type> { using type = typename T::mapped_type; }; template using mapped_type = typename mapped_of::type; template struct is_integral_ext : std::false_type {}; template struct is_integral_ext< T, typename std::enable_if::value>::type> : std::true_type {}; #if __INT128_DEFINED__ template <> struct is_integral_ext<__int128_t> : std::true_type {}; template <> struct is_integral_ext<__uint128_t> : std::true_type {}; #endif #if __cplusplus >= 201402 template constexpr static bool is_integral_ext_v = is_integral_ext::value; #endif template struct multiplicable_uint { using type = uint_least32_t; }; template struct multiplicable_uint< T, typename std::enable_if<(2 < sizeof(T)) && (!__INT128_DEFINED__ || sizeof(T) <= 4)>::type> { using type = uint_least64_t; }; #if __INT128_DEFINED__ template struct multiplicable_uint::type> { using type = __uint128_t; }; #endif template struct multiplicable_int { using type = typename std::make_signed::type>::type; }; } // namespace workspace #line 15 "Library/src/utils/io/istream.hpp" namespace workspace { namespace internal { template struct istream_helper { istream_helper(std::istream &is, Tp &x) { if constexpr (has_begin::value) for (auto &&e : x) istream_helper::type>(is, e); else static_assert(has_begin::value, "istream unsupported type."); } }; template struct istream_helper< Tp, decltype(std::declval() >> std::declval())>>(), nullptr)> { istream_helper(std::istream &is, Tp &x) { is >> x; } }; #ifdef __SIZEOF_INT128__ template <> struct istream_helper<__int128_t, std::nullptr_t> { istream_helper(std::istream &is, __int128_t &x) { std::string s; is >> s; bool negative = s.front() == '-' ? s.erase(s.begin()), true : false; x = 0; for (char e : s) x = x * 10 + e - '0'; if (negative) x = -x; } }; template <> struct istream_helper<__uint128_t, std::nullptr_t> { istream_helper(std::istream &is, __uint128_t &x) { std::string s; is >> s; bool negative = s.front() == '-' ? s.erase(s.begin()), true : false; x = 0; for (char e : s) x = x * 10 + e - '0'; if (negative) x = -x; } }; #endif // INT128 template struct istream_helper> { istream_helper(std::istream &is, std::pair &x) { istream_helper(is, x.first), istream_helper(is, x.second); } }; template struct istream_helper> { istream_helper(std::istream &is, std::tuple &x) { iterate(is, x); } private: template void iterate(std::istream &is, Tp &x) { if constexpr (N == std::tuple_size::value) return; else istream_helper::type>(is, std::get(x)), iterate(is, x); } }; } // namespace internal /** * @brief A wrapper class for std::istream. */ class istream : public std::istream { public: /** * @brief Wrapped operator. */ template istream &operator>>(Tp &x) { internal::istream_helper(*this, x); if (std::istream::fail()) { static auto once = atexit([] { std::cerr << "\n\033[43m\033[30mwarning: failed to read \'" << abi::__cxa_demangle(typeid(Tp).name(), 0, 0, 0) << "\'.\033[0m\n\n"; }); assert(!once); } return *this; } }; namespace internal { auto *const cin_ptr = (istream *)&std::cin; } auto &cin = *internal::cin_ptr; } // namespace workspace #line 2 "Library/src/utils/io/ostream.hpp" /* * @file ostream.hpp * @brief Output Stream */ #line 10 "Library/src/utils/io/ostream.hpp" namespace workspace { template std::ostream &operator<<(std::ostream &os, const std::pair &p) { return os << p.first << ' ' << p.second; } template struct tuple_os { static std::ostream &apply(std::ostream &os, const tuple_t &t) { tuple_os::apply(os, t); return os << ' ' << std::get(t); } }; template struct tuple_os { static std::ostream &apply(std::ostream &os, const tuple_t &t) { return os << std::get<0>(t); } }; template struct tuple_os { static std::ostream &apply(std::ostream &os, const tuple_t &t) { return os; } }; template std::ostream &operator<<(std::ostream &os, const std::tuple &t) { return tuple_os, std::tuple_size>::value - 1>::apply(os, t); } template ()))> typename std::enable_if< !std::is_same::type, std::string>::value && !std::is_same::type, char *>::value, std::ostream &>::type operator<<(std::ostream &os, const Container &cont) { bool head = true; for (auto &&e : cont) head ? head = 0 : (os << ' ', 0), os << e; return os; } } // namespace workspace #line 9 "Library/lib/utils" // #include "src/utils/io/read.hpp" #line 2 "Library/src/utils/io/setup.hpp" /* * @file setup.hpp * @brief I/O Setup */ #line 10 "Library/src/utils/io/setup.hpp" namespace workspace { /* * @fn io_setup * @brief Setup I/O. * @param precision Standard output precision */ void io_setup(int precision) { std::ios::sync_with_stdio(false); std::cin.tie(nullptr); std::cout << std::fixed << std::setprecision(precision); #ifdef _buffer_check atexit([] { char bufc; if (std::cin >> bufc) std::cerr << "\n\033[43m\033[30mwarning: buffer not empty.\033[0m\n\n"; }); #endif } } // namespace workspace #line 2 "Library/src/utils/iterator/category.hpp" /* * @file category.hpp * @brief Iterator Category */ #line 10 "Library/src/utils/iterator/category.hpp" namespace workspace { /* * @tparam Tuple Tuple of iterator types */ template ::value - 1> struct common_iterator_category { using type = typename std::common_type< typename common_iterator_category::type, typename std::iterator_traits::type>::iterator_category>::type; }; template struct common_iterator_category { using type = typename std::iterator_traits< typename std::tuple_element<0, Tuple>::type>::iterator_category; }; } // namespace workspace #line 12 "Library/lib/utils" // #include "src/utils/iterator/reverse.hpp" #line 2 "Library/src/utils/make_vector.hpp" /* * @file make_vector.hpp * @brief Multi-dimensional Vector */ #if __cplusplus >= 201703L #include #include namespace workspace { /* * @brief Make a multi-dimensional vector. * @tparam Tp type of the elements * @tparam N dimension * @tparam S integer type * @param sizes The size of each dimension * @param init The initial value */ template constexpr auto make_vector([[maybe_unused]] S *sizes, Tp const &init = Tp()) { static_assert(std::is_convertible_v); if constexpr (N) return std::vector(*sizes, make_vector(std::next(sizes), init)); else return init; } /* * @brief Make a multi-dimensional vector. * @param sizes The size of each dimension * @param init The initial value */ template constexpr auto make_vector(const S (&sizes)[N], Tp const &init = Tp()) { return make_vector((S *)sizes, init); } /* * @brief Make a multi-dimensional vector. * @param sizes The size of each dimension * @param init The initial value */ template constexpr auto make_vector([[maybe_unused]] std::array const &sizes, Tp const &init = Tp()) { static_assert(std::is_convertible_v); if constexpr (I == N) return init; else return std::vector(sizes[I], make_vector(sizes, init)); } /* * @brief Make a multi-dimensional vector. * @param sizes The size of each dimension * @param init The initial value */ template constexpr auto make_vector([[maybe_unused]] std::tuple const &sizes, Tp const &init = Tp()) { using tuple_type = std::tuple; if constexpr (I == std::tuple_size_v || I == N) return init; else { static_assert( std::is_convertible_v, size_t>); return std::vector(std::get(sizes), make_vector(sizes, init)); } } /* * @brief Make a multi-dimensional vector. * @param sizes The size of each dimension * @param init The initial value */ template constexpr auto make_vector(std::pair const &sizes, Tp const &init = Tp()) { static_assert(std::is_convertible_v); static_assert(std::is_convertible_v); return make_vector({(size_t)sizes.first, (size_t)sizes.second}, init); } } // namespace workspace #endif #line 14 "Library/lib/utils" // #include "src/utils/py-like/enumerate.hpp" #line 2 "Library/src/utils/py-like/range.hpp" /** * @file range.hpp * @brief Range */ #line 9 "Library/src/utils/py-like/range.hpp" #line 2 "Library/src/utils/iterator/reverse.hpp" /* * @file reverse_iterator.hpp * @brief Reverse Iterator */ #if __cplusplus >= 201703L #include #include namespace workspace { /* * @class reverse_iterator * @brief Wrapper class for `std::reverse_iterator`. * @see http://gcc.gnu.org/PR51823 */ template class reverse_iterator : public std::reverse_iterator { using base_std = std::reverse_iterator; std::optional deref; public: using base_std::reverse_iterator; constexpr typename base_std::reference operator*() noexcept { if (!deref) { Iterator tmp = base_std::current; deref = *--tmp; } return deref.value(); } constexpr reverse_iterator &operator++() noexcept { base_std::operator++(); deref.reset(); return *this; } constexpr reverse_iterator &operator--() noexcept { base_std::operator++(); deref.reset(); return *this; } constexpr reverse_iterator operator++(int) noexcept { base_std::operator++(); deref.reset(); return *this; } constexpr reverse_iterator operator--(int) noexcept { base_std::operator++(); deref.reset(); return *this; } }; } // namespace workspace #endif #line 2 "Library/src/utils/py-like/reversed.hpp" /** * @file reversed.hpp * @brief Reversed */ #include #line 10 "Library/src/utils/py-like/reversed.hpp" namespace workspace { namespace internal { template class reversed { Container cont; public: constexpr reversed(Container &&cont) : cont(cont) {} constexpr auto begin() { return std::rbegin(cont); } constexpr auto end() { return std::rend(cont); } }; } // namespace internal template constexpr auto reversed(Container &&cont) noexcept { return internal::reversed{std::forward(cont)}; } template constexpr auto reversed(std::initializer_list &&cont) noexcept { return internal::reversed>{ std::forward>(cont)}; } } // namespace workspace #line 12 "Library/src/utils/py-like/range.hpp" #if __cplusplus >= 201703L namespace workspace { template class range { Index first, last; public: class iterator { Index current; public: using difference_type = std::ptrdiff_t; using value_type = Index; using reference = typename std::add_const::type &; using pointer = iterator; using iterator_category = std::bidirectional_iterator_tag; constexpr iterator(Index const &__i = Index()) noexcept : current(__i) {} constexpr bool operator==(iterator const &rhs) const noexcept { return current == rhs.current; } constexpr bool operator!=(iterator const &rhs) const noexcept { return current != rhs.current; } constexpr iterator &operator++() noexcept { ++current; return *this; } constexpr iterator &operator--() noexcept { --current; return *this; } constexpr reference operator*() const noexcept { return current; } }; constexpr range(Index first, Index last) noexcept : first(first), last(last) {} constexpr range(Index last) noexcept : first(), last(last) {} constexpr iterator begin() const noexcept { return iterator{first}; } constexpr iterator end() const noexcept { return iterator{last}; } constexpr reverse_iterator rbegin() const noexcept { return reverse_iterator(end()); } constexpr reverse_iterator rend() const noexcept { return reverse_iterator(begin()); } }; template constexpr auto rrange(Args &&...args) noexcept { return internal::reversed(range(std::forward(args)...)); } } // namespace workspace #endif #line 16 "Library/lib/utils" // #include "src/utils/py-like/reversed.hpp" // #include "src/utils/py-like/zip.hpp" // #include "src/utils/rand/rng.hpp" // #include "src/utils/rand/shuffle.hpp" // #include "src/utils/round_div.hpp" // #include "src/utils/sfinae.hpp" #line 11 "other/a.cc" signed main() { using namespace workspace; io_setup(15); /* given case_info.read(); //*/ /* unspecified case_info.total = -1; //*/ return case_info.iterate(); } #include using namespace std; template struct FormalPowerSeries : vector { using vector::vector; using FPS = FormalPowerSeries; FPS &operator+=(const FPS &r) { if (r.size() > this->size()) this->resize(r.size()); for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i]; return *this; } FPS &operator+=(const mint &r) { if (this->empty()) this->resize(1); (*this)[0] += r; return *this; } FPS &operator-=(const FPS &r) { if (r.size() > this->size()) this->resize(r.size()); for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i]; return *this; } FPS &operator-=(const mint &r) { if (this->empty()) this->resize(1); (*this)[0] -= r; return *this; } FPS &operator*=(const mint &v) { for (int k = 0; k < (int)this->size(); k++) (*this)[k] *= v; return *this; } FPS &operator/=(const FPS &r) { if (this->size() < r.size()) { this->clear(); return *this; } int n = this->size() - r.size() + 1; if ((int)r.size() <= 64) { FPS f(*this), g(r); g.shrink(); mint coeff = g.back().inverse(); for (auto &x : g) x *= coeff; int deg = (int)f.size() - (int)g.size() + 1; int gs = g.size(); FPS quo(deg); for (int i = deg - 1; i >= 0; i--) { quo[i] = f[i + gs - 1]; for (int j = 0; j < gs; j++) f[i + j] -= quo[i] * g[j]; } *this = quo * coeff; this->resize(n, mint(0)); return *this; } return *this = ((*this).rev().pre(n) * r.rev().inv(n)).pre(n).rev(); } FPS &operator%=(const FPS &r) { *this -= *this / r * r; shrink(); return *this; } FPS operator+(const FPS &r) const { return FPS(*this) += r; } FPS operator+(const mint &v) const { return FPS(*this) += v; } FPS operator-(const FPS &r) const { return FPS(*this) -= r; } FPS operator-(const mint &v) const { return FPS(*this) -= v; } FPS operator*(const FPS &r) const { return FPS(*this) *= r; } FPS operator*(const mint &v) const { return FPS(*this) *= v; } FPS operator/(const FPS &r) const { return FPS(*this) /= r; } FPS operator%(const FPS &r) const { return FPS(*this) %= r; } FPS operator-() const { FPS ret(this->size()); for (int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i]; return ret; } void shrink() { while (this->size() && this->back() == mint(0)) this->pop_back(); } FPS rev() const { FPS ret(*this); reverse(begin(ret), end(ret)); return ret; } FPS dot(FPS r) const { FPS ret(min(this->size(), r.size())); for (int i = 0; i < (int)ret.size(); i++) ret[i] = (*this)[i] * r[i]; return ret; } FPS pre(int sz) const { return FPS(begin(*this), begin(*this) + min((int)this->size(), sz)); } FPS operator>>(int sz) const { if ((int)this->size() <= sz) return {}; FPS ret(*this); ret.erase(ret.begin(), ret.begin() + sz); return ret; } FPS operator<<(int sz) const { FPS ret(*this); ret.insert(ret.begin(), sz, mint(0)); return ret; } FPS diff() const { const int n = (int)this->size(); FPS ret(max(0, n - 1)); mint one(1), coeff(1); for (int i = 1; i < n; i++) { ret[i - 1] = (*this)[i] * coeff; coeff += one; } return ret; } FPS integral() const { const int n = (int)this->size(); FPS ret(n + 1); ret[0] = mint(0); if (n > 0) ret[1] = mint(1); auto mod = mint::get_mod(); for (int i = 2; i <= n; i++) ret[i] = (-ret[mod % i]) * (mod / i); for (int i = 0; i < n; i++) ret[i + 1] *= (*this)[i]; return ret; } mint eval(mint x) const { mint r = 0, w = 1; for (auto &v : *this) r += w * v, w *= x; return r; } FPS log(int deg = -1) const { assert((*this)[0] == mint(1)); if (deg == -1) deg = (int)this->size(); return (this->diff() * this->inv(deg)).pre(deg - 1).integral(); } FPS pow(int64_t k, int deg = -1) const { const int n = (int)this->size(); if (deg == -1) deg = n; for (int i = 0; i < n; i++) { if ((*this)[i] != mint(0)) { if (i * k > deg) return FPS(deg, mint(0)); mint rev = mint(1) / (*this)[i]; FPS ret = (((*this * rev) >> i).log() * k).exp() * ((*this)[i].pow(k)); ret = (ret << (i * k)).pre(deg); if ((int)ret.size() < deg) ret.resize(deg, mint(0)); return ret; } } return FPS(deg, mint(0)); } static void *ntt_ptr; static void set_fft(); FPS &operator*=(const FPS &r); void ntt(); void intt(); void ntt_doubling(); static int ntt_pr(); FPS inv(int deg = -1) const; FPS exp(int deg = -1) const; }; template void *FormalPowerSeries::ntt_ptr = nullptr; template struct LazyMontgomeryModInt { using mint = LazyMontgomeryModInt; using i32 = int32_t; using u32 = uint32_t; using u64 = uint64_t; static constexpr u32 get_r() { u32 ret = mod; for (i32 i = 0; i < 4; ++i) ret *= 2 - mod * ret; return ret; } static constexpr u32 r = get_r(); static constexpr u32 n2 = -u64(mod) % mod; static_assert(r * mod == 1, "invalid, r * mod != 1"); static_assert(mod < (1 << 30), "invalid, mod >= 2 ^ 30"); static_assert((mod & 1) == 1, "invalid, mod % 2 == 0"); u32 a; constexpr LazyMontgomeryModInt() : a(0) {} constexpr LazyMontgomeryModInt(const int64_t &b) : a(reduce(u64(b % mod + mod) * n2)){}; static constexpr u32 reduce(const u64 &b) { return (b + u64(u32(b) * u32(-r)) * mod) >> 32; } constexpr mint &operator+=(const mint &b) { if (i32(a += b.a - 2 * mod) < 0) a += 2 * mod; return *this; } constexpr mint &operator-=(const mint &b) { if (i32(a -= b.a) < 0) a += 2 * mod; return *this; } constexpr mint &operator*=(const mint &b) { a = reduce(u64(a) * b.a); return *this; } constexpr mint &operator/=(const mint &b) { *this *= b.inverse(); return *this; } constexpr mint operator+(const mint &b) const { return mint(*this) += b; } constexpr mint operator-(const mint &b) const { return mint(*this) -= b; } constexpr mint operator*(const mint &b) const { return mint(*this) *= b; } constexpr mint operator/(const mint &b) const { return mint(*this) /= b; } constexpr bool operator==(const mint &b) const { return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a); } constexpr bool operator!=(const mint &b) const { return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a); } constexpr mint operator-() const { return mint() - mint(*this); } constexpr mint pow(u64 n) const { mint ret(1), mul(*this); while (n > 0) { if (n & 1) ret *= mul; mul *= mul; n >>= 1; } return ret; } constexpr mint inverse() const { return pow(mod - 2); } friend ostream &operator<<(ostream &os, const mint &b) { return os << b.get(); } friend istream &operator>>(istream &is, mint &b) { int64_t t; is >> t; b = LazyMontgomeryModInt(t); return (is); } constexpr u32 get() const { u32 ret = reduce(a); return ret >= mod ? ret - mod : ret; } static constexpr u32 get_mod() { return mod; } }; // __attribute__((target("sse4.2"))) __attribute__((always_inline)) __m128i my128_mullo_epu32(const __m128i &a, const __m128i &b) { return _mm_mullo_epi32(a, b); } __attribute__((target("sse4.2"))) __attribute__((always_inline)) __m128i my128_mulhi_epu32(const __m128i &a, const __m128i &b) { __m128i a13 = _mm_shuffle_epi32(a, 0xF5); __m128i b13 = _mm_shuffle_epi32(b, 0xF5); __m128i prod02 = _mm_mul_epu32(a, b); __m128i prod13 = _mm_mul_epu32(a13, b13); __m128i prod = _mm_unpackhi_epi64(_mm_unpacklo_epi32(prod02, prod13), _mm_unpackhi_epi32(prod02, prod13)); return prod; } __attribute__((target("sse4.2"))) __attribute__((always_inline)) __m128i montgomery_mul_128(const __m128i &a, const __m128i &b, const __m128i &r, const __m128i &m1) { return _mm_sub_epi32( _mm_add_epi32(my128_mulhi_epu32(a, b), m1), my128_mulhi_epu32(my128_mullo_epu32(my128_mullo_epu32(a, b), r), m1)); } __attribute__((target("sse4.2"))) __attribute__((always_inline)) __m128i montgomery_add_128(const __m128i &a, const __m128i &b, const __m128i &m2, const __m128i &m0) { __m128i ret = _mm_sub_epi32(_mm_add_epi32(a, b), m2); return _mm_add_epi32(_mm_and_si128(_mm_cmpgt_epi32(m0, ret), m2), ret); } __attribute__((target("sse4.2"))) __attribute__((always_inline)) __m128i montgomery_sub_128(const __m128i &a, const __m128i &b, const __m128i &m2, const __m128i &m0) { __m128i ret = _mm_sub_epi32(a, b); return _mm_add_epi32(_mm_and_si128(_mm_cmpgt_epi32(m0, ret), m2), ret); } __attribute__((target("avx2"))) __attribute__((always_inline)) __m256i my256_mullo_epu32(const __m256i &a, const __m256i &b) { return _mm256_mullo_epi32(a, b); } __attribute__((target("avx2"))) __attribute__((always_inline)) __m256i my256_mulhi_epu32(const __m256i &a, const __m256i &b) { __m256i a13 = _mm256_shuffle_epi32(a, 0xF5); __m256i b13 = _mm256_shuffle_epi32(b, 0xF5); __m256i prod02 = _mm256_mul_epu32(a, b); __m256i prod13 = _mm256_mul_epu32(a13, b13); __m256i prod = _mm256_unpackhi_epi64(_mm256_unpacklo_epi32(prod02, prod13), _mm256_unpackhi_epi32(prod02, prod13)); return prod; } __attribute__((target("avx2"))) __attribute__((always_inline)) __m256i montgomery_mul_256(const __m256i &a, const __m256i &b, const __m256i &r, const __m256i &m1) { return _mm256_sub_epi32( _mm256_add_epi32(my256_mulhi_epu32(a, b), m1), my256_mulhi_epu32(my256_mullo_epu32(my256_mullo_epu32(a, b), r), m1)); } __attribute__((target("avx2"))) __attribute__((always_inline)) __m256i montgomery_add_256(const __m256i &a, const __m256i &b, const __m256i &m2, const __m256i &m0) { __m256i ret = _mm256_sub_epi32(_mm256_add_epi32(a, b), m2); return _mm256_add_epi32(_mm256_and_si256(_mm256_cmpgt_epi32(m0, ret), m2), ret); } __attribute__((target("avx2"))) __attribute__((always_inline)) __m256i montgomery_sub_256(const __m256i &a, const __m256i &b, const __m256i &m2, const __m256i &m0) { __m256i ret = _mm256_sub_epi32(a, b); return _mm256_add_epi32(_mm256_and_si256(_mm256_cmpgt_epi32(m0, ret), m2), ret); } constexpr int SZ = 1 << 19; uint32_t buf1_[SZ * 2] __attribute__((aligned(64))); uint32_t buf2_[SZ * 2] __attribute__((aligned(64))); template struct NTT { static constexpr uint32_t get_pr() { uint32_t mod = mint::get_mod(); using u64 = uint64_t; u64 ds[32] = {}; int idx = 0; u64 m = mod - 1; for (u64 i = 2; i * i <= m; ++i) { if (m % i == 0) { ds[idx++] = i; while (m % i == 0) m /= i; } } if (m != 1) ds[idx++] = m; uint32_t pr = 2; while (1) { int flg = 1; for (int i = 0; i < idx; ++i) { u64 a = pr, b = (mod - 1) / ds[i], r = 1; while (b) { if (b & 1) r = r * a % mod; a = a * a % mod; b >>= 1; } if (r == 1) { flg = 0; break; } } if (flg == 1) break; ++pr; } return pr; }; static constexpr uint32_t mod = mint::get_mod(); static constexpr uint32_t pr = get_pr(); static constexpr int level = __builtin_ctzll(mod - 1); mint dw[level], dy[level]; mint *buf1, *buf2; constexpr NTT() { setwy(level); buf1 = reinterpret_cast(::buf1_); buf2 = reinterpret_cast(::buf2_); } constexpr void setwy(int k) { mint w[level], y[level]; w[k - 1] = mint(pr).pow((mod - 1) / (1 << k)); y[k - 1] = w[k - 1].inverse(); for (int i = k - 2; i > 0; --i) w[i] = w[i + 1] * w[i + 1], y[i] = y[i + 1] * y[i + 1]; dw[0] = dy[0] = w[1] * w[1]; dw[1] = w[1], dy[1] = y[1], dw[2] = w[2], dy[2] = y[2]; for (int i = 3; i < k; ++i) { dw[i] = dw[i - 1] * y[i - 2] * w[i]; dy[i] = dy[i - 1] * w[i - 2] * y[i]; } } __attribute__((target("avx2"))) void ntt(mint *a, int n) { int k = n ? __builtin_ctz(n) : 0; if (k == 0) return; if (k == 1) { mint a1 = a[1]; a[1] = a[0] - a[1]; a[0] = a[0] + a1; return; } if (k & 1) { int v = 1 << (k - 1); if (v < 8) { for (int j = 0; j < v; ++j) { mint ajv = a[j + v]; a[j + v] = a[j] - ajv; a[j] += ajv; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m2 = _mm256_set1_epi32(mod + mod); int j0 = 0; int j1 = v; for (; j0 < v; j0 += 8, j1 += 8) { __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); __m256i naj = montgomery_add_256(T0, T1, m2, m0); __m256i najv = montgomery_sub_256(T0, T1, m2, m0); _mm256_storeu_si256((__m256i *)(a + j0), naj); _mm256_storeu_si256((__m256i *)(a + j1), najv); } } } int u = 1 << (2 + (k & 1)); int v = 1 << (k - 2 - (k & 1)); mint one = mint(1); mint imag = dw[1]; while (v) { if (v == 1) { mint ww = one, xx = one, wx = one; for (int jh = 0; jh < u;) { ww = xx * xx, wx = ww * xx; mint t0 = a[jh + 0], t1 = a[jh + 1] * xx; mint t2 = a[jh + 2] * ww, t3 = a[jh + 3] * wx; mint t0p2 = t0 + t2, t1p3 = t1 + t3; mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag; a[jh + 0] = t0p2 + t1p3, a[jh + 1] = t0p2 - t1p3; a[jh + 2] = t0m2 + t1m3, a[jh + 3] = t0m2 - t1m3; xx *= dw[__builtin_ctz((jh += 4))]; } } else if (v == 4) { const __m128i m0 = _mm_set1_epi32(0); const __m128i m1 = _mm_set1_epi32(mod); const __m128i m2 = _mm_set1_epi32(mod + mod); const __m128i r = _mm_set1_epi32(mint::r); const __m128i Imag = _mm_set1_epi32(imag.a); mint ww = one, xx = one, wx = one; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = j1 + v; int j3 = j2 + v; int je = v; for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i T0P2 = montgomery_add_128(T0, T2, m2, m0); const __m128i T1P3 = montgomery_add_128(T1, T3, m2, m0); const __m128i T0M2 = montgomery_sub_128(T0, T2, m2, m0); const __m128i T1M3 = montgomery_mul_128( montgomery_sub_128(T1, T3, m2, m0), Imag, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_sub_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j2), montgomery_add_128(T0M2, T1M3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j3), montgomery_sub_128(T0M2, T1M3, m2, m0)); } } else { ww = xx * xx, wx = ww * xx; const __m128i WW = _mm_set1_epi32(ww.a); const __m128i WX = _mm_set1_epi32(wx.a); const __m128i XX = _mm_set1_epi32(xx.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i MT1 = montgomery_mul_128(T1, XX, r, m1); const __m128i MT2 = montgomery_mul_128(T2, WW, r, m1); const __m128i MT3 = montgomery_mul_128(T3, WX, r, m1); const __m128i T0P2 = montgomery_add_128(T0, MT2, m2, m0); const __m128i T1P3 = montgomery_add_128(MT1, MT3, m2, m0); const __m128i T0M2 = montgomery_sub_128(T0, MT2, m2, m0); const __m128i T1M3 = montgomery_mul_128( montgomery_sub_128(MT1, MT3, m2, m0), Imag, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_sub_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j2), montgomery_add_128(T0M2, T1M3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j3), montgomery_sub_128(T0M2, T1M3, m2, m0)); } } xx *= dw[__builtin_ctz((jh += 4))]; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m1 = _mm256_set1_epi32(mod); const __m256i m2 = _mm256_set1_epi32(mod + mod); const __m256i r = _mm256_set1_epi32(mint::r); const __m256i Imag = _mm256_set1_epi32(imag.a); mint ww = one, xx = one, wx = one; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = j1 + v; int j3 = j2 + v; int je = v; for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i T0P2 = montgomery_add_256(T0, T2, m2, m0); const __m256i T1P3 = montgomery_add_256(T1, T3, m2, m0); const __m256i T0M2 = montgomery_sub_256(T0, T2, m2, m0); const __m256i T1M3 = montgomery_mul_256( montgomery_sub_256(T1, T3, m2, m0), Imag, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_sub_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j2), montgomery_add_256(T0M2, T1M3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j3), montgomery_sub_256(T0M2, T1M3, m2, m0)); } } else { ww = xx * xx, wx = ww * xx; const __m256i WW = _mm256_set1_epi32(ww.a); const __m256i WX = _mm256_set1_epi32(wx.a); const __m256i XX = _mm256_set1_epi32(xx.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i MT1 = montgomery_mul_256(T1, XX, r, m1); const __m256i MT2 = montgomery_mul_256(T2, WW, r, m1); const __m256i MT3 = montgomery_mul_256(T3, WX, r, m1); const __m256i T0P2 = montgomery_add_256(T0, MT2, m2, m0); const __m256i T1P3 = montgomery_add_256(MT1, MT3, m2, m0); const __m256i T0M2 = montgomery_sub_256(T0, MT2, m2, m0); const __m256i T1M3 = montgomery_mul_256( montgomery_sub_256(MT1, MT3, m2, m0), Imag, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_sub_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j2), montgomery_add_256(T0M2, T1M3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j3), montgomery_sub_256(T0M2, T1M3, m2, m0)); } } xx *= dw[__builtin_ctz((jh += 4))]; } } u <<= 2; v >>= 2; } } __attribute__((target("avx2"))) void intt(mint *a, int n, int normalize = true) { int k = n ? __builtin_ctz(n) : 0; if (k == 0) return; if (k == 1) { mint a1 = a[1]; a[1] = a[0] - a[1]; a[0] = a[0] + a1; if (normalize) { a[0] *= mint(2).inverse(); a[1] *= mint(2).inverse(); } return; } int u = 1 << (k - 2); int v = 1; mint one = mint(1); mint imag = dy[1]; while (u) { if (v == 1) { mint ww = one, xx = one, yy = one; u <<= 2; for (int jh = 0; jh < u;) { ww = xx * xx, yy = xx * imag; mint t0 = a[jh + 0], t1 = a[jh + 1]; mint t2 = a[jh + 2], t3 = a[jh + 3]; mint t0p1 = t0 + t1, t2p3 = t2 + t3; mint t0m1 = (t0 - t1) * xx, t2m3 = (t2 - t3) * yy; a[jh + 0] = t0p1 + t2p3, a[jh + 2] = (t0p1 - t2p3) * ww; a[jh + 1] = t0m1 + t2m3, a[jh + 3] = (t0m1 - t2m3) * ww; xx *= dy[__builtin_ctz(jh += 4)]; } } else if (v == 4) { const __m128i m0 = _mm_set1_epi32(0); const __m128i m1 = _mm_set1_epi32(mod); const __m128i m2 = _mm_set1_epi32(mod + mod); const __m128i r = _mm_set1_epi32(mint::r); const __m128i Imag = _mm_set1_epi32(imag.a); mint ww = one, xx = one, yy = one; u <<= 2; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = v + v; int j3 = j2 + v; for (; j0 < v; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i T0P1 = montgomery_add_128(T0, T1, m2, m0); const __m128i T2P3 = montgomery_add_128(T2, T3, m2, m0); const __m128i T0M1 = montgomery_sub_128(T0, T1, m2, m0); const __m128i T2M3 = montgomery_mul_128( montgomery_sub_128(T2, T3, m2, m0), Imag, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P1, T2P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j2), montgomery_sub_128(T0P1, T2P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_add_128(T0M1, T2M3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j3), montgomery_sub_128(T0M1, T2M3, m2, m0)); } } else { ww = xx * xx, yy = xx * imag; const __m128i WW = _mm_set1_epi32(ww.a); const __m128i XX = _mm_set1_epi32(xx.a); const __m128i YY = _mm_set1_epi32(yy.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i T0P1 = montgomery_add_128(T0, T1, m2, m0); const __m128i T2P3 = montgomery_add_128(T2, T3, m2, m0); const __m128i T0M1 = montgomery_mul_128( montgomery_sub_128(T0, T1, m2, m0), XX, r, m1); __m128i T2M3 = montgomery_mul_128( montgomery_sub_128(T2, T3, m2, m0), YY, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P1, T2P3, m2, m0)); _mm_storeu_si128( (__m128i *)(a + j2), montgomery_mul_128(montgomery_sub_128(T0P1, T2P3, m2, m0), WW, r, m1)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_add_128(T0M1, T2M3, m2, m0)); _mm_storeu_si128( (__m128i *)(a + j3), montgomery_mul_128(montgomery_sub_128(T0M1, T2M3, m2, m0), WW, r, m1)); } } xx *= dy[__builtin_ctz(jh += 4)]; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m1 = _mm256_set1_epi32(mod); const __m256i m2 = _mm256_set1_epi32(mod + mod); const __m256i r = _mm256_set1_epi32(mint::r); const __m256i Imag = _mm256_set1_epi32(imag.a); mint ww = one, xx = one, yy = one; u <<= 2; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = v + v; int j3 = j2 + v; for (; j0 < v; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i T0P1 = montgomery_add_256(T0, T1, m2, m0); const __m256i T2P3 = montgomery_add_256(T2, T3, m2, m0); const __m256i T0M1 = montgomery_sub_256(T0, T1, m2, m0); const __m256i T2M3 = montgomery_mul_256( montgomery_sub_256(T2, T3, m2, m0), Imag, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P1, T2P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j2), montgomery_sub_256(T0P1, T2P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_add_256(T0M1, T2M3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j3), montgomery_sub_256(T0M1, T2M3, m2, m0)); } } else { ww = xx * xx, yy = xx * imag; const __m256i WW = _mm256_set1_epi32(ww.a); const __m256i XX = _mm256_set1_epi32(xx.a); const __m256i YY = _mm256_set1_epi32(yy.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i T0P1 = montgomery_add_256(T0, T1, m2, m0); const __m256i T2P3 = montgomery_add_256(T2, T3, m2, m0); const __m256i T0M1 = montgomery_mul_256( montgomery_sub_256(T0, T1, m2, m0), XX, r, m1); const __m256i T2M3 = montgomery_mul_256( montgomery_sub_256(T2, T3, m2, m0), YY, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P1, T2P3, m2, m0)); _mm256_storeu_si256( (__m256i *)(a + j2), montgomery_mul_256(montgomery_sub_256(T0P1, T2P3, m2, m0), WW, r, m1)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_add_256(T0M1, T2M3, m2, m0)); _mm256_storeu_si256( (__m256i *)(a + j3), montgomery_mul_256(montgomery_sub_256(T0M1, T2M3, m2, m0), WW, r, m1)); } } xx *= dy[__builtin_ctz(jh += 4)]; } } u >>= 4; v <<= 2; } if (k & 1) { v = 1 << (k - 1); if (v < 8) { for (int j = 0; j < v; ++j) { mint ajv = a[j] - a[j + v]; a[j] += a[j + v]; a[j + v] = ajv; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m2 = _mm256_set1_epi32(mod + mod); int j0 = 0; int j1 = v; for (; j0 < v; j0 += 8, j1 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); __m256i naj = montgomery_add_256(T0, T1, m2, m0); __m256i najv = montgomery_sub_256(T0, T1, m2, m0); _mm256_storeu_si256((__m256i *)(a + j0), naj); _mm256_storeu_si256((__m256i *)(a + j1), najv); } } } if (normalize) { mint invn = mint(n).inverse(); for (int i = 0; i < n; i++) a[i] *= invn; } } __attribute__((target("avx2"))) void inplace_multiply( int l1, int l2, int zero_padding = true) { int l = l1 + l2 - 1; int M = 4; while (M < l) M <<= 1; if (zero_padding) { for (int i = l1; i < M; i++) buf1_[i] = 0; for (int i = l2; i < M; i++) buf2_[i] = 0; } const __m256i m0 = _mm256_set1_epi32(0); const __m256i m1 = _mm256_set1_epi32(mod); const __m256i r = _mm256_set1_epi32(mint::r); const __m256i N2 = _mm256_set1_epi32(mint::n2); for (int i = 0; i < l1; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(buf1_ + i)); __m256i b = montgomery_mul_256(a, N2, r, m1); _mm256_storeu_si256((__m256i *)(buf1_ + i), b); } for (int i = 0; i < l2; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(buf2_ + i)); __m256i b = montgomery_mul_256(a, N2, r, m1); _mm256_storeu_si256((__m256i *)(buf2_ + i), b); } ntt(buf1, M); ntt(buf2, M); for (int i = 0; i < M; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(buf1_ + i)); __m256i b = _mm256_loadu_si256((__m256i *)(buf2_ + i)); __m256i c = montgomery_mul_256(a, b, r, m1); _mm256_storeu_si256((__m256i *)(buf1_ + i), c); } intt(buf1, M, false); const __m256i INVM = _mm256_set1_epi32((mint(M).inverse()).a); for (int i = 0; i < l; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(buf1_ + i)); __m256i b = montgomery_mul_256(a, INVM, r, m1); __m256i c = my256_mulhi_epu32(my256_mullo_epu32(b, r), m1); __m256i d = _mm256_and_si256(_mm256_cmpgt_epi32(c, m0), m1); __m256i e = _mm256_sub_epi32(d, c); _mm256_storeu_si256((__m256i *)(buf1_ + i), e); } } void ntt(vector &a) { int M = (int)a.size(); for (int i = 0; i < M; i++) buf1[i].a = a[i].a; ntt(buf1, M); for (int i = 0; i < M; i++) a[i].a = buf1[i].a; } void intt(vector &a) { int M = (int)a.size(); for (int i = 0; i < M; i++) buf1[i].a = a[i].a; intt(buf1, M, true); for (int i = 0; i < M; i++) a[i].a = buf1[i].a; } vector multiply(const vector &a, const vector &b) { if (a.size() == 0 && b.size() == 0) return vector{}; int l = a.size() + b.size() - 1; if (min(a.size(), b.size()) <= 40) { vector s(l); for (int i = 0; i < (int)a.size(); ++i) for (int j = 0; j < (int)b.size(); ++j) s[i + j] += a[i] * b[j]; return s; } int M = 4; while (M < l) M <<= 1; for (int i = 0; i < (int)a.size(); ++i) buf1[i].a = a[i].a; for (int i = (int)a.size(); i < M; ++i) buf1[i].a = 0; for (int i = 0; i < (int)b.size(); ++i) buf2[i].a = b[i].a; for (int i = (int)b.size(); i < M; ++i) buf2[i].a = 0; ntt(buf1, M); ntt(buf2, M); for (int i = 0; i < M; ++i) buf1[i].a = mint::reduce(uint64_t(buf1[i].a) * buf2[i].a); intt(buf1, M, false); vector s(l); mint invm = mint(M).inverse(); for (int i = 0; i < l; ++i) s[i] = buf1[i] * invm; return s; } void ntt_doubling(vector &a) { int M = (int)a.size(); for (int i = 0; i < M; i++) buf1[i].a = a[i].a; intt(buf1, M); mint r = 1, zeta = mint(pr).pow((mint::get_mod() - 1) / (M << 1)); for (int i = 0; i < M; i++) buf1[i] *= r, r *= zeta; ntt(buf1, M); a.resize(2 * M); for (int i = 0; i < M; i++) a[M + i].a = buf1[i].a; } }; template void FormalPowerSeries::set_fft() { if (!ntt_ptr) ntt_ptr = new NTT; } template FormalPowerSeries &FormalPowerSeries::operator*=( const FormalPowerSeries &r) { if (this->empty() || r.empty()) { this->clear(); return *this; } set_fft(); auto ret = static_cast *>(ntt_ptr)->multiply(*this, r); return *this = FormalPowerSeries(ret.begin(), ret.end()); } template void FormalPowerSeries::ntt() { set_fft(); static_cast *>(ntt_ptr)->ntt(*this); } template void FormalPowerSeries::intt() { set_fft(); static_cast *>(ntt_ptr)->intt(*this); } template void FormalPowerSeries::ntt_doubling() { set_fft(); static_cast *>(ntt_ptr)->ntt_doubling(*this); } template int FormalPowerSeries::ntt_pr() { set_fft(); return static_cast *>(ntt_ptr)->pr; } template FormalPowerSeries FormalPowerSeries::inv(int deg) const { assert((*this)[0] != mint(0)); if (deg == -1) deg = (int)this->size(); FormalPowerSeries res(deg); res[0] = {mint(1) / (*this)[0]}; for (int d = 1; d < deg; d <<= 1) { FormalPowerSeries f(2 * d), g(2 * d); for (int j = 0; j < min((int)this->size(), 2 * d); j++) f[j] = (*this)[j]; for (int j = 0; j < d; j++) g[j] = res[j]; f.ntt(); g.ntt(); for (int j = 0; j < 2 * d; j++) f[j] *= g[j]; f.intt(); for (int j = 0; j < d; j++) f[j] = 0; f.ntt(); for (int j = 0; j < 2 * d; j++) f[j] *= g[j]; f.intt(); for (int j = d; j < min(2 * d, deg); j++) res[j] = -f[j]; } return res.pre(deg); } namespace workspace { using mint = LazyMontgomeryModInt<998244353>; using poly = FormalPowerSeries; void main() { // start here! int n; cin >> n; poly p(n); mint f = 1; for (int i = 0; i < n; ++i) { p[i] = f *= i + 1; } p = p.inv(); cout << mint{1} - std::accumulate(begin(p), end(p), mint{0}) << "\n"; } } // namespace workspace