Submit Info #41725

Problem Lang User Status Time Memory
Log of Formal Power Series cpp jell AC 326 ms 16.43 MiB

ケース詳細
Name Status Time Memory
example_00 AC 2 ms 0.67 MiB
max_all_zero_00 AC 165 ms 8.61 MiB
max_random_00 AC 326 ms 16.41 MiB
max_random_01 AC 325 ms 16.40 MiB
max_random_02 AC 324 ms 16.41 MiB
max_random_03 AC 323 ms 16.43 MiB
max_random_04 AC 324 ms 16.41 MiB
near_262144_00 AC 160 ms 8.78 MiB
near_262144_01 AC 160 ms 8.84 MiB
near_262144_02 AC 161 ms 8.78 MiB
random_00 AC 296 ms 15.16 MiB
random_01 AC 314 ms 16.03 MiB
random_02 AC 37 ms 2.62 MiB
random_03 AC 305 ms 15.53 MiB
random_04 AC 264 ms 13.91 MiB
small_degree_00 AC 1 ms 0.67 MiB
small_degree_01 AC 1 ms 0.65 MiB
small_degree_02 AC 1 ms 0.59 MiB
small_degree_03 AC 1 ms 0.58 MiB
small_degree_04 AC 2 ms 0.67 MiB
small_degree_05 AC 1 ms 0.65 MiB
small_degree_06 AC 0 ms 0.62 MiB
small_degree_07 AC 1 ms 0.67 MiB
small_degree_08 AC 1 ms 0.67 MiB
small_degree_09 AC 1 ms 0.68 MiB

#line 1 "atcoder-workspace/2.cc" // #undef _GLIBCXX_DEBUG // #define NDEBUG #include <bits/extc++.h> #line 2 "Library/lib/alias" /** * @file alias * @brief Alias */ #line 13 "Library/lib/alias" #line 2 "Library/lib/bit" #if __cplusplus > 201703L #include <bit> #else #ifndef _GLIBCXX_BIT #define _GLIBCXX_BIT 1 #include <limits> #include <type_traits> namespace std { template <typename _Tp> 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 <typename _Tp> 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 <typename _Tp> 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<unsigned long long>::digits; constexpr auto _Nd_ul = numeric_limits<unsigned long>::digits; constexpr auto _Nd_u = numeric_limits<unsigned>::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<unsigned long long>::max(); unsigned long long __low = __x & __max_ull; return (_Nd - _Nd_ull) + __builtin_clzll(__low); } } template <typename _Tp> 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 <typename _Tp> 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<unsigned long long>::digits; constexpr auto _Nd_ul = numeric_limits<unsigned long>::digits; constexpr auto _Nd_u = numeric_limits<unsigned>::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<unsigned long long>::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 <typename _Tp> constexpr int __countr_one(_Tp __x) noexcept { if (__x == numeric_limits<_Tp>::max()) return numeric_limits<_Tp>::digits; return __countr_zero((_Tp)~__x); } template <typename _Tp> constexpr int __popcount(_Tp __x) noexcept { constexpr auto _Nd = numeric_limits<_Tp>::digits; if (__x == 0) return 0; constexpr auto _Nd_ull = numeric_limits<unsigned long long>::digits; constexpr auto _Nd_ul = numeric_limits<unsigned long>::digits; constexpr auto _Nd_u = numeric_limits<unsigned>::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<unsigned long long>::max(); unsigned long long __low = __x & __max_ull; unsigned long long __high = __x >> _Nd_ull; return __builtin_popcountll(__low) + __builtin_popcountll(__high); } } template <typename _Tp> constexpr bool __has_single_bit(_Tp __x) noexcept { return __popcount(__x) == 1; } template <typename _Tp> 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 <typename _Tp> 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 <typename _Tp> constexpr _Tp __bit_width(_Tp __x) noexcept { constexpr auto _Nd = numeric_limits<_Tp>::digits; return _Nd - __countl_zero(__x); } } // namespace std #endif #endif #line 2 "Library/lib/limits" #line 4 "Library/lib/limits" namespace workspace { template <class _Tp> struct numeric_limits : std::numeric_limits<_Tp> {}; #ifdef __SIZEOF_INT128__ template <> struct numeric_limits<__uint128_t> { constexpr static __uint128_t max() { return ~__uint128_t(0); } constexpr static __uint128_t min() { return 0; } }; template <> struct numeric_limits<__int128_t> { constexpr static __int128_t max() { return numeric_limits<__uint128_t>::max() >> 1; } constexpr static __int128_t min() { return -max() - 1; } }; #endif } // namespace workspace #line 16 "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 <class T, class Comp = less<T>> using priority_queue = std::priority_queue<T, vector<T>, Comp>; template <class T> using stack = std::stack<T, vector<T>>; template <class _T1, class _T2, typename = decltype(std::declval<const _T2 &>() < std::declval<const _T1 &>())> constexpr typename std::conditional<std::is_same<_T1, _T2>::value, const _T1 &, typename std::common_type<_T1, _T2>::type>::type min(const _T1 &__x, const _T2 &__y) noexcept { return __y < __x ? __y : __x; } template <class _T1, class _T2, class _Compare, typename = decltype(std::declval<_Compare>()( std::declval<const _T2 &>(), std::declval<const _T1 &>()))> constexpr typename std::conditional<std::is_same<_T1, _T2>::value, const _T1 &, typename std::common_type<_T1, _T2>::type>::type min(const _T1 &__x, const _T2 &__y, _Compare __comp) noexcept { return __comp(__y, __x) ? __y : __x; } template <class _Tp, typename = decltype(std::declval<const _Tp &>() < std::declval<const _Tp &>())> constexpr _Tp min(std::initializer_list<_Tp> __x) noexcept { return *std::min_element(__x.begin(), __x.end()); } template <class _Tp, class _Compare, typename = decltype(std::declval<_Compare>()( std::declval<const _Tp &>(), std::declval<const _Tp &>()))> constexpr _Tp min(std::initializer_list<_Tp> __x, _Compare __comp) noexcept { return *std::min_element(__x.begin(), __x.end(), __comp); } template <class _T1, class _T2, typename = decltype(std::declval<const _T1 &>() < std::declval<const _T2 &>())> constexpr typename std::conditional<std::is_same<_T1, _T2>::value, const _T1 &, typename std::common_type<_T1, _T2>::type>::type max(const _T1 &__x, const _T2 &__y) noexcept { return __x < __y ? __y : __x; } template <class _T1, class _T2, class _Compare, typename = decltype(std::declval<_Compare>()( std::declval<const _T1 &>(), std::declval<const _T2 &>()))> constexpr typename std::conditional<std::is_same<_T1, _T2>::value, const _T1 &, typename std::common_type<_T1, _T2>::type>::type max(const _T1 &__x, const _T2 &__y, _Compare __comp) noexcept { return __comp(__x, __y) ? __y : __x; } template <class _Tp, typename = decltype(std::declval<const _Tp &>() < std::declval<const _Tp &>())> constexpr _Tp max(std::initializer_list<_Tp> __x) noexcept { return *std::max_element(__x.begin(), __x.end()); } template <class _Tp, class _Compare, typename = decltype(std::declval<_Compare>()( std::declval<const _Tp &>(), std::declval<const _Tp &>()))> constexpr _Tp max(std::initializer_list<_Tp> __x, _Compare __comp) noexcept { return *std::max_element(__x.begin(), __x.end(), __comp); } template <typename _Tp> constexpr _Tp __bsf(_Tp __x) noexcept { return std::__countr_zero(__x); } template <typename _Tp> constexpr _Tp __bsr(_Tp __x) noexcept { return std::__bit_width(__x) - 1; } } // namespace workspace #line 2 "Library/lib/cxx20" /* * @file cxx20 * @brief C++20 Features */ #line 9 "Library/lib/cxx20" #if __cplusplus <= 201703L #include <algorithm> #include <vector> namespace std { /* * @fn erase_if * @brief Erase the elements of a container that do not satisfy the condition. * @param __cont Container. * @param __pred Predicate. * @return Number of the erased elements. */ template <typename _Tp, typename _Alloc, typename _Predicate> inline typename vector<_Tp, _Alloc>::size_type erase_if( vector<_Tp, _Alloc>& __cont, _Predicate __pred) { const auto __osz = __cont.size(); __cont.erase(remove_if(__cont.begin(), __cont.end(), __pred), __cont.end()); return __osz - __cont.size(); } /* * @fn erase * @brief Erase the elements of a container that are equal to the given value. * @param __cont Container. * @param __value Value. * @return Number of the erased elements. */ template <typename _Tp, typename _Alloc, typename _Up> inline typename vector<_Tp, _Alloc>::size_type erase( vector<_Tp, _Alloc>& __cont, const _Up& __value) { const auto __osz = __cont.size(); __cont.erase(remove(__cont.begin(), __cont.end(), __value), __cont.end()); return __osz - __cont.size(); } } // namespace std #endif #line 2 "Library/lib/direct" /* * @file direct * @brief Pragma Directive */ #ifdef ONLINE_JUDGE #pragma GCC optimize("O3") #pragma GCC target("avx,avx2") #pragma GCC optimize("unroll-loops") #endif #line 2 "Library/src/opt/binary_search.hpp" /** * @file binary_search.hpp * @brief Binary Search */ #line 12 "Library/src/opt/binary_search.hpp" namespace workspace { namespace _impl { template <class _Iter, class _Pred> _Iter discrete_binary_search(_Iter __first, _Iter __last, _Pred &&__pred, bool __crit) { if (__first < __last) while (1 < __last - __first) { _Iter __mid = __first + ((__last - __first) >> 1); (__pred(__mid) == __crit ? __first : __last) = __mid; } else while (1 < __first - __last) { _Iter __mid = __last + ((__first - __last) >> 1); (__pred(__mid) == __crit ? __first : __last) = __mid; } return __first; } template <class _Real, class _Pred> _Real continuous_binary_search(_Real __first, _Real __last, _Real __eps, _Pred &&__pred, bool __crit) { if (__first < __last) for (auto __loop = std::numeric_limits<_Real>::digits; __first + __eps < __last && __loop; --__loop) { _Real __mid = (__first + __last) / 2; (__pred(__mid) == __crit ? __first : __last) = __mid; } else for (auto __loop = std::numeric_limits<_Real>::digits; __last + __eps < __first && __loop; --__loop) { _Real __mid = (__first + __last) / 2; (__pred(__mid) == __crit ? __first : __last) = __mid; } return __first; } } // namespace _impl /** * @brief Binary search on a discrete range. * * @param __first Inclusive endpoint * @param __last Exclusive endpoint * @param __pred Predicate * @return Return __x s.t. `__pred(__x) == __pred(__first)` and `__pred(__x) != * __pred(__x +/- 1)`. The sign is `+` iff `__first < __last`. */ template <class _Iter, class _Pred> typename std::enable_if< std::is_convertible<decltype(std::declval<_Pred>()(std::declval<_Iter>())), bool>::value, _Iter>::type binary_search(_Iter __first, _Iter __last, _Pred &&__pred) { assert(__first != __last); return _impl::discrete_binary_search( __first, __last, std::forward<_Pred>(__pred), __pred(__first)); } /** * @brief Binary search on the real number line. * @param __first __pred(__first) is true * @param __last __pred(__last) is false * @param eps the error tolerance * @param __pred the predicate * @return the boundary point */ template <class _Real, class _Pred> typename std::enable_if< std::is_convertible<decltype(std::declval<_Pred>()(std::declval<_Real>())), bool>::value, _Real>::type binary_search(_Real __first, _Real __last, _Real __eps, _Pred &&__pred) { assert(__first != __last); return _impl::continuous_binary_search( __first, __last, __eps, std::forward<_Pred>(__pred), __pred(__first)); } /** * @brief Parallel binary search on discrete ranges. * @param ends a vector of pairs; __pred(first) is true, __pred(second) is false * @param __pred the predicate * @return the closest points to (second) where __pred is true */ template <class Array, class _Iter = typename std::decay< decltype(std::get<0>(std::declval<Array>()[0]))>::type, class _Pred> typename std::enable_if< std::is_convertible< decltype(std::declval<_Pred>()(std::declval<std::vector<_Iter>>())[0]), bool>::value, std::vector<_Iter>>::type parallel_binary_search(Array ends, _Pred __pred) { std::vector<_Iter> mids(std::size(ends)); for (;;) { bool all_found = true; for (size_t i{}; i != std::size(ends); ++i) { const _Iter &__first = std::get<0>(ends[i]); const _Iter &__last = std::get<1>(ends[i]); const _Iter mid( __first + typename std::make_signed<decltype(__last - __first)>::type( __last - __first) / 2); if (mids[i] != mid) { all_found = false; mids[i] = mid; } } if (all_found) break; const auto res = __pred(mids); for (size_t i{}; i != std::size(ends); ++i) { (res[i] ? std::get<0>(ends[i]) : std::get<1>(ends[i])) = mids[i]; } } return mids; } /** * @brief Parallel binary search on the real number line. * @param ends a vector of pairs; __pred(first) is true, __pred(second) is false * @param eps the error tolerance * @param __pred the predicate * @return the boundary points */ template <class Array, class _Real = typename std::decay< decltype(std::get<0>(std::declval<Array>()[0]))>::type, class _Pred> typename std::enable_if< std::is_convertible< decltype(std::declval<_Pred>()(std::declval<std::vector<_Real>>())[0]), bool>::value, std::vector<_Real>>::type parallel_binary_search(Array ends, const _Real eps, _Pred __pred) { std::vector<_Real> mids(std::size(ends)); for (auto loops = 0; loops != std::numeric_limits<_Real>::digits; ++loops) { bool all_found = true; for (size_t i{}; i != std::size(ends); ++i) { const _Real __first = std::get<0>(ends[i]); const _Real __last = std::get<1>(ends[i]); if (__first + eps < __last || __last + eps < __first) { all_found = false; mids[i] = (__first + __last) / 2; } } if (all_found) break; const auto res = __pred(mids); for (size_t i{}; i != std::size(ends); ++i) { (res[i] ? std::get<0>(ends[i]) : std::get<1>(ends[i])) = mids[i]; } } return mids; } } // namespace workspace #line 2 "Library/src/opt/exponential_search.hpp" /** * @file exponential_search.hpp * @brief Exponential Search */ #line 9 "Library/src/opt/exponential_search.hpp" namespace workspace { /** * @brief Exponential search on a discrete range. * * @param __first Inclusive endpoint * @param __last Exclusive endpoint * @param __pred Predicate * @return Return __x s.t. `__pred(__x) == __pred(__first)` and `__pred(__x) != * __pred(__x +/- 1)`. The sign is `+` iff `__first < __last`. */ template <class _Iter, class _Pred> typename std::enable_if< std::is_convertible<decltype(std::declval<_Pred>()(std::declval<_Iter>())), bool>::value, _Iter>::type exponential_search(_Iter __first, _Iter __last, _Pred&& __pred) { assert(__first != __last); bool __crit = __pred(__first); if (__first < __last) { for (decltype(__last - __first) __step = 1; __step < __last - __first; __step <<= 1) if (__pred(__first += __step) != __crit) return _impl::discrete_binary_search( __first - __step, __first, std::forward<_Pred>(__pred), __crit); } else { for (decltype(__first - __last) __step = 1; __step < __first - __last; __step <<= 1) if (__pred(__first -= __step) != __crit) return _impl::discrete_binary_search( __first + __step, __first, std::forward<_Pred>(__pred), __crit); } return _impl::discrete_binary_search(__first, __last, std::forward<_Pred>(__pred), __crit); } /** * @brief Exponential search on the real number line. * * @param __first Endpoint * @param __last Endpoint * @param __eps Error tolerance * @param __pred Predicate * @return A boundary point `__x` s.t. `__pred(__x) == __pred(__first)`. */ template <class _Real, class _Pred> typename std::enable_if< std::is_convertible<decltype(std::declval<_Pred>()(std::declval<_Real>())), bool>::value, _Real>::type exponential_search(_Real __first, _Real __last, _Real __eps, _Pred&& __pred) { bool __crit = __pred(__first); if (__first < __last) { for (_Real __step = __eps; __step < __last - __first; __step += __step) if (__pred(__first += __step) != __crit) return _impl::continuous_binary_search(__first - __step, __first, __eps, std::forward<_Pred>(__pred), __crit); } else { for (_Real __step = __eps; __step < __first - __last; __step += __step) if (__pred(__first -= __step) != __crit) return _impl::continuous_binary_search(__first + __step, __first, __eps, std::forward<_Pred>(__pred), __crit); } return _impl::continuous_binary_search(__first, __last, __eps, std::forward<_Pred>(__pred), __crit); } } // namespace workspace #line 2 "Library/src/opt/trinary_search.hpp" /* * @file trinary_search.hpp * @brief Trinary Search */ #line 9 "Library/src/opt/trinary_search.hpp" #include <type_traits> namespace workspace { /* * @brief Trinary search on discrete range. * @param first Left end, inclusive * @param last Right end, exclusive * @param comp Compare function * @return Local minimal point. */ template <class Iter, class Comp> typename std::enable_if< std::is_convertible<decltype(std::declval<Comp>()(std::declval<Iter>(), std::declval<Iter>())), bool>::value, Iter>::type trinary_search(Iter first, Iter last, Comp comp) { assert(first < last); typename std::make_signed<decltype(last - first)>::type dist(last - first); while (2 < dist) { Iter left(first + dist / 3), right(first + dist * 2 / 3); if (comp(left, right)) last = right, dist = (dist + dist) / 3; else first = left, dist -= dist / 3; } if (1 < dist && comp(first + 1, first)) ++first; return first; } /* * @brief Trinary search on discrete range. * @param first Left end, inclusive * @param last Right end, exclusive * @param func Function * @return Local minimal point. */ template <class Iter, class Func> typename std::enable_if< std::is_same<decltype(std::declval<Func>()(std::declval<Iter>()), nullptr), std::nullptr_t>::value, Iter>::type trinary_search(Iter const &first, Iter const &last, Func func) { return trinary_search(first, last, [&](Iter const &__i, Iter const &__j) { return func(__i) < func(__j); }); } /* * @brief Trinary search on the real number line. * @param first Left end * @param last Right end * @param eps Error tolerance * @param comp Compare function * @return Local minimal point. */ template <class Real, class Comp> typename std::enable_if< std::is_convertible<decltype(std::declval<Comp>()(std::declval<Real>(), std::declval<Real>())), bool>::value, Real>::type trinary_search(Real first, Real last, Real const &eps, Comp comp) { assert(first < last); while (eps < last - first) { Real left{(first * 2 + last) / 3}, right{(first + last * 2) / 3}; if (comp(left, right)) last = right; else first = left; } return first; } /* * @brief Trinary search on the real number line. * @param first Left end * @param last Right end * @param eps Error tolerance * @param func Function * @return Local minimal point. */ template <class Real, class Func> typename std::enable_if< std::is_same<decltype(std::declval<Func>()(std::declval<Real>()), nullptr), std::nullptr_t>::value, Real>::type trinary_search(Real const &first, Real const &last, Real const &eps, Func func) { return trinary_search( first, last, eps, [&](Real const &__i, Real const &__j) { return func(__i) < func(__j); }); } } // namespace workspace #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<milliseconds>(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 <class Tp> 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/cached.hpp" /** * @file cached.hpp * @brief Cached * @date 2021-02-25 * * */ #line 13 "Library/src/utils/cached.hpp" #line 2 "Library/src/utils/fixed_point.hpp" /** * @file fixed_point.hpp * @brief Fixed Point Combinator */ #line 9 "Library/src/utils/fixed_point.hpp" namespace workspace { /** * @brief Fixed Point Combinator */ template <class _F> class fixed_point { _F __fn; public: /** * @brief Construct a new fixed point object * * @param __fn 1st argument callable with the rest of its arguments. * Return type specified. */ fixed_point(_F __fn) noexcept : __fn(std::forward<_F>(__fn)) {} /** * @brief Apply *this to 1st argument. * @param __args Rest of arguments. */ template <class... _Args> decltype(auto) operator()(_Args &&...__args) const noexcept { return __fn(*this, std::forward<_Args>(__args)...); } }; } // namespace workspace #line 15 "Library/src/utils/cached.hpp" namespace workspace { namespace cached_impl { // Convert keys to tuple. template <class... _Args> struct get_tuple { using type = decltype( std::tuple_cat(std::declval<std::tuple<std::conditional_t< std::is_convertible<std::decay_t<_Args>, _Args>::value, std::decay_t<_Args>, _Args>>>()...)); }; // Associative array. template <class _Value, class... _Keys> struct assoc : std::integral_constant<int, !std::is_void<_Value>::value>, std::conditional_t<std::is_void<_Value>::value, std::set<typename get_tuple<_Keys...>::type>, std::map<typename get_tuple<_Keys...>::type, _Value>> { }; // Non-resursive lambda type. template <class _F, class = void> struct is_recursive : std::false_type {}; // Resursive lambda type. template <class _F> struct is_recursive< _F, std::__void_t<decltype(&_F::template operator()<fixed_point<_F> &>)>> : std::true_type {}; // Recursive ver. template <class _F> class _recursive { template <class...> struct _cache; template <class _G, class _R, class _H, class... _Args> struct _cache<_R (_G::*)(_H, _Args...)> : assoc<_R, _Args...> {}; template <class _G, class _R, class _H, class... _Args> struct _cache<_R (_G::*)(_H, _Args...) const> : assoc<_R, _Args...> {}; public: using cache = _cache<decltype(&_F::template operator()<_recursive<_F> &>)>; _recursive(_F __x) noexcept : __fn(__x), __cptr(new cache) {} /** * @brief Apply `*this` to 1st argument of the lambda. * @param __args Rest of arguments. */ template <class... _Args> decltype(auto) operator()(_Args &&...__args) noexcept { typename cache::key_type __key{__args...}; if constexpr (cache::value) { if (auto __i = __cptr->lower_bound(__key); __i != __cptr->end() && __i->first == __key) return __i->second; else return __cptr ->emplace_hint(__i, std::move(__key), __fn(*this, std::forward<_Args>(__args)...)) ->second; } else if (auto __i = __cptr->lower_bound(__key); __i == __cptr->end() || *__i != __key) __cptr->emplace_hint(__i, std::move(__key)), __fn(*this, std::forward<_Args>(__args)...); } private: _F __fn; std::shared_ptr<cache> __cptr; }; // Non-recursive ver. template <class _F> class _non_recursive { template <class _T, class = void> struct _get_func { using type = _T; }; template <class _T> struct _get_func<_T, std::__void_t<decltype(&_T::operator())>> { using type = decltype(&_T::operator()); }; template <class...> struct _cache; template <class _R, class... _Args> struct _cache<_R(_Args...)> : assoc<_R, _Args...> {}; template <class _R, class... _Args> struct _cache<_R (*)(_Args...)> : assoc<_R, _Args...> {}; template <class _G, class _R, class... _Args> struct _cache<_R (_G::*)(_Args...)> : assoc<_R, _Args...> {}; template <class _G, class _R, class... _Args> struct _cache<_R (_G::*)(_Args...) const> : assoc<_R, _Args...> {}; public: using cache = _cache<typename _get_func<_F>::type>; _non_recursive(_F __x) noexcept : __fn(__x), __cptr(new cache) {} /** * @param __args */ template <class... _Args> decltype(auto) operator()(_Args &&...__args) noexcept { typename cache::key_type __key{__args...}; if constexpr (cache::value) { if (auto __i = __cptr->lower_bound(__key); __i != __cptr->end() && __i->first == __key) return __i->second; else return __cptr ->emplace_hint(__i, std::move(__key), __fn(std::forward<_Args>(__args)...)) ->second; } else if (auto __i = __cptr->lower_bound(__key); __i == __cptr->end() || *__i != __key) __cptr->emplace_hint(__i, std::move(__key)), __fn(std::forward<_Args>(__args)...); } private: _F __fn; std::shared_ptr<cache> __cptr; }; template <class _F> using _cached = std::conditional_t<is_recursive<_F>::value, _recursive<_F>, _non_recursive<_F>>; } // namespace cached_impl /** * @brief Cached caller of function */ template <class _F> class cached : public cached_impl::_cached<_F> { public: /** * @brief Construct a new cached object */ cached() noexcept : cached_impl::_cached<_F>(_F{}) {} /** * @brief Construct a new cached object * * @param __x Function */ cached(_F __x) noexcept : cached_impl::_cached<_F>(__x) {} }; } // namespace workspace #line 2 "Library/src/utils/cat.hpp" /** * @file cat.hpp * @brief Cat */ #line 9 "Library/src/utils/cat.hpp" namespace workspace { /** * @brief Concatenate two sequences. * * @param __c1 * @param __c2 * @return Concatenated sequence. */ template <class _C1, class _C2> constexpr decltype(auto) cat(_C1 &&__c1, _C2 &&__c2) noexcept { auto __c = std::forward<_C1>(__c1); if constexpr (std::is_rvalue_reference<decltype(__c2)>::value) __c.insert(std::end(__c), std::move_iterator(std::begin(__c2)), std::move_iterator(std::end(__c2))); else __c.insert(std::end(__c), std::cbegin(__c2), std::cend(__c2)); return __c; } /** * @return Concatenated sequence. */ template <class _C1, class _C2, class... _Args> constexpr decltype(auto) cat(_C1 &&__c1, _C2 &&__c2, _Args &&...__args) noexcept { return cat(cat(std::forward<_C1>(__c1), std::forward<_C2>(__c2)), std::forward<_Args>(__args)...); } } // 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 { /** * @brief Substitute __y for __x if __y < __x. * @param __x Reference * @param __y Comparison target * @return Whether or not __x is updated. */ template <class _T1, class _T2, typename = decltype(std::declval<_T2>() < std::declval<_T1 &>())> typename std::enable_if<std::is_assignable<_T1 &, _T2>::value, bool>::type chle( _T1 &__x, _T2 &&__y) noexcept { return __y < __x ? __x = std::forward<_T2>(__y), true : false; } /** * @brief Substitute __y for __x if __x < __y. * @param __x Reference * @param __y Comparison target * @return Whether or not __x is updated. */ template <class _T1, class _T2, typename = decltype(std::declval<_T1 &>() < std::declval<_T2>())> typename std::enable_if<std::is_assignable<_T1 &, _T2>::value, bool>::type chge( _T1 &__x, _T2 &&__y) noexcept { return __x < __y ? __x = std::forward<_T2>(__y), true : false; } /** * @brief Substitute __y for __x if __comp(__y, __x) is true. * @param __x Reference * @param __y Comparison target * @param __comp Compare function object * @return Whether or not __x is updated. */ template <class _T1, class _T2, class _Compare, typename = decltype(std::declval<_Compare>()(std::declval<_T2>(), std::declval<_T1 &>()))> typename std::enable_if<std::is_assignable<_T1 &, _T2>::value, bool>::type chle( _T1 &__x, _T2 &&__y, _Compare __comp) noexcept { return __comp(__y, __x) ? __x = std::forward<_T2>(__y), true : false; } /** * @brief Substitute __y for __x if __comp(__x, __y) is true. * @param __x Reference * @param __y Comparison target * @param __comp Compare function object * @return Whether or not __x is updated. */ template <class _T1, class _T2, class _Compare, typename = decltype(std::declval<_Compare>()(std::declval<_T1 &>(), std::declval<_T2>()))> typename std::enable_if<std::is_assignable<_T1 &, _T2>::value, bool>::type chge( _T1 &__x, _T2 &&__y, _Compare __comp) noexcept { return __comp(__x, __y) ? __x = std::forward<_T2>(__y), true : false; } } // namespace workspace #line 6 "Library/lib/utils" // #include "src/utils/grid.hpp" #line 2 "Library/src/utils/hash.hpp" #line 8 "Library/src/utils/hash.hpp" #line 2 "Library/src/utils/sfinae.hpp" /** * @file sfinae.hpp * @brief SFINAE */ #line 11 "Library/src/utils/sfinae.hpp" #ifndef __INT128_DEFINED__ #ifdef __SIZEOF_INT128__ #define __INT128_DEFINED__ 1 #else #define __INT128_DEFINED__ 0 #endif #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 Tp, class... Args> struct variadic_front { using type = Tp; }; template <class... Args> struct variadic_back; template <class Tp> struct variadic_back<Tp> { using type = Tp; }; template <class Tp, class... Args> struct variadic_back<Tp, Args...> { using type = typename variadic_back<Args...>::type; }; template <class type, template <class> class trait> using enable_if_trait_type = typename std::enable_if<trait<type>::value>::type; /** * @brief Return type of subscripting ( @c [] ) access. */ template <class _Tp> using subscripted_type = typename std::decay<decltype(std::declval<_Tp&>()[0])>::type; template <class Container> using element_type = typename std::decay<decltype( *std::begin(std::declval<Container&>()))>::type; template <class _Tp, class = std::nullptr_t> struct has_begin : std::false_type {}; template <class _Tp> struct has_begin<_Tp, decltype(std::begin(std::declval<_Tp>()), nullptr)> : std::true_type {}; template <class _Tp, class = std::nullptr_t> struct has_mod : std::false_type {}; template <class _Tp> struct has_mod<_Tp, decltype(_Tp::mod, nullptr)> : std::true_type {}; template <class _Tp, class = void> struct is_integral_ext : std::false_type {}; template <class _Tp> struct is_integral_ext< _Tp, typename std::enable_if<std::is_integral<_Tp>::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 <class _Tp> constexpr static bool is_integral_ext_v = is_integral_ext<_Tp>::value; #endif template <typename _Tp, typename = void> struct multiplicable_uint { using type = uint_least32_t; }; template <typename _Tp> struct multiplicable_uint< _Tp, typename std::enable_if<(2 < sizeof(_Tp)) && (!__INT128_DEFINED__ || sizeof(_Tp) <= 4)>::type> { using type = uint_least64_t; }; #if __INT128_DEFINED__ template <typename _Tp> struct multiplicable_uint<_Tp, typename std::enable_if<(4 < sizeof(_Tp))>::type> { using type = __uint128_t; }; #endif template <typename _Tp> struct multiplicable_int { using type = typename std::make_signed<typename multiplicable_uint<_Tp>::type>::type; }; } // namespace workspace #line 10 "Library/src/utils/hash.hpp" namespace workspace { template <class T, class = void> struct hash : std::hash<T> {}; template <class _Tp> struct hash<_Tp *> : std::hash<_Tp *> {}; #if __cplusplus >= 201703L template <class Unique_bits_type> struct hash<Unique_bits_type, enable_if_trait_type<Unique_bits_type, std::has_unique_object_representations>> { size_t operator()(uint64_t x) const { static const uint64_t m = std::random_device{}(); x ^= x >> 23; x ^= m; x ^= x >> 47; return x - (x >> 32); } }; #endif template <class Key> size_t hash_combine(const size_t &seed, const Key &key) { return seed ^ (hash<Key>()(key) + 0x9e3779b9 /* + (seed << 6) + (seed >> 2) */); } template <class T1, class T2> struct hash<std::pair<T1, T2>> { size_t operator()(const std::pair<T1, T2> &pair) const { return hash_combine(hash<T1>()(pair.first), pair.second); } }; template <class... T> class hash<std::tuple<T...>> { template <class Tuple, size_t index = std::tuple_size<Tuple>::value - 1> struct tuple_hash { static uint64_t apply(const Tuple &t) { return hash_combine(tuple_hash<Tuple, index - 1>::apply(t), std::get<index>(t)); } }; template <class Tuple> struct tuple_hash<Tuple, size_t(-1)> { static uint64_t apply(const Tuple &t) { return 0; } }; public: uint64_t operator()(const std::tuple<T...> &t) const { return tuple_hash<std::tuple<T...>>::apply(t); } }; template <class hash_table> struct hash_table_wrapper : hash_table { using key_type = typename hash_table::key_type; size_t count(const key_type &key) const { return hash_table::find(key) != hash_table::end(); } template <class... Args> auto emplace(Args &&... args) { return hash_table::insert(typename hash_table::value_type(args...)); } }; template <class Key, class Mapped = __gnu_pbds::null_type> using cc_hash_table = hash_table_wrapper<__gnu_pbds::cc_hash_table<Key, Mapped, hash<Key>>>; template <class Key, class Mapped = __gnu_pbds::null_type> using gp_hash_table = hash_table_wrapper<__gnu_pbds::gp_hash_table<Key, Mapped, hash<Key>>>; template <class Key, class Mapped> using unordered_map = std::unordered_map<Key, Mapped, hash<Key>>; template <class Key> using unordered_set = std::unordered_set<Key, hash<Key>>; } // namespace workspace #line 2 "Library/src/utils/io/istream.hpp" /** * @file istream.hpp * @brief Input Stream */ #include <cxxabi.h> #line 13 "Library/src/utils/io/istream.hpp" #line 15 "Library/src/utils/io/istream.hpp" namespace workspace { namespace internal { template <class Tp, typename = std::nullptr_t> struct istream_helper { istream_helper(std::istream &is, Tp &x) { if constexpr (has_begin<Tp>::value) for (auto &&e : x) istream_helper<typename std::decay<decltype(e)>::type>(is, e); else static_assert(has_begin<Tp>::value, "istream unsupported type."); } }; template <class Tp> struct istream_helper< Tp, decltype(std::declval<std::decay<decltype(std::declval<std::istream &>() >> std::declval<Tp &>())>>(), nullptr)> { istream_helper(std::istream &is, Tp &x) { is >> x; } }; #ifdef __SIZEOF_INT128__ template <> struct istream_helper<__uint128_t, std::nullptr_t> { istream_helper(std::istream &__is, __uint128_t &__x) { std::string __s; __is >> __s; bool __neg = false; if (__s.front() == '-') __neg = true, __s.erase(__s.begin()); __x = 0; for (char __d : __s) { __x *= 10; __d -= '0'; if (__neg) __x -= __d; else __x += __d; } } }; template <> struct istream_helper<__int128_t, std::nullptr_t> { istream_helper(std::istream &__is, __int128_t &__x) { std::string __s; __is >> __s; bool __neg = false; if (__s.front() == '-') __neg = true, __s.erase(__s.begin()); __x = 0; for (char __d : __s) { __x *= 10; __d -= '0'; if (__neg) __x -= __d; else __x += __d; } } }; #endif // INT128 template <class T1, class T2> struct istream_helper<std::pair<T1, T2>> { istream_helper(std::istream &is, std::pair<T1, T2> &x) { istream_helper<T1>(is, x.first), istream_helper<T2>(is, x.second); } }; template <class... Tps> struct istream_helper<std::tuple<Tps...>> { istream_helper(std::istream &is, std::tuple<Tps...> &x) { iterate(is, x); } private: template <class Tp, size_t N = 0> void iterate(std::istream &is, Tp &x) { if constexpr (N == std::tuple_size<Tp>::value) return; else istream_helper<typename std::tuple_element<N, Tp>::type>(is, std::get<N>(x)), iterate<Tp, N + 1>(is, x); } }; } // namespace internal /** * @brief A wrapper class for std::istream. */ class istream : public std::istream { public: /** * @brief Wrapped operator. */ template <typename Tp> istream &operator>>(Tp &x) { internal::istream_helper<Tp>(*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; } }; decltype(auto) cin = static_cast<istream &>(std::cin); } // namespace workspace #line 2 "Library/src/utils/io/ostream.hpp" /** * @file ostream.hpp * @brief Output Stream */ #line 9 "Library/src/utils/io/ostream.hpp" namespace workspace { template <class _Os> struct is_ostream { template <typename... _Args> static std::true_type __test(std::basic_ostream<_Args...> *); static std::false_type __test(void *); constexpr static bool value = decltype(__test(std::declval<_Os *>()))::value; }; template <class _Os> using ostream_ref = typename std::enable_if<is_ostream<_Os>::value, _Os &>::type; /** * @brief Stream insertion operator for C-style array. * * @param __os Output stream * @param __a Array * @return Reference to __os. */ template <class _Os, class _Tp, size_t _Nm> typename std::enable_if<bool(sizeof(_Tp) > 2), ostream_ref<_Os>>::type operator<<(_Os &__os, const _Tp (&__a)[_Nm]) { if constexpr (_Nm) { __os << *__a; for (auto __i = __a + 1, __e = __a + _Nm; __i != __e; ++__i) __os << ' ' << *__i; } return __os; } /** * @brief Stream insertion operator for std::pair. * * @param __os Output stream * @param __p Pair * @return Reference to __os. */ template <class _Os, class _T1, class _T2> ostream_ref<_Os> operator<<(_Os &__os, const std::pair<_T1, _T2> &__p) { return __os << __p.first << ' ' << __p.second; } /** * @brief Stream insertion operator for std::tuple. * * @param __os Output stream * @param __t Tuple * @return Reference to __os. */ template <class _Os, class _Tp, size_t _Nm = 0> typename std::enable_if<bool(std::tuple_size<_Tp>::value + 1), ostream_ref<_Os>>::type operator<<(_Os &__os, const _Tp &__t) { if constexpr (_Nm != std::tuple_size<_Tp>::value) { if constexpr (_Nm) __os << ' '; __os << std::get<_Nm>(__t); operator<<<_Os, _Tp, _Nm + 1>(__os, __t); } return __os; } template <class _Os, class _Container, typename = decltype(std::begin(std::declval<_Container>()))> typename std::enable_if< !std::is_same<typename std::decay<_Container>::type, std::string>::value && !std::is_same<typename std::decay<_Container>::type, char *>::value, ostream_ref<_Os>>::type operator<<(_Os &__os, const _Container &__cont) { bool __h = true; for (auto &&__e : __cont) __h ? __h = 0 : (__os << ' ', 0), __os << __e; return __os; } #ifdef __SIZEOF_INT128__ /** * @brief Stream insertion operator for __int128_t. * * @param __os Output Stream * @param __x 128-bit integer * @return Reference to __os. */ template <class _Os> ostream_ref<_Os> operator<<(_Os &__os, __int128_t __x) { if (!__x) return __os << '0'; if (__x < 0) __os << '-'; char __s[40], *__p = __s; while (__x) { auto __d = __x % 10; *__p++ = '0' + (__x < 0 ? -__d : __d); __x /= 10; } *__p = 0; for (char *__t = __s; __t < --__p; ++__t) *__t ^= *__p ^= *__t ^= *__p; return __os << __s; } /** * @brief Stream insertion operator for __uint128_t. * * @param __os Output Stream * @param __x 128-bit unsigned integer * @return Reference to __os. */ template <class _Os> ostream_ref<_Os> operator<<(_Os &__os, __uint128_t __x) { if (!__x) return __os << '0'; char __s[40], *__p = __s; while (__x) *__p++ = '0' + __x % 10, __x /= 10; *__p = 0; for (char *__t = __s; __t < --__p; ++__t) *__t ^= *__p ^= *__t ^= *__p; return __os << __s; } #endif } // namespace workspace #line 10 "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 { /** * @brief Setup I/O. * @param __n Standard output precision */ void io_setup(int __n) { std::cin.tie(0)->sync_with_stdio(0); std::cout << std::fixed << std::setprecision(__n); #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 <class Tuple, size_t N = std::tuple_size<Tuple>::value - 1> struct common_iterator_category { using type = typename std::common_type< typename common_iterator_category<Tuple, N - 1>::type, typename std::iterator_traits<typename std::tuple_element< N, Tuple>::type>::iterator_category>::type; }; template <class Tuple> struct common_iterator_category<Tuple, 0> { using type = typename std::iterator_traits< typename std::tuple_element<0, Tuple>::type>::iterator_category; }; } // namespace workspace #line 2 "Library/src/utils/iterator/reverse.hpp" /* * @file reverse_iterator.hpp * @brief Reverse Iterator */ #if __cplusplus >= 201703L #include <iterator> #include <optional> namespace workspace { /* * @class reverse_iterator * @brief Wrapper class for `std::reverse_iterator`. * @see http://gcc.gnu.org/PR51823 */ template <class Iterator> class reverse_iterator : public std::reverse_iterator<Iterator> { using base_std = std::reverse_iterator<Iterator>; std::optional<typename base_std::value_type> 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/make_vector.hpp" /* * @file make_vector.hpp * @brief Multi-dimensional Vector */ #if __cplusplus >= 201703L #include <tuple> #include <vector> 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 <typename Tp, size_t N, typename S> constexpr auto make_vector([[maybe_unused]] S* sizes, Tp const& init = Tp()) { static_assert(std::is_convertible_v<S, size_t>); if constexpr (N) return std::vector(*sizes, make_vector<Tp, N - 1, S>(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 <typename Tp, size_t N, typename S> constexpr auto make_vector(const S (&sizes)[N], Tp const& init = Tp()) { return make_vector<Tp, N, S>((S*)sizes, init); } /* * @brief Make a multi-dimensional vector. * @param sizes The size of each dimension * @param init The initial value */ template <typename Tp, size_t N, typename S, size_t I = 0> constexpr auto make_vector([[maybe_unused]] std::array<S, N> const& sizes, Tp const& init = Tp()) { static_assert(std::is_convertible_v<S, size_t>); if constexpr (I == N) return init; else return std::vector(sizes[I], make_vector<Tp, N, S, I + 1>(sizes, init)); } /* * @brief Make a multi-dimensional vector. * @param sizes The size of each dimension * @param init The initial value */ template <typename Tp, size_t N = SIZE_MAX, size_t I = 0, class... Args> constexpr auto make_vector([[maybe_unused]] std::tuple<Args...> const& sizes, Tp const& init = Tp()) { using tuple_type = std::tuple<Args...>; if constexpr (I == std::tuple_size_v<tuple_type> || I == N) return init; else { static_assert( std::is_convertible_v<std::tuple_element_t<I, tuple_type>, size_t>); return std::vector(std::get<I>(sizes), make_vector<Tp, N, I + 1>(sizes, init)); } } /* * @brief Make a multi-dimensional vector. * @param sizes The size of each dimension * @param init The initial value */ template <typename Tp, class Fst, class Snd> constexpr auto make_vector(std::pair<Fst, Snd> const& sizes, Tp const& init = Tp()) { static_assert(std::is_convertible_v<Fst, size_t>); static_assert(std::is_convertible_v<Snd, size_t>); return make_vector({(size_t)sizes.first, (size_t)sizes.second}, init); } } // namespace workspace #endif #line 2 "Library/src/utils/py-like/enumerate.hpp" /* * @file enumerate.hpp * @brief Enumerate */ #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/py-like/reversed.hpp" /** * @file reversed.hpp * @brief Reversed */ #include <initializer_list> #line 10 "Library/src/utils/py-like/reversed.hpp" namespace workspace { namespace internal { template <class Container> 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 <class Container> constexpr auto reversed(Container &&cont) noexcept { return internal::reversed<Container>{std::forward<Container>(cont)}; } template <class Tp> constexpr auto reversed(std::initializer_list<Tp> &&cont) noexcept { return internal::reversed<std::initializer_list<Tp>>{ std::forward<std::initializer_list<Tp>>(cont)}; } } // namespace workspace #line 12 "Library/src/utils/py-like/range.hpp" #if __cplusplus >= 201703L namespace workspace { template <class Index> 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<Index>::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<iterator> rbegin() const noexcept { return reverse_iterator<iterator>(end()); } constexpr reverse_iterator<iterator> rend() const noexcept { return reverse_iterator<iterator>(begin()); } }; template <class... Args> constexpr auto rrange(Args &&... args) noexcept { return internal::reversed(range(std::forward<Args>(args)...)); } } // namespace workspace #endif #line 2 "Library/src/utils/py-like/zip.hpp" /** * @file zip.hpp * @brief Zip */ #line 11 "Library/src/utils/py-like/zip.hpp" #line 14 "Library/src/utils/py-like/zip.hpp" #if __cplusplus >= 201703L namespace workspace { namespace internal { template <class> struct zipped_iterator; template <class...> struct zipped_iterator_tuple; template <class... Args> class zipped { using ref_tuple = std::tuple<Args...>; ref_tuple args; template <size_t N = 0> constexpr auto begin_cat() const noexcept { if constexpr (N != std::tuple_size<ref_tuple>::value) { return std::tuple_cat(std::tuple(std::begin(std::get<N>(args))), begin_cat<N + 1>()); } else return std::tuple<>(); } template <size_t N = 0> constexpr auto end_cat() const noexcept { if constexpr (N != std::tuple_size<ref_tuple>::value) { return std::tuple_cat(std::tuple(std::end(std::get<N>(args))), end_cat<N + 1>()); } else return std::tuple<>(); } public: constexpr zipped(Args &&... args) noexcept : args(args...) {} class iterator { using base_tuple = typename zipped_iterator_tuple<Args...>::type; public: using iterator_category = typename common_iterator_category<base_tuple>::type; using difference_type = std::ptrdiff_t; using value_type = zipped_iterator<base_tuple>; using reference = zipped_iterator<base_tuple> &; using pointer = iterator; protected: value_type current; template <size_t N = 0> constexpr bool equal(const iterator &rhs) const noexcept { if constexpr (N != std::tuple_size<base_tuple>::value) { return std::get<N>(current) == std::get<N>(rhs.current) || equal<N + 1>(rhs); } else return false; } template <size_t N = 0> constexpr void increment() noexcept { if constexpr (N != std::tuple_size<base_tuple>::value) { ++std::get<N>(current); increment<N + 1>(); } } template <size_t N = 0> constexpr void decrement() noexcept { if constexpr (N != std::tuple_size<base_tuple>::value) { --std::get<N>(current); decrement<N + 1>(); } } template <size_t N = 0> constexpr void advance(difference_type __d) noexcept { if constexpr (N != std::tuple_size<base_tuple>::value) { std::get<N>(current) += __d; advance<N + 1>(__d); } } public: constexpr iterator() noexcept = default; constexpr iterator(base_tuple const &current) noexcept : current(current) {} constexpr bool operator==(const iterator &rhs) const noexcept { return equal(rhs); } constexpr bool operator!=(const iterator &rhs) const noexcept { return !equal(rhs); } constexpr iterator &operator++() noexcept { increment(); return *this; } constexpr iterator &operator--() noexcept { decrement(); return *this; } constexpr bool operator<(const iterator &rhs) const noexcept { return std::get<0>(current) < std::get<0>(rhs.current); } constexpr bool operator<=(const iterator &rhs) const noexcept { return std::get<0>(current) <= std::get<0>(rhs.current); } constexpr iterator &operator+=(difference_type __d) noexcept { advance(__d); return *this; } constexpr iterator &operator-=(difference_type __d) noexcept { advance(-__d); return *this; } constexpr iterator operator+(difference_type __d) const noexcept { return iterator{*this} += __d; } constexpr iterator operator-(difference_type __d) const noexcept { return iterator{*this} -= __d; } constexpr difference_type operator-(const iterator &rhs) const noexcept { return std::get<0>(current) - std::get<0>(rhs.current); } constexpr reference operator*() noexcept { return current; } }; constexpr iterator begin() const noexcept { return iterator{begin_cat()}; } constexpr iterator end() const noexcept { return iterator{end_cat()}; } constexpr reverse_iterator<iterator> rbegin() const noexcept { return reverse_iterator<iterator>{end()}; } constexpr reverse_iterator<iterator> rend() const noexcept { return reverse_iterator<iterator>{begin()}; } }; template <class Tp, class... Args> struct zipped_iterator_tuple<Tp, Args...> { using type = decltype(std::tuple_cat( std::declval<std::tuple<decltype(std::begin(std::declval<Tp>()))>>(), std::declval<typename zipped_iterator_tuple<Args...>::type>())); }; template <> struct zipped_iterator_tuple<> { using type = std::tuple<>; }; template <class Iter_tuple> struct zipped_iterator : Iter_tuple { constexpr zipped_iterator(Iter_tuple const &__t) noexcept : Iter_tuple::tuple(__t) {} constexpr zipped_iterator(zipped_iterator const &__t) = default; constexpr zipped_iterator &operator=(zipped_iterator const &__t) = default; // Avoid move initialization. constexpr zipped_iterator(zipped_iterator &&__t) : zipped_iterator(static_cast<zipped_iterator const &>(__t)) {} // Avoid move assignment. zipped_iterator &operator=(zipped_iterator &&__t) { return operator=(static_cast<zipped_iterator const &>(__t)); } template <size_t N> friend constexpr auto &get(zipped_iterator<Iter_tuple> const &__z) noexcept { return *std::get<N>(__z); } template <size_t N> friend constexpr auto get(zipped_iterator<Iter_tuple> &&__z) noexcept { return *std::get<N>(__z); } }; } // namespace internal } // namespace workspace namespace std { template <size_t N, class Iter_tuple> struct tuple_element<N, workspace::internal::zipped_iterator<Iter_tuple>> { using type = typename remove_reference<typename iterator_traits< typename tuple_element<N, Iter_tuple>::type>::reference>::type; }; template <class Iter_tuple> struct tuple_size<workspace::internal::zipped_iterator<Iter_tuple>> : tuple_size<Iter_tuple> {}; } // namespace std namespace workspace { template <class... Args> constexpr auto zip(Args &&... args) noexcept { return internal::zipped<Args...>(std::forward<Args>(args)...); } template <class... Args> constexpr auto zip(std::initializer_list<Args> const &... args) noexcept { return internal::zipped<const std::initializer_list<Args>...>(args...); } } // namespace workspace #endif #line 10 "Library/src/utils/py-like/enumerate.hpp" #if __cplusplus >= 201703L namespace workspace { constexpr size_t min_size() noexcept { return SIZE_MAX; } template <class Container, class... Args> constexpr size_t min_size(Container const &cont, Args &&... args) noexcept { return std::min(std::size(cont), min_size(std::forward<Args>(args)...)); } template <class... Args> constexpr auto enumerate(Args &&... args) noexcept { return zip(range(min_size(args...)), std::forward<Args>(args)...); } template <class... Args> constexpr auto enumerate(std::initializer_list<Args> const &... args) noexcept { return zip(range(min_size(args...)), std::vector(args)...); } } // namespace workspace #endif #line 17 "Library/lib/utils" // #include "src/utils/py-like/reversed.hpp" // #include "src/utils/py-like/zip.hpp" #line 2 "Library/src/utils/rand/rng.hpp" /** * @file rng.hpp * @brief Random Number Generator */ #line 9 "Library/src/utils/rand/rng.hpp" namespace workspace { template <typename _Arithmetic> using uniform_distribution = typename std::conditional< std::is_integral<_Arithmetic>::value, std::uniform_int_distribution<_Arithmetic>, std::uniform_real_distribution<_Arithmetic>>::type; template <typename _Arithmetic, class _Engine = std::mt19937> class random_number_generator : uniform_distribution<_Arithmetic> { using base = uniform_distribution<_Arithmetic>; _Engine __engine; public: random_number_generator(_Arithmetic __min, _Arithmetic __max) : base(__min, __max), __engine(std::random_device{}()) {} random_number_generator(_Arithmetic __max = 1) : random_number_generator(0, __max) {} random_number_generator(typename base::param_type const& __param) : base(__param), __engine(std::random_device{}()) {} decltype(auto) operator()() noexcept { return base::operator()(__engine); } }; } // namespace workspace #line 2 "Library/src/utils/rand/shuffle.hpp" /** * @file shuffle.hpp * @brief Shuffle */ #line 10 "Library/src/utils/rand/shuffle.hpp" namespace workspace { template <class _RAIter, class _Engine = std::mt19937> void shuffle(_RAIter __first, _RAIter __last) { static _Engine __engine(std::random_device{}()); std::shuffle(__first, __last, __engine); } } // namespace workspace #line 2 "Library/src/utils/round_div.hpp" /* * @file round_div.hpp * @brief Round Integer Division */ #line 9 "Library/src/utils/round_div.hpp" #line 11 "Library/src/utils/round_div.hpp" namespace workspace { /* * @fn floor_div * @brief floor of fraction. * @param x the numerator * @param y the denominator * @return maximum integer z s.t. z <= x / y * @note y must be nonzero. */ template <typename T1, typename T2> constexpr typename std::enable_if<(is_integral_ext<T1>::value && is_integral_ext<T2>::value), typename std::common_type<T1, T2>::type>::type floor_div(T1 x, T2 y) { assert(y != 0); if (y < 0) x = -x, y = -y; return x < 0 ? (x - y + 1) / y : x / y; } /* * @fn ceil_div * @brief ceil of fraction. * @param x the numerator * @param y the denominator * @return minimum integer z s.t. z >= x / y * @note y must be nonzero. */ template <typename T1, typename T2> constexpr typename std::enable_if<(is_integral_ext<T1>::value && is_integral_ext<T2>::value), typename std::common_type<T1, T2>::type>::type ceil_div(T1 x, T2 y) { assert(y != 0); if (y < 0) x = -x, y = -y; return x < 0 ? x / y : (x + y - 1) / y; } } // namespace workspace #line 22 "Library/lib/utils" // #include "src/utils/sfinae.hpp" #line 11 "atcoder-workspace/2.cc" signed main() { using namespace workspace; io_setup(15); /* given case_info.read(); //*/ /* unspecified case_info.total = -1; //*/ return case_info.iterate(); } #line 2 "Library/src/graph/directed/flow/Dinic.hpp" /** * @file Dinic.hpp * @brief Dinic's Algorithm * @date 2021-01-15 * * */ #line 12 "Library/src/graph/directed/flow/Dinic.hpp" #line 2 "Library/src/graph/directed/flow/base.hpp" /** * @file base.hpp * @brief Flow Graph * @date 2021-01-15 * * */ #line 15 "Library/src/graph/directed/flow/base.hpp" namespace workspace { template <class Cap, class Cost = void> class flow_graph { protected: class adjacency_impl; public: using container_type = std::vector<adjacency_impl>; using size_type = typename container_type::size_type; class unweighted_edge { public: size_type src; // Source size_type dst; // Destination Cap cap; // Capacity Cap flow = 0; // Flow unweighted_edge(size_type __s, size_type __d, const Cap &__u = 1) : src(__s), dst(__d), cap(__u) { assert(!(cap < static_cast<Cap>(0))); } /** * @brief Source, Destination, Capacity, Flow */ template <class Os> friend Os &operator<<(Os &__os, unweighted_edge const &__e) { return __os << __e.src << ' ' << __e.dst << ' ' << __e.cap << ' ' << __e.flow; } protected: unweighted_edge() = default; unweighted_edge(size_type __s, size_type __d, const Cap &__u, const Cap &__f) : src(__s), dst(__d), cap(__u), flow(__f) {} unweighted_edge make_rev() const { return {dst, src, flow, cap}; } }; class weighted_edge : public unweighted_edge { public: Cost cost; // Cost weighted_edge(const unweighted_edge &__e, const Cost &__c = 0) : unweighted_edge(__e), cost(__c) {} weighted_edge(size_type __s, size_type __d, const Cap &__u = 1, const Cost &__c = 0) : unweighted_edge(__s, __d, __u), cost(__c) {} /** * @brief Source, Destination, Capacity, Flow, Cost */ template <class Os> friend Os &operator<<(Os &__os, weighted_edge const &__e) { return __os << static_cast<unweighted_edge>(__e) << ' ' << __e.cost; } protected: weighted_edge() = default; weighted_edge make_rev() const { return {unweighted_edge::make_rev(), -cost}; } }; using edge = typename std::conditional<std::is_void<Cost>::value, unweighted_edge, weighted_edge>::type; protected: struct edge_impl : edge { bool aux = false; edge_impl *rev = nullptr; edge_impl() = default; edge_impl(const edge_impl &__e) = default; edge_impl &operator=(const edge_impl &__e) = default; edge_impl(edge_impl &&__e) = default; edge_impl &operator=(edge_impl &&__e) = default; edge_impl(const edge &__e) : edge(__e) {} edge_impl(edge &&__e) : edge(__e) {} void push(Cap __f) { edge::cap -= __f; edge::flow += __f; if (rev) { rev->cap += __f; rev->flow -= __f; } } edge_impl make_rev() { edge_impl __e = edge::make_rev(); __e.aux = true; __e.rev = this; return __e; } }; public: class adjacency { public: using value_type = edge; using reference = edge &; using const_reference = edge const &; using pointer = edge *; using const_pointer = const edge *; class iterator { edge_impl *__p; public: iterator(edge_impl *__p = nullptr) : __p(__p) {} bool operator!=(iterator const &__x) const { return __p != __x.__p; } iterator &operator++() { do ++__p; while (__p->rev && __p->aux); return *this; } pointer operator->() const { return __p; } reference operator*() const { return *__p; } }; class const_iterator { const edge_impl *__p; public: const_iterator(const edge_impl *__p = nullptr) : __p(__p) {} bool operator!=(const_iterator const &__x) const { return __p != __x.__p; } const_iterator &operator++() { do ++__p; while (__p->rev && __p->aux); return *this; } const_pointer operator->() const { return __p; } const_reference operator*() const { return *__p; } }; adjacency() : first(new edge_impl[2]), last(first + 1), __s(first), __t(first) {} ~adjacency() { delete[] first; } const_reference operator[](size_type i) const { assert(i < size()); return *(first + i); } size_type size() const { return __t - first; } auto begin() { return iterator{__s}; } auto begin() const { return const_iterator{__s}; } auto end() { return iterator{__t}; } auto end() const { return const_iterator{__t}; } /** * @brief Construct a new adjacency object * * @param __x Rvalue reference to another object */ adjacency(adjacency &&__x) : first(nullptr) { operator=(std::move(__x)); } /** * @brief Assignment operator. * * @param __x Rvalue reference to another object * @return Reference to this object. */ adjacency &operator=(adjacency &&__x) { std::swap(first, __x.first); last = __x.last; __s = __x.__s; __t = __x.__t; return *this; } protected: edge_impl *first, *last, *__s, *__t; }; using value_type = adjacency; using reference = adjacency &; using const_reference = adjacency const &; protected: class adjacency_impl : public adjacency { public: using base = adjacency; using base::__s; using base::__t; using base::first; using base::last; using iterator = edge_impl *; iterator _push(edge_impl &&__e) { if (__t == last) { size_type __n(last - first); iterator loc = new edge_impl[__n << 1 | 1]; __s += loc - first; __t = loc; for (iterator __p{first}; __p != last; ++__p, ++__t) { *__t = *__p; if (__p->rev) __p->rev->rev = __t; } delete[] first; first = loc; last = __t + __n; } *__t = std::move(__e); if (__s->aux) ++__s; return __t++; } iterator begin() const { return first; } iterator end() const { return __t; } }; /** * @brief The only member variable. */ container_type graph; public: /** * @brief Construct a new flow graph object * * @param __n Number of vertices */ flow_graph(size_type __n = 0) : graph(__n) {} /** * @brief Construct a new flow graph object * * @param __x Const reference to another object */ flow_graph(const flow_graph &__x) : graph(__x.size()) { for (auto &&__adj : __x) for (auto &&__e : __adj) add_edge(__e); } /** * @brief Construct a new flow graph object * * @param __x Rvalue reference to another object */ flow_graph(flow_graph &&__x) : graph(std::move(__x.graph)) {} /** * @brief Assignment operator. * * @param __x Const reference to another object * @return Reference to this object. */ flow_graph &operator=(const flow_graph &__x) { return operator=(std::move(flow_graph{__x})); } /** * @brief Assignment operator. * * @param __x Rvalue reference to another object * @return Reference to this object. */ flow_graph &operator=(flow_graph &&__x) { graph = std::move(__x.graph); return *this; } /** * @return Whether the graph is empty. */ bool empty() const { return graph.empty(); } /** * @return Number of nodes. */ size_type size() const { return graph.size(); } /** * @param node Node * @return Referece to the adjacency list of the node. */ reference operator[](size_type node) { assert(node < size()); return graph[node]; } /** * @param node Node * @return Const referece to the adjacency list of the node. */ const_reference operator[](size_type node) const { assert(node < size()); return graph[node]; } class iterator : public container_type::iterator { using base = typename container_type::iterator; public: using reference = adjacency &; using pointer = adjacency *; iterator(base const &__i) : base(__i) {} pointer operator->() const { return base::operator->(); } reference operator*() const { return base::operator*(); } }; class const_iterator : public container_type::const_iterator { using base = typename container_type::const_iterator; public: using const_reference = const adjacency &; using const_pointer = const adjacency *; const_iterator(base const &__i) : base(__i) {} const_pointer operator->() const { return base::operator->(); } const_reference operator*() const { return base::operator*(); } }; auto begin() { return iterator{graph.begin()}; } auto begin() const { return const_iterator{graph.begin()}; } auto end() { return iterator{graph.end()}; } auto end() const { return const_iterator{graph.end()}; } /** * @brief Add a node to the graph. * * @return Index of the node. */ size_type add_node() { return add_nodes(1).front(); } /** * @brief Add some nodes to the graph. * * @param __n Number of nodes added * @return List of indices of the nodes. */ virtual std::vector<size_type> add_nodes(size_type __n) { std::vector<size_type> __nds(__n); std::iota(__nds.begin(), __nds.end(), graph.size()); __n += graph.size(); if (__n > graph.capacity()) { flow_graph __x(__n); for (auto &&adj : graph) for (auto &&__e : adj) if (!__e.aux) __x.add_edge(__e); graph = std::move(__x.graph); } else graph.resize(__n); return __nds; } /** * @brief Add a directed edge to the graph. * * @return Reference to the edge. */ template <class... Args> typename std::enable_if<std::is_constructible<edge, Args...>::value, edge &>::type add_edge(Args &&... __args) { edge_impl __e = edge(std::forward<Args>(__args)...); assert(__e.src < size()); assert(__e.dst < size()); edge_impl *__p = graph[__e.src]._push(std::move(__e)); // Careful with a self loop. if (__e.src != __e.dst) __p->rev = graph[__e.dst]._push(__p->make_rev()); return *__p; } /** * @brief Add a directed edge to the graph. * * @return Reference to the edge. */ template <class Tp> typename std::enable_if< (std::tuple_size<typename std::decay<Tp>::type>::value >= 0), edge &>::type add_edge(Tp &&__t) { return _unpack_directed(std::forward<Tp>(__t)); } /** * @brief Add an undirected edge to the graph. Its cost must be non-negative. * * @return Reference to the edge. */ template <class... Args> edge &add_undirected_edge(Args &&... __args) { edge_impl __e = edge(std::forward<Args>(__args)...); assert(__e.src < size()); assert(__e.dst < size()); (__e.flow += __e.flow) += __e.cap; edge_impl *__p = graph[__e.src]._push(std::move(__e)); // Careful with a self loop. if (__e.src != __e.dst) { edge_impl __r = __p->make_rev(); __r.aux = false; __p->rev = graph[__e.dst]._push(std::move(__r)); } return *__p; } /** * @brief Add an undirected edge to the graph. Its cost must be non-negative. * * @return Reference to the edge. */ template <class Tp> typename std::enable_if< (std::tuple_size<typename std::decay<Tp>::type>::value >= 0), edge &>::type add_undirected_edge(Tp &&__t) { return _unpack_undirected(std::forward<Tp>(__t)); } protected: // internal template <class Tp, size_t N = 0, class... Args> decltype(auto) _unpack_directed(Tp &&__t, Args &&... __args) { if constexpr (N == std::tuple_size<typename std::decay<Tp>::type>::value) return add_edge(std::forward<Args>(__args)...); else return _unpack_directed<Tp, N + 1>(std::forward<Tp>(__t), std::forward<Args>(__args)..., std::get<N>(__t)); } // internal template <class Tp, size_t N = 0, class... Args> decltype(auto) _unpack_undirected(Tp &&__t, Args &&... __args) { if constexpr (N == std::tuple_size<typename std::decay<Tp>::type>::value) return add_undirected_edge(std::forward<Args>(__args)...); else return _unpack_undirected<Tp, N + 1>(std::forward<Tp>(__t), std::forward<Args>(__args)..., std::get<N>(__t)); } template <class Os> friend Os &operator<<(Os &__os, flow_graph const &__g) { for (const auto &adj : __g) for (const auto &e : adj) __os << e << "\n"; return __os; } }; } // namespace workspace #line 14 "Library/src/graph/directed/flow/Dinic.hpp" namespace workspace { /** * @brief Compute the maximum flow. * * @tparam Cap Capacity type */ template <class Cap> class Dinic : public flow_graph<Cap> { using base = flow_graph<Cap>; public: using size_type = typename base::size_type; protected: constexpr static size_type nil = -1; std::vector<size_type> level; std::vector<typename base::container_type::value_type::iterator> iter; Cap dfs(size_type __s, size_type __d, Cap bound) { if (__s == __d) return bound; Cap flow(0); for (auto &e{iter[__d]}; e != base::graph[__d].end(); ++e) if (static_cast<Cap>(0) < e->flow && level[e->dst] < level[__d]) if (Cap achv = dfs(__s, e->dst, std::min(bound, e->flow)); static_cast<Cap>(0) < achv) { e->push(-achv); flow += achv, bound -= achv; if (bound == static_cast<Cap>(0)) break; } return flow; } public: /** * @brief Construct a new Dinic object * * @param __n Number of nodes */ Dinic(size_type __n = 0) : base::flow_graph(__n), level(__n, nil), iter(__n) {} /** * @brief Add some nodes to the graph. * * @param __n Number of nodes added * @return List of indices of the nodes. */ std::vector<size_type> add_nodes(size_type __n) override { auto __nds = base::add_nodes(__n); level.resize(base::size(), nil); iter.resize(base::size()); return __nds; } /** * @brief Run Dinic's algorithm. * * @param __s Source * @param __d Destination * @return Maximum flow. */ Cap run(size_type __s, size_type __d) { assert(__s < base::size()); assert(__d < base::size()); assert(__s != __d); Cap flow = 0, bound = std::numeric_limits<Cap>::max(); for (std::vector<size_type> que(base::size());; std::fill(level.begin(), level.end(), nil)) { level[que.front() = __s] = 0; for (auto ql{que.begin()}, qr{std::next(ql)}; level[__d] == nil && ql != qr; ++ql) for (const auto &e : base::graph[*ql]) if (static_cast<Cap>(0) < e.cap && level[e.dst] == nil) level[ *qr++ = e.dst] = level[*ql] + 1; if (level[__d] == nil) break; for (size_type node{}; node != base::size(); ++node) iter[node] = base::graph[node].begin(); flow += dfs(__s, __d, bound); } return flow; } }; } // namespace workspace #line 2 "Library/src/graph/directed/flow/min_cost_flow.hpp" /** * @file min_cost_flow.hpp * @brief Minimum Cost Flow * @date 2021-01-15 * * */ #line 13 "Library/src/graph/directed/flow/min_cost_flow.hpp" #line 16 "Library/src/graph/directed/flow/min_cost_flow.hpp" namespace workspace { /** * @brief Capacity Scaling Algorithm. * * @tparam Cap Capacity type * @tparam Cost Cost type */ template <class Cap, class Cost = Cap> class min_cost_flow : public flow_graph<Cap, Cost> { using base = flow_graph<Cap, Cost>; using edge_impl = typename base::edge_impl; public: using edge = typename base::edge; using size_type = typename base::size_type; /** * @brief Construct a new min_cost_flow object * * @param __n Number of vertices */ min_cost_flow(size_type __n = 0) : base::flow_graph(__n), b(__n) {} std::vector<size_type> add_nodes(size_type __n) override { b.resize(b.size() + __n); return base::add_nodes(__n); } using base::add_edge; /** * @brief Add a directed edge to the graph. * * @param __s Source * @param __d Destination * @param __l Lower bound of flow * @param __u Upper bound of flow * @param __c Cost * @return Reference to the edge. */ edge &add_edge(size_type __s, size_type __d, Cap __l, Cap __u, Cost __c) { assert(!(__u < __l)); b[__s] -= __l; b[__d] += __l; auto &__e = base::add_edge(__s, __d, __u - __l, __c); __e.flow = __l; return __e; } /** * @brief Add an undirected edge to the graph. * * @return Reference to the edge. */ template <class... Args> edge &add_undirected_edge(Args &&... __args) { auto &__e = static_cast<edge_impl &>( base::add_undirected_edge(std::forward<Args>(__args)...)); assert(!(__e.cost < 0)); __e.rev->cost = __e.cost; return __e; } /** * @brief Increase the balance of a node. * * @param node * @param __f Default: 1 */ void supply(size_type node, Cap __f = 1) { assert(node < b.size()); b[node] += __f; } /** * @brief Decrease the balance of a node. * * @param node * @param __f Default: 1 */ void demand(size_type node, Cap __f = 1) { assert(node < b.size()); b[node] -= __f; } /** * @return Balance. */ const auto &balance() const { return b; } /** * @param node Node * @return Balance of the node. */ Cap balance(size_type node) const { return b[node]; } /** * @return Potential. The dual solution. */ const auto &potential() const { return p; } /** * @param node Node * @return Potential of the node. */ Cost potential(size_type node) const { return p[node]; } /** * @return Cost of current flow. */ Cost cost() const { return current; } /** * @brief Run Capacity Scaling Algorithm. * * @return Whether a balanced flow exists. */ bool run() { p.resize(b.size()); Cap delta = 0; for (auto &&__adj : base::graph) for (auto &&__e : __adj) delta = std::max(delta, __e.cap); if (delta == static_cast<Cap>(0)) return std::all_of(b.begin(), b.end(), [](Cap __x) { return __x == static_cast<Cap>(0); }); parent.resize(b.size()); while (static_cast<Cap>(0) < delta) { delta /= 2; for (auto &&__adj : base::graph) for (auto &&__e : __adj) if (delta < __e.cap && __e.cost < p[__e.dst] - p[__e.src]) { b[__e.src] -= __e.cap; b[__e.dst] += __e.cap; __e.push(__e.cap); } sources.clear(); sinks.clear(); for (size_type __v = 0; __v != b.size(); ++__v) if (delta < b[__v]) sources.emplace_back(__v); else if (b[__v] < -delta) sinks.emplace_back(__v); while (dual(delta)) { primal(delta); sources.erase( std::remove_if(sources.begin(), sources.end(), [&](auto __v) { return !(delta < b[__v]); }), sources.end()); sinks.erase( std::remove_if(sinks.begin(), sinks.end(), [&](auto __v) { return !(b[__v] < -delta); }), sinks.end()); } } current = 0; for (auto &&__adj : base::graph) for (auto &&__e : __adj) if (!__e.aux) current += __e.cost * __e.flow; return sources.empty() && sinks.empty(); } protected: // Cost of flow. Cost current{}; // Balance std::vector<Cap> b; // The dual solution. std::vector<Cost> p; std::vector<edge_impl *> parent; std::vector<size_type> sources, sinks; // Augment along the dual solution. void primal(Cap delta) { for (auto __t : sinks) if (parent[__t]) { auto __f = -b[__t]; auto __s = __t; while (parent[__s]) __f = std::min(__f, parent[__s]->cap), __s = parent[__s]->src; if (delta < b[__s]) { __f = std::min(__f, b[__s]); b[__s] -= __f; b[__t] += __f; for (auto *__p = parent[__t]; __p; __p = parent[__p->src]) { __p->push(__f); parent[__p->dst] = nullptr; } } } } // Improve the dual solution. bool dual(Cap delta) { std::fill(parent.begin(), parent.end(), nullptr); size_type reachable = 0; struct state { size_type __v; Cost __d; state(size_type __v, Cost __d) : __v(__v), __d(__d) {} bool operator<(const state &__x) const { return __x.__d < __d; } }; std::priority_queue<state> __q; decltype(p) __nx(p.size(), numeric_limits<Cost>::max()); Cost __ld = 0; for (auto __v : sources) { __nx[__v] = p[__v]; __q.emplace(__v, 0); } while (!__q.empty()) { auto [__v, __d] = __q.top(); __q.pop(); if (__d + p[__v] != __nx[__v]) continue; __ld = __d; if (b[__v] < -delta && ++reachable == sinks.size()) break; for (auto &__e : base::graph[__v]) if (delta < __e.cap && (__d = __nx[__v] + __e.cost) < __nx[__e.dst]) { __q.emplace(__e.dst, (__nx[__e.dst] = __d) - p[__e.dst]); parent[__e.dst] = &__e; } } for (size_type __v = 0; __v != p.size(); ++__v) p[__v] = std::min(__nx[__v], p[__v] += __ld); return reachable; } }; // namespace workspace /** * @brief Capacity Scaling Algorithm. * * @tparam Cap Capacity type * @tparam Gain Gain type */ template <class Cap, class Gain = Cap> class max_gain_flow : public min_cost_flow<Cap, Gain> { using base = min_cost_flow<Cap, Gain>; using base::cost; public: using base::min_cost_flow; using edge = typename base::edge; /** * @brief Add a directed edge to the graph. The default capacity is 1. * * @return Reference to the edge. */ template <class... Args> decltype(auto) add_edge(Args &&... __args) { return add_edge(std::tuple<Args...>{std::forward<Args>(__args)...}); } /** * @brief Add a directed edge to the graph. The default capacity is 1. * * @return Reference to the edge. */ template <class Tp> typename std::enable_if< (std::tuple_size<typename std::decay<Tp>::type>::value >= 0), const edge &>::type add_edge(Tp __t) { std::get<std::tuple_size<decltype(__t)>::value - 1>(__t) *= -1; // Flip the sign of cost. return base::add_edge(std::move(__t)); } /** * @return Gain of current flow. */ Gain gain() const { return -cost(); } }; } // namespace workspace #line 2 "Library/src/graph/directed/strongly_connected_components.hpp" /* * @file strongly_connected_components.hpp * @brief Strongly Connected Components */ #line 11 "Library/src/graph/directed/strongly_connected_components.hpp" namespace workspace { struct strongly_connected_components { strongly_connected_components(size_t n) : graph(n) {} /* * @fn add_edge * @brief Add an edge. * @param src Source * @param dst Destination */ void add_edge(size_t src, size_t dst) { assert(src < size()); assert(dst < size()); graph[src].emplace_back(dst); } /* * @fn count * @return The number of components. */ size_t count() { assert(made()); return dag.size(); } /* * @fn size * @return The number of vertices. */ size_t size() const { return graph.size(); } /* * @fn operator[] * @param v Vertex * @return The component including given vertex. */ size_t operator[](size_t v) { assert(made()); return low[v]; } /* * @fn shrinked_dag * @return Directed Acyclic Graph consisting of components. */ const std::vector<std::vector<size_t>> &shrinked_dag() { assert(made()); return dag; } /* * @fn make * @brief Run SCC decomposition in linear time. */ void make() { low.assign(size(), 0); dag.clear(); size_t *ptr = new size_t[size()]; bool *const used = new bool[size()]; for (size_t v{}, c{}; v != size(); ++v) affix(v, c, ptr, used + size()); delete[] ptr; delete[] used; for (auto &e : low) e += dag.size(); reverse(begin(dag), end(dag)); for (auto &arcs : dag) for (auto &dst : arcs) dst += dag.size(); } protected: std::vector<std::vector<size_t>> graph, dag; std::vector<size_t> low; bool made() const { return !low.empty(); } /* * @fn affix * @param src Vertex * @param c Counter * @param ptr Pointer to a stack * @param used Negative indexed * @return Low-link number of the vertex. */ size_t affix(size_t src, size_t &c, size_t *&ptr, bool *const used) { if (low[src]) return low[src]; const size_t idx = ++c; low[src] = idx; *ptr++ = src; for (const size_t dst : graph[src]) low[src] = std::min(low[src], affix(dst, c, ptr, used)); if (low[src] == idx) { dag.push_back({}); const size_t ccnt = dag.size(); used[-ccnt] = true; size_t *srcp = ptr; do low[*--srcp] = -ccnt; while (*srcp != src); while (ptr != srcp) { const size_t now = *--ptr; for (const size_t dst : graph[now]) { if (!used[(int)low[dst]]) { dag.back().emplace_back(low[dst]); used[(int)low[dst]] = true; } } } for (const size_t dst : dag.back()) used[dst] = false; used[-ccnt] = false; return idx; } return low[src]; } }; } // namespace workspace #line 3 "Library/src/graph/undirected/tree/centroid_decomposition.hpp" class centroid_decomposition { std::vector<std::vector<size_t>> tree, sub; std::vector<size_t> siz; std::vector<bool> chkd; bool dfs(size_t nod, size_t pre, size_t all, size_t rot) { siz[nod] = 1; bool cent = true; for(size_t to : tree[nod]) { if(to == pre || chkd[to]) continue; if(dfs(to, nod, all, rot)) return true; siz[nod] += siz[to]; if(siz[to] > all / 2) cent = false; } if(all - siz[nod] > all / 2) cent = false; if(!cent) return false; chkd[nod] = true; if(~rot) sub[rot].emplace_back(nod); else root = nod; for(size_t to : tree[nod]) { if(chkd[to]) continue; if(to == pre) assert(dfs(to, nod, all - siz[nod], nod)); else assert(dfs(to, nod, siz[to], nod)); } return true; } public: size_t root; centroid_decomposition(size_t _n = 0) : tree(_n), sub(_n), siz(_n), chkd(_n) {} size_t size() const { return tree.size(); } void add_edge(size_t u, size_t v) { assert(u < size()); assert(v < size()); tree[u].emplace_back(v); tree[v].emplace_back(u); } size_t make() { assert(dfs(0, 0, size(), -1)); return root; } decltype(sub)::const_reference operator[](size_t nod) const { return sub[nod]; } }; // class centroid_decomposition #line 2 "Library/src/graph/undirected/tree/diameter.hpp" /* * @file diameter.hpp * @brief Diameter */ #line 9 "Library/src/graph/undirected/tree/diameter.hpp" namespace workspace { /* * @fn diameter * @brief Calculate a diameter of given tree in O(|V|) time. * @param tree Adjacency list of the tree * @return A longest shortest path between any pair of vertices. */ std::vector<size_t> diameter(std::vector<std::vector<size_t>> const& tree) { std::vector<size_t> diam; if (!tree.empty()) { } return diam; } } // namespace workspace #line 2 "Library/src/graph/undirected/tree/heavy_light_decomposition.hpp" /* * @file heavy_light_decomposition.hpp * @brief Heavy-Light Decomposition */ #line 11 "Library/src/graph/undirected/tree/heavy_light_decomposition.hpp" class heavy_light_decomposition { std::vector<std::vector<size_t>> tree; std::vector<size_t> sorted, in, out, head, __dpt; size_t sort_children(size_t node, size_t prev) { size_t sum = 0, max_size = 0; for (size_t &to : tree[node]) { if (to == prev) continue; size_t child_size = sort_children(to, node); sum += child_size; if (max_size < child_size) { max_size = child_size; std::swap(tree[node].front(), to); } } return sum + 1; } void traverse(size_t node, size_t prev) { in[node] = sorted.size(); sorted.emplace_back(node); if (tree[node].front() != prev) { for (const size_t to : tree[node]) if (to != prev) head[to] = node + size(), __dpt[to] = __dpt[node] + 1; head[tree[node].front()] = head[node] < size() ? head[node] : node; for (const size_t to : tree[node]) if (to != prev) traverse(to, node); } out[node] = sorted.size(); } bool made() const { return !sorted.empty(); } public: using interval = std::pair<size_t, size_t>; heavy_light_decomposition(size_t n = 0) : tree(n), in(n), out(n), head(n), __dpt(n) {} /** * @return The size of the tree. */ size_t size() const { return tree.size(); } /** * @param node The root of the subtree * @return The size of the subtree. */ size_t size(size_t node) const { assert(made()); return out[node] - in[node]; } void add_edge(size_t u, size_t v) { assert(u < size()); assert(v < size()); tree[u].emplace_back(v); tree[v].emplace_back(u); } const decltype(tree) &graph() const { return tree; } /** * @brief Run HLD with given root in linear time. * @param root The root node. */ void make(size_t root) { sorted.clear(); sort_children(root, root); __dpt[root] = 0; traverse(root, root); } size_t prev_node(size_t node) const { assert(made()); assert(in[node]); return sorted[in[node] - 1]; } size_t next_node(size_t node) const { assert(made()); assert(in[node] + 1 < size()); return sorted[in[node] + 1]; } size_t index(size_t node) const { assert(made()); return in[node]; } size_t node(size_t index) const { assert(made()); return sorted[index]; } /** * @return The current root of the tree. */ size_t root() const { assert(made()); return sorted.front(); } /** * @param root The root of the subtree. * @return The interval representing the subtree. */ interval subtree(size_t root) const { assert(made()); return {in[root], out[root]}; } size_t parent(size_t node) const { assert(made()); return head[node] < size() ? prev_node(node) : head[node] - size(); } size_t top(size_t node) const { assert(made()); return head[node] < size() ? head[node] : node; } /** * @brief Get LCA in O(log(size)) time. * @param u 1st node * @param v 2nd node * @return Lowest Common Ancestor of the two. */ size_t lca(size_t u, size_t v) const { assert(made()); if (in[v] < in[u]) std::swap(u, v); if (in[v] < out[u]) return u; while (in[u] < in[v]) v = parent(top(v)); return v; } size_t ancestor(size_t __v, size_t __k) const { assert(made()); while (__k) { assert(in[__v]); auto __t = top(__v); auto __d = in[__v] - in[__t]; if (__d < __k) { __k -= __d + 1; __v = head[__t] - size(); } else { __v = sorted[in[__v] - __k]; __k = 0; } } return __v; } size_t depth(size_t __v) const { return __dpt[__v]; } // O(log(n)) std::pair<std::vector<interval>, std::vector<interval>> path_decomposition( size_t u, size_t v) const { assert(made()); if (in[v] < in[u]) std::swap(u, v); std::vector<interval> left, right; size_t utop = top(u), vtop = top(v); while (utop != vtop) { left.emplace_back(in[vtop], in[v] + 1); vtop = top(v = parent(vtop)); if (in[v] < in[u]) { std::swap(u, v); std::swap(utop, vtop); std::swap(left, right); } } left.emplace_back(in[u], in[v] + 1); return std::make_pair(left, right); } // O(log(n)) std::vector<interval> path_decomposition(size_t node) const { assert(made()); auto [left, right] = path_decomposition(root(), node); assert(left.size() == 1); right.insert(right.begin(), left.front()); return right; } }; #line 2 "Library/src/graph/undirected/tree/lowest_common_ancestor.hpp" /** * @file lowest_common_ancestor.hpp * @brief Lowest Common Ancestor */ #line 10 "Library/src/graph/undirected/tree/lowest_common_ancestor.hpp" template <class Weight = void> class lowest_common_ancestor { struct unweighted_edge { size_t dest; }; struct weighted_edge : unweighted_edge { Weight weight; }; public: constexpr static bool is_weighted = !std::is_void<Weight>::value; using edge = typename std::conditional<is_weighted, weighted_edge, unweighted_edge>::type; std::vector<std::vector<edge>> __tree; std::vector<std::vector<size_t>> table; std::vector<size_t> index; void traverse(size_t node, size_t prev) { index[node] = table.front().size(); table.front().emplace_back(node); for (edge const &e : __tree[node]) { if (e.dest != prev) { traverse(e.dest, node); table.front().emplace_back(node); } } } void make_table() { for (size_t w = 1; w < size(); w <<= 1) { auto &curr(table.back()), next(curr); for (size_t i = 0, j = w; j != curr.size(); ++i, ++j) if (index[curr[j]] < index[curr[i]]) next[i] = curr[j]; table.emplace_back(next); } } public: lowest_common_ancestor(size_t n = 0) : __tree(n), table(1), index(n) {} /** * @return The number of nodes. */ size_t size() const { return __tree.size(); } const auto &tree() const { return __tree; } template <class _Tp> void add_edge(size_t __u, size_t __v, _Tp const &__w) { static_assert(is_weighted); assert(__u < size()); assert(__v < size()); __tree[__u].push_back({__v, __w}); __tree[__v].push_back({__u, __w}); } void add_edge(size_t __u, size_t __v) { static_assert(!is_weighted); assert(__u < size()); assert(__v < size()); __tree[__u].push_back({__v}); __tree[__v].push_back({__u}); } // O(n log(n)) void make(size_t root) { assert(root < size()); traverse(root, root); make_table(); } // O(1) size_t get(size_t u, size_t v) const { assert(u < size()); assert(v < size()); if (u == v) return u; if ((v = index[v]) < (u = index[u])) std::swap(u, v); size_t h = std::__lg(v - u); return std::min(table[h][u], table[h][v - (1 << h)]); } Weight fold(size_t __u, size_t __v) const { assert(__u < size()); assert(__v < size()); Weight __w; // if(__u != __v) {} return __w; } }; #line 4 "Library/src/graph/undirected/two_edge_connected_components.hpp" // instance: an undirected and not necessarily simple graph class two_edge_connected_component { static constexpr size_t nil = -1; std::vector<size_t> stack, low, comp; std::vector<std::vector<size_t>> graph, tree, memb; void make(size_t now, size_t pre) { size_t ord = low[now] = stack.size(); stack.emplace_back(now); std::vector<size_t> brid; for (size_t to : graph[now]) { if (to == pre) { pre = nil; continue; } if (low[to] == nil) make(to, now); if (low[to] > ord) { brid.emplace_back(to); graph[to].emplace_back(now); } else if (low[now] > low[to]) low[now] = low[to]; } brid.swap(graph[now]); if (ord == low[now]) { auto pos = stack.end(); tree.resize(tree.size() + 1); auto &adjc = tree.back(); do { --pos; comp[*pos] = memb.size(); for (size_t u : graph[*pos]) adjc.emplace_back(comp[u]); } while (*pos != now); memb.emplace_back(pos, stack.end()); stack.erase(pos, stack.end()); } } public: two_edge_connected_component(size_t n) : comp(n), graph(n) { stack.reserve(n), tree.reserve(n), memb.reserve(n); } void add_edge(size_t u, size_t v) { assert(u < size()), assert(v < size()); graph[u].emplace_back(v), graph[v].emplace_back(u); } void make() { low.assign(size(), nil); for (size_t v = 0; v != size(); ++v) if (low[v] == nil) make(v, nil); } size_t size() const { return graph.size(); } size_t size(size_t i) const { assert(i < count()); return memb[i].size(); } size_t count() const { return memb.size(); } size_t operator[](size_t v) const { assert(v < size()); return comp[v]; } const std::vector<size_t> &bridge(size_t v) const { assert(v < size()); return graph[v]; } const std::vector<size_t> &component(size_t i) const { assert(i < count()); return memb[i]; } const std::vector<std::vector<size_t>> &bridge_tree() const { return tree; } }; // class two_edge_connected_component #line 2 "Library/src/modular/inverse.hpp" /* * @file inverse.hpp * @brief Inverse Table */ #line 9 "Library/src/modular/inverse.hpp" #line 2 "Library/src/modular/modint.hpp" /** * @file modint.hpp * * @brief Modular Arithmetic */ #line 12 "Library/src/modular/modint.hpp" #line 14 "Library/src/modular/modint.hpp" namespace workspace { namespace internal { /** * @brief Base of modular arithmetic. * * @tparam Mod identifier, which represents modulus if positive * @tparam Storage Reserved size for inverse calculation */ template <auto Mod, unsigned Storage> struct modint_base { static_assert(is_integral_ext<decltype(Mod)>::value, "Mod must be integral type."); using mod_type = typename std::make_signed<typename std::conditional< 0 < Mod, typename std::add_const<decltype(Mod)>::type, decltype(Mod)>::type>::type; using value_type = typename std::decay<mod_type>::type; using mul_type = typename multiplicable_uint<value_type>::type; // Modulus static mod_type mod; static unsigned storage; constexpr static void reserve(unsigned __n) noexcept { storage = __n; } value_type value = 0; public: constexpr modint_base() noexcept = default; template <class _Tp, typename std::enable_if< is_integral_ext<_Tp>::value>::type * = nullptr> constexpr modint_base(_Tp __n) noexcept : value((__n %= mod) < 0 ? __n += mod : __n) {} constexpr modint_base(bool __n) noexcept : value(__n) {} constexpr operator value_type() const noexcept { return value; } template <class _Tp, typename = typename std::is_floating_point<_Tp>::value> constexpr explicit operator _Tp() const noexcept { return value; } constexpr static modint_base one() noexcept { return 1; } // unary operators {{ constexpr modint_base operator++(int) noexcept { modint_base __t{*this}; operator++(); return __t; } constexpr modint_base operator--(int) noexcept { modint_base __t{*this}; operator--(); return __t; } constexpr modint_base &operator++() noexcept { if (++value == mod) value = 0; return *this; } constexpr modint_base &operator--() noexcept { if (!value) value = mod; --value; return *this; } constexpr modint_base operator-() const noexcept { modint_base __t; __t.value = value ? mod - value : 0; return __t; } // }} unary operators // operator+= {{ constexpr modint_base &operator+=(const modint_base &__x) noexcept { if ((value += __x.value) >= mod) value -= mod; return *this; } template <class _Tp> constexpr typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type & operator+=(_Tp const &__x) noexcept { if (((value += __x) %= mod) < 0) value += mod; return *this; } // }} operator+= // operator+ {{ template <class _Tp> constexpr typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type operator+(_Tp const &__x) const noexcept { return modint_base{*this} += __x; } constexpr modint_base operator+(modint_base __x) const noexcept { return __x += *this; } template <class _Tp> constexpr friend typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type operator+(_Tp const &__x, modint_base __y) noexcept { return __y += __x; } // }} operator+ // operator-= {{ constexpr modint_base &operator-=(const modint_base &__x) noexcept { if ((value -= __x.value) < 0) value += mod; return *this; } template <class _Tp> constexpr typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type & operator-=(_Tp __x) noexcept { if (((value -= __x) %= mod) < 0) value += mod; return *this; } // }} operator-= // operator- {{ template <class _Tp> constexpr typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type operator-(_Tp const &__x) const noexcept { return modint_base{*this} -= __x; } constexpr modint_base operator-(const modint_base &__x) const noexcept { modint_base __t; if (((__t.value = value) -= __x.value) < 0) __t.value += mod; return __t; } template <class _Tp> constexpr friend typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type operator-(_Tp __x, const modint_base &__y) noexcept { if (((__x -= __y.value) %= mod) < 0) __x += mod; modint_base __t; __t.value = __x; return __t; } // }} operator- // operator*= {{ constexpr modint_base &operator*=(const modint_base &__x) noexcept { value = static_cast<value_type>(value * static_cast<mul_type>(__x.value) % mod); return *this; } template <class _Tp> constexpr typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type & operator*=(_Tp __x) noexcept { value = static_cast<value_type>( value * mul_type((__x %= mod) < 0 ? __x + mod : __x) % mod); return *this; } // }} operator*= // operator* {{ constexpr modint_base operator*(const modint_base &__x) const noexcept { modint_base __t; __t.value = static_cast<mul_type>(value) * __x.value % mod; return __t; } template <class _Tp> constexpr typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type operator*(_Tp __x) const noexcept { if ((__x %= mod) < 0) __x += mod; modint_base __t; __t.value = static_cast<mul_type>(value) * __x % mod; return __t; } template <class _Tp> constexpr friend typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type operator*(_Tp __x, const modint_base &__y) noexcept { if (!__x or !__y.value) return {}; if ((__x %= mod) < 0) __x += mod; mul_type __r(__x); modint_base __t; __t.value = static_cast<value_type &&>((__r *= __y.value) %= mod); return __t; } // }} operator* protected: static value_type _mem(value_type __x) { static std::vector<value_type> __m{0, 1}; static value_type __i = (__m.reserve(Storage), 1); while (__i < __x) { ++__i; __m.emplace_back(mod - mul_type(mod / __i) * __m[mod % __i] % mod); } return __m[__x]; } template <class _Tp> constexpr static typename std::enable_if<is_integral_ext<_Tp>::value, value_type>::type _div(mul_type __r, _Tp __x) noexcept { assert(__x != _Tp(0)); if (!__r) return 0; std::make_signed_t<_Tp> __v{}; bool __neg = __x < 0 ? __x = -__x, true : false; if (static_cast<decltype(storage)>(__x) < storage) __v = _mem(__x); else { decltype(__v) __y{mod}, __u{1}, __t; while (__x) __t = __y / __x, __y ^= __x ^= (__y -= __t * __x) ^= __x, __v ^= __u ^= (__v -= __t * __u) ^= __u; if (__y < 0) __neg ^= 1; } if (__neg) __v = 0 < __v ? mod - __v : -__v; else if (__v < 0) __v += mod; return __r == mul_type(1) ? static_cast<value_type>(__v) : static_cast<value_type>(__r * __v % mod); } public: // operator/= {{ constexpr modint_base &operator/=(const modint_base &__x) noexcept { if (value) value = _div(value, __x.value); return *this; } template <class _Tp> constexpr typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type & operator/=(_Tp __x) noexcept { if (value) value = _div(value, __x %= mod); return *this; } // }} operator/= // operator/ {{ constexpr modint_base operator/(const modint_base &__x) const noexcept { if (!value) return {}; modint_base __t; __t.value = _div(value, __x.value); return __t; } template <class _Tp> constexpr typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type operator/(_Tp __x) const noexcept { if (!value) return {}; modint_base __t; __t.value = _div(value, __x %= mod); return __t; } template <class _Tp> constexpr friend typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type operator/(_Tp __x, const modint_base &__y) noexcept { if (!__x) return {}; if ((__x %= mod) < 0) __x += mod; modint_base __t; __t.value = _div(__x, __y.value); return __t; } // }} operator/ constexpr modint_base inv() const noexcept { return _div(1, value); } template <class _Tp> friend constexpr typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type pow(modint_base b, _Tp e) noexcept { if (e < 0) { e = -e; b.value = _div(1, b.value); } modint_base __r; for (__r.value = 1; e; e >>= 1, b *= b) if (e & 1) __r *= b; return __r; } template <class _Tp> constexpr typename std::enable_if<is_integral_ext<_Tp>::value, modint_base>::type pow(_Tp __e) const noexcept { modint_base __r, b; __r.value = 1; for (b.value = __e < 0 ? __e = -__e, _div(1, value) : value; __e; __e >>= 1, b *= b) if (__e & 1) __r *= b; return __r; } template <class _Os> friend _Os &operator<<(_Os &__os, const modint_base &__x) noexcept { return __os << __x.value; } friend std::istream &operator>>(std::istream &__is, modint_base &__x) noexcept { intmax_t value; __x = (__is >> value, value); return __is; } }; template <auto Mod, unsigned Storage> typename modint_base<Mod, Storage>::mod_type modint_base<Mod, Storage>::mod = Mod > 0 ? Mod : 0; template <auto Mod, unsigned Storage> unsigned modint_base<Mod, Storage>::storage = Storage; } // namespace internal /** * @brief Modular arithmetic. * * @tparam Mod modulus * @tparam Storage Reserved size for inverse calculation */ template <auto Mod, unsigned Storage = 0, typename std::enable_if<(Mod > 0)>::type * = nullptr> using modint = internal::modint_base<Mod, Storage>; /** * @brief Runtime modular arithmetic. * * @tparam type_id uniquely assigned * @tparam Storage Reserved size for inverse calculation */ template <unsigned type_id = 0, unsigned Storage = 0> using modint_runtime = internal::modint_base<-(signed)type_id, Storage>; // #define modint_newtype modint_runtime<__COUNTER__> } // namespace workspace #line 11 "Library/src/modular/inverse.hpp" namespace workspace { // Modulus must be prime. template <class Modint> struct inverse_table { static_assert(std::is_same<std::nullptr_t, decltype((void *)Modint::mod, nullptr)>::value); using value_type = Modint; constexpr value_type operator()(int n) const { constexpr int_fast64_t mod = value_type::mod; assert(n %= mod); if (n < 0) n += mod; if (inv.empty()) inv = {1, mod != 1}; for (int m(inv.size()); m <= n; ++m) inv.emplace_back(mod / m * -inv[mod % m]); return inv[n]; } private: static std::vector<value_type> inv; }; template <class Modint> std::vector<Modint> inverse_table<Modint>::inv; } // namespace workspace #line 2 "Library/src/algebra/fps.hpp" /** * @file fps.hpp * @brief Formal Power Series */ #line 2 "Library/src/algebra/polynomial.hpp" /** * @file polynomial.hpp * @brief Polynomial */ #line 11 "Library/src/algebra/polynomial.hpp" #line 2 "Library/src/algebra/ntt.hpp" /** * @file ntt.hpp * @brief Number Theoretic Transform * @date 2021-02-20 * * */ #line 2 "Library/src/number_theory/ext_gcd.hpp" /** * @file ext_gcd.hpp * @brief Extended Euclidean Algorithm */ #line 9 "Library/src/number_theory/ext_gcd.hpp" #line 11 "Library/src/number_theory/ext_gcd.hpp" namespace workspace { /** * @param __a Integer * @param __b Integer * @return Pair of integers (x, y) s.t. ax + by = g = gcd(a, b), 0 <= x < * |b/g|, -|a/g| < y <= 0. Return (0, 0) if (a, b) = (0, 0). */ template <typename _T1, typename _T2> constexpr auto ext_gcd(_T1 __a, _T2 __b) { static_assert(is_integral_ext<_T1>::value); static_assert(is_integral_ext<_T2>::value); using result_type = typename std::make_signed< typename std::common_type<_T1, _T2>::type>::type; result_type a{__a}, b{__b}, p{1}, q{}, r{}, s{1}; // Euclidean algorithm while (b) { result_type t = a / b; r ^= p ^= r ^= p -= t * r; s ^= q ^= s ^= q -= t * s; b ^= a ^= b ^= a -= t * b; } // Normalize if (a < 0) p = -p, q = -q; if (p < 0) p += __b / a, q -= __a / a; return std::make_pair(p, q); } } // namespace workspace #line 2 "Library/src/number_theory/primitive_root.hpp" /** * @file primitive_root.hpp * @brief Primitive Root * @date 2020-12-28 */ #line 10 "Library/src/number_theory/primitive_root.hpp" namespace workspace { /** * @brief Compile time primitive root. * * @tparam __mod Positive integer * @return Minimum positive one if it exists. Otherwise 0. */ template <class Tp> constexpr typename std::enable_if<(is_integral_ext<Tp>::value), Tp>::type primitive_root(const Tp __mod) noexcept { assert(__mod > 0); using int_type = typename multiplicable_uint<Tp>::type; int_type __r = __mod, __p[16] = {}, *__q = __p; for (int_type __i = 2; __i <= __r / __i; ++__i) { if (__r % __i) continue; *__q++ = __i; while (!(__r % __i)) __r /= __i; } if (__r != 1) *__q++ = __r; int_type __tot = __mod; for (__q = __p; *__q; *__q++ = 0) (__tot /= *__q) *= *__q - 1; __r = __tot, __q = __p + 1, __p[0] = 1; for (int_type __i = 2; __i <= __r / __i; ++__i) { if (__r % __i) continue; *__q++ = __i; while (!(__r % __i)) __r /= __i; } if (__r != 1) *__q++ = __r; for (Tp __r = 1; __r != __mod; ++__r) { auto __cnt = 0; for (__q = __p; *__q; ++__q) { int_type __w = 1; for (int_type __e = __tot / *__q, __x = __r; __e; __e >>= 1, (__x *= __x) %= __mod) if (__e & 1) (__w *= __x) %= __mod; if (__w == 1 && ++__cnt > 1) break; } if (__cnt == 1) return __r; } return 0; }; } // namespace workspace #line 13 "Library/src/algebra/ntt.hpp" namespace workspace { namespace ntt_impl { /** * @see * https://github.com/atcoder/ac-library/blob/master/atcoder/convolution.hpp */ template <class _Tp> struct __coef { _Tp sum_e[30]; // sum_e[i] = ies[0] * ... * ies[i - 1] * es[i] constexpr __coef() : sum_e{} { if (_Tp::mod < 2) return; int cnt2 = __builtin_ctz(_Tp::mod - 1); _Tp e = 1; { auto p = (_Tp::mod - 1) >> cnt2; _Tp w = primitive_root(_Tp::mod); while (p) { if (p & 1) e *= w; p >>= 1; w *= w; } } _Tp ie = ext_gcd(decltype(_Tp::mod)(e), _Tp::mod).first; _Tp es[30] = {}, ies[30] = {}; // es[i]^(2^(2+i)) == 1 for (int i = cnt2; i >= 2; i--) { // e^(2^i) == 1 es[i - 2] = e; ies[i - 2] = ie; e *= e; ie *= ie; } _Tp now = 1; for (int i = 0; i <= cnt2 - 2; i++) { sum_e[i] = es[i] * now; now *= ies[i]; } } }; template <class _Tp> struct __icoef { _Tp sum_ie[30]; // sum_e[i] = ies[0] * ... * ies[i - 1] * es[i] constexpr __icoef() : sum_ie{} { if (_Tp::mod < 2) return; int cnt2 = __builtin_ctz(_Tp::mod - 1); _Tp e = 1; { auto p = (_Tp::mod - 1) >> cnt2; _Tp w = primitive_root(_Tp::mod); while (p) { if (p & 1) e *= w; p >>= 1; w *= w; } } _Tp ie = ext_gcd(decltype(_Tp::mod)(e), _Tp::mod).first; _Tp es[30] = {}, ies[30] = {}; // es[i]^(2^(2+i)) == 1 for (int i = cnt2; i >= 2; i--) { // e^(2^i) == 1 es[i - 2] = e; ies[i - 2] = ie; e *= e; ie *= ie; } _Tp now = 1; for (int i = 0; i <= cnt2 - 2; i++) { sum_ie[i] = ies[i] * now; now *= es[i]; } } }; template <class _Tp> struct __ipow2 { _Tp __ip2[30]; constexpr __ipow2() : __ip2{1, (1 + _Tp::mod) / 2} { for (size_t __i = 1; __i + 1 != std::size(__ip2); ++__i) __ip2[__i + 1] = __ip2[__i] * __ip2[1]; } }; template <class _FIter> constexpr void ntt(_FIter __first, _FIter __last) noexcept { using value_type = typename std::decay<decltype(*__first)>::type; constexpr __coef<value_type> _; auto __h = __builtin_ctz(std::distance(__first, __last)); for (ptrdiff_t __p = 1 << __h; __p >>= 1;) { value_type now = -1; auto __l = __first; for (size_t __s = 1 << __h; __l != __last; now *= _.sum_e[__builtin_ctz(--__s)]) { auto __r = __l + __p; for (auto __mid = __r; __l != __mid; ++__l, ++__r) { auto __tmp = *__l; *__l -= *__r *= now; *__r += __tmp; } __l = __r; } } } template <class _A> constexpr void ntt(_A &a) noexcept { ntt(std::begin(a), std::end(a)); } template <class _FIter> constexpr void intt(_FIter __first, _FIter __last) noexcept { using value_type = typename std::decay<decltype(*__first)>::type; constexpr __icoef<value_type> _; auto __h = __builtin_ctz(std::distance(__first, __last)); for (ptrdiff_t __p = 1; __p >> __h ^ 1; __p <<= 1) { value_type inow = 1; auto __l = __first; for (size_t __s = 1 << __h; __l != __last; inow *= _.sum_ie[__builtin_ctz(--__s)]) { auto __r = __l + __p; for (auto __mid = __r; __l != __mid; ++__l, ++__r) { auto __tmp = (*__l - *__r) * inow; *__l += *__r; *__r = __tmp; } __l = __r; } } constexpr __ipow2<value_type> __; while (__first != __last) *--__last *= __.__ip2[__h]; } // namespace ntt_impl template <class _A> constexpr void intt(_A &a) noexcept { intt(std::begin(a), std::end(a)); } } // namespace ntt_impl using ntt_impl::intt; using ntt_impl::ntt; } // namespace workspace #line 14 "Library/src/algebra/polynomial.hpp" namespace workspace { /** * @brief Polynomial class. * * @tparam _Tp Ring structure * @tparam _Conv_threshold Threshold for convolution method */ template <class _Tp, std::size_t _Conv_threshold = 64> class polynomial : public std::vector<_Tp> { using vec = std::vector<_Tp>; using poly = polynomial; public: using vec::vec; using size_type = typename vec::size_type; protected: void _erase_leading_zeros() noexcept { auto __i = vec::_M_impl._M_finish; while (__i != vec::_M_impl._M_start && *(__i - 1) == _Tp(0)) --__i; vec::_M_erase_at_end(__i); } template <class _Iter> void _dft(_Iter __first, _Iter __last) const noexcept { if constexpr (has_mod<_Tp>::value) ntt(__first, __last); else { // fft(__first, __last); assert(0); // Not implemented! } } template <class _Iter> void _idft(_Iter __first, _Iter __last) const noexcept { if constexpr (has_mod<_Tp>::value) intt(__first, __last); else { // ifft(__first, __last); assert(0); // Not implemented! } } void _conv_naive(const poly& __x) noexcept { if (vec::_M_impl._M_start == vec::_M_impl._M_finish) return; if (__x._M_impl._M_start == __x._M_impl._M_finish) { vec::_M_erase_at_end(vec::_M_impl._M_start); return; } vec::_M_default_append(__x._M_impl._M_finish - __x._M_impl._M_start - 1); for (auto __i = vec::_M_impl._M_finish; __i-- != vec::_M_impl._M_start;) { auto __j = __i, __k = __x._M_impl._M_start; *__i *= *__k++; while (__j != vec::_M_impl._M_start && __k != __x._M_impl._M_finish) *__i += *--__j * *__k++; } } void _conv_dft(poly&& __x) noexcept { if constexpr (has_mod<_Tp>::value) _conv_ntt(std::move(__x)); else { // _conv_fft(std::move(__x)); assert(0); // Not implemented! } } void _conv_fft(poly&& __x) noexcept; void _conv_ntt(poly&& __x) noexcept { size_type __n = vec::_M_impl._M_finish - vec::_M_impl._M_start, __m = __x._M_impl._M_finish - __x._M_impl._M_start, __len = 1 << (32 - __builtin_clz(__n + __m - 1)); vec::_M_default_append(__len - __n); __x._M_default_append(__len - __m); ntt(vec::_M_impl._M_start, vec::_M_impl._M_finish); ntt(__x._M_impl._M_start, __x._M_impl._M_finish); for (auto __i = vec::_M_impl._M_start, __j = __x._M_impl._M_start; __i != vec::_M_impl._M_finish; ++__i, ++__j) *__i *= std::move(*__j); intt(vec::_M_impl._M_start, vec::_M_impl._M_finish); vec::_M_erase_at_end(vec::_M_impl._M_start + __n + __m - 1); } /** * @brief * * @param __x * @return Degree of __x. */ size_type _divmod_naive(const poly& __x) { auto __xfin = __x._M_impl._M_finish; auto __xlen = __x.size(); while (__xfin != __x._M_impl._M_start && *(__xfin - 1) == _Tp(0)) --__xfin, --__xlen; assert(__xlen != 0); _erase_leading_zeros(); auto __p = vec::_M_impl._M_finish; while (size_type(__p - vec::_M_impl._M_start) >= __xlen) { --__p; auto __src = __xfin; auto __dst = __p; *__dst /= *--__src; while (__src != __x._M_impl._M_start) *--__dst -= *--__src * *__p; } return std::min<size_type>(__xlen - 1, __p - vec::_M_impl._M_start); } void _div_naive(const poly& __x) { operator>>=(_divmod_naive(__x)); } void _div_doubling(poly&& __x) noexcept { _erase_leading_zeros(); __x._erase_leading_zeros(); auto __n = vec::_M_impl._M_finish - vec::_M_impl._M_start; auto __m = __x._M_impl._M_finish - __x._M_impl._M_start; if (__n < __m) vec::clear(); else { assert(__m != 0); std::reverse(__x._M_impl._M_start, __x._M_impl._M_finish); __x = __x.inv(__n - __m + 1); std::reverse(vec::_M_impl._M_start, vec::_M_impl._M_finish); vec::_M_erase_at_end(vec::_M_impl._M_finish - (__m - 1)); operator*=(__x).resize(__n - __m + 1); std::reverse(vec::_M_impl._M_start, vec::_M_impl._M_finish); } } public: /** * @return Degree of %polynomial. Return -1 if it equals zero. */ size_type deg() const noexcept { return vec::size() - 1; } /** * @param __i Not exceeding the degree. * @return Coefficient of x^i. */ typename vec::reference operator[](size_type __i) noexcept { assert(__i < vec::size()); return *(vec::_M_impl._M_start + __i); } /** * @param __i Not exceeding the degree. * @return Coefficient of x^i. */ typename vec::const_reference operator[](size_type __i) const noexcept { assert(__i < vec::size()); return *(vec::_M_impl._M_start + __i); } /** * @brief Evaluate at given point. */ _Tp eval(const _Tp& __a) const noexcept { _Tp __v(0), __p(1); for (auto __i = vec::_M_impl._M_start; __i != vec::_M_impl._M_finish; ++__i, __p *= __a) __v += *__i * __p; return __v; } /** * @brief In-place multipoint evaluation. */ template <class _Iter> _Iter eval(_Iter __first, _Iter __last) const noexcept { return eval(__first, __last, __first); } /** * @brief Multipoint evaluation. */ template <class _IIter, class _OIter> _OIter eval(_IIter __first, _IIter __last, _OIter __result) const noexcept { size_type __n = std::distance(__first, __last); if (!__n) return __result; auto __tree = new poly[__n << 1]; for (auto __p = __tree + __n; __first != __last; ++__p, ++__first) *__p = {-*__first, 1}; for (size_type __i = __n; --__i;) __tree[__i] = __tree[__i << 1] * __tree[__i << 1 | 1]; __tree[1] = operator%(std::move(__tree[1])); for (size_type __i = 2; __i != __n << 1; __i += 2) __tree[__i] = __tree[__i >> 1] % std::move(__tree[__i]), __tree[__i | 1] = std::move(__tree[__i >> 1] %= std::move(__tree[__i | 1])); for (size_type __i = 0; __i != __n; ++__i) *__result++ = std::move(*__tree[__n + __i]._M_impl._M_start); delete[] __tree; return __result; } /** * @brief Conversion to bool. * * @return Whether the %polynomial is not zero. */ operator bool() const noexcept { auto __first = vec::_M_impl._M_start, __last = vec::_M_impl._M_finish; while (__first != __last) if (*__first++ != _Tp(0)) return true; return false; } bool operator==(const poly& __x) const noexcept { auto __first1 = vec::_M_impl._M_start, __last1 = vec::_M_impl._M_finish; auto __first2 = __x._M_impl._M_start, __last2 = __x._M_impl._M_finish; if (__last1 - __first1 < __last2 - __first2) { while (__first1 != __last1) if (*__first1++ != *__first2++) return false; while (__first2 != __last2) if (*__first2++ != _Tp(0)) return false; } else { while (__first2 != __last2) if (*__first1++ != *__first2++) return false; while (__first1 != __last1) if (*__first1++ != _Tp(0)) return false; } return true; } bool operator!=(const poly& __x) const noexcept { return !operator==(__x); } /** * @brief Multiply by x^i. */ poly& operator<<=(size_type __i) noexcept { vec::insert(vec::begin(), __i, _Tp(0)); return *this; } /** * @brief Divide by x^i. */ poly& operator>>=(size_type __i) noexcept { vec::_M_erase_at_end( std::move(vec::_M_impl._M_start + std::min(__i, vec::size()), vec::_M_impl._M_finish, vec::_M_impl._M_start)); return *this; } /** * @brief Multiply by x^i. */ poly operator<<(size_type __i) const noexcept { return poly(*this).operator<<=(__i); } /** * @brief Divide by x^i. */ poly operator>>(size_type __i) const noexcept { return poly(*this).operator>>=(__i); } poly operator+() const noexcept { return *this; } poly operator-() const noexcept { poly __x = *this; for (auto __i = __x._M_impl._M_start; __i != __x._M_impl._M_finish; ++__i) *__i = -*__i; return __x; } poly& operator+=(const poly& __x) noexcept { if (vec::size() < __x.size()) vec::_M_default_append(__x.size() - vec::size()); for (auto __i = vec::_M_impl._M_start, __j = __x._M_impl._M_start; __j != __x._M_impl._M_finish; ++__i, ++__j) *__i += *__j; _erase_leading_zeros(); return *this; } poly& operator+=(const _Tp& __c) noexcept { if (__c != static_cast<_Tp>(0)) { if (vec::_M_impl._M_start == vec::_M_impl._M_finish) vec::emplace_back(__c); else *vec::_M_impl._M_start += __c, _erase_leading_zeros(); } return *this; } poly& operator-=(const poly& __x) noexcept { if (vec::size() < __x.size()) vec::_M_default_append(__x.size() - vec::size()); for (auto __i = vec::_M_impl._M_start, __j = __x._M_impl._M_start; __j != __x._M_impl._M_finish; ++__i, ++__j) *__i -= *__j; _erase_leading_zeros(); return *this; } poly& operator-=(const _Tp& __c) noexcept { if (__c != static_cast<_Tp>(0)) { if (vec::_M_impl._M_start == vec::_M_impl._M_finish) vec::emplace_back(-__c); else *vec::_M_impl._M_start -= __c, _erase_leading_zeros(); } return *this; } poly& operator*=(const poly& __x) noexcept { std::min(vec::size(), __x.size()) > _Conv_threshold ? _conv_dft(poly(__x)) : _conv_naive(this == std::addressof(__x) ? poly(__x) : __x); return *this; } poly& operator*=(poly&& __x) noexcept { std::min(vec::size(), __x.size()) > _Conv_threshold ? _conv_dft(std::move(__x)) : _conv_naive(__x); return *this; } poly& operator*=(const _Tp& __c) noexcept { if (__c == static_cast<_Tp>(0)) vec::_M_erase_at_end(vec::_M_impl._M_start); else for (auto __i = vec::_M_impl._M_start; __i != vec::_M_impl._M_finish; ++__i) *__i *= __c; return *this; } poly& operator/=(const _Tp& __c) noexcept { assert(__c != static_cast<_Tp>(0)); for (auto __i = vec::_M_impl._M_start; __i != vec::_M_impl._M_finish; ++__i) *__i /= __c; return *this; } poly rev() const noexcept { return rev(vec::size()); } poly rev(size_type __n) const noexcept { poly __r(__n); auto __src = vec::_M_impl._M_start; auto __dst = __r._M_impl._M_finish; for (size_type __i = std::min(__n, vec::size()); __i; --__i) *--__dst = *__src++; return __r; } poly inv() const noexcept { return inv(vec::size()); } /** * @brief Multiplicative inverse modulo x^n. * * @param __n Degree of modulus * @return */ poly inv(size_type __n) const noexcept { if (!__n) return {}; assert(*vec::_M_impl._M_start != _Tp(0)); size_type __len = 1; while (__len < __n) __len <<= 1; poly __y(__len); auto __xp = new _Tp[__len], __yp = __y._M_impl._M_start, __zp = new _Tp[__len]; *__yp = _Tp(1) / *vec::_M_impl._M_start; for (size_type __i = 1; __i != __len; __i <<= 1) { std::fill(std::copy_n(__yp, __i, __zp), __zp + (__i << 1), _Tp(0)); _dft(__zp, __zp + (__i << 1)); std::fill(std::copy_n(vec::_M_impl._M_start, std::min(__i << 1, vec::size()), __xp), __xp + (__i << 1), _Tp(0)); _dft(__xp, __xp + (__i << 1)); for (size_type __j = 0; __j != (__i << 1); ++__j) __xp[__j] *= -__zp[__j]; _idft(__xp, __xp + (__i << 1)); std::fill(std::move(__xp + __i, __xp + (__i << 1), __xp), __xp + (__i << 1), _Tp(0)); _dft(__xp, __xp + (__i << 1)); for (size_type __j = 0; __j != (__i << 1); ++__j) __xp[__j] *= static_cast<_Tp&&>(__zp[__j]); _idft(__xp, __xp + (__i << 1)); std::move(__xp, __xp + __i, __yp + __i); } delete[] __xp; delete[] __zp; __y._M_erase_at_end(__yp + __n); return __y; } poly& operator/=(const poly& __x) noexcept { if (__x.size() > _Conv_threshold) _div_doubling(poly(__x)); else _div_naive(__x); return *this; } poly& operator/=(poly&& __x) noexcept { if (__x.size() > _Conv_threshold) _div_doubling(std::move(__x)); else _div_naive(__x); return *this; } poly& operator%=(const poly& __x) noexcept { if (__x.size() > _Conv_threshold) return operator-=(__x.operator*(operator/(__x))); vec::_M_erase_at_end(vec::_M_impl._M_start + _divmod_naive(__x)); return *this; } template <class _T> poly operator+(_T&& __x) const noexcept { return poly(*this).operator+=(std::forward<_T>(__x)); } template <class _T> poly operator-(_T&& __x) const noexcept { return poly(*this).operator-=(std::forward<_T>(__x)); } template <class _T> poly operator*(_T&& __x) const noexcept { return poly(*this).operator*=(std::forward<_T>(__x)); } template <class _T> poly operator/(_T&& __x) const noexcept { return poly(*this).operator/=(std::forward<_T>(__x)); } template <class _T> poly operator%(_T&& __x) const noexcept { return poly(*this).operator%=(std::forward<_T>(__x)); } std::pair<poly, poly> divmod(const poly& __x) const { if (__x.size() > _Conv_threshold) return {operator/(__x), operator%(__x)}; poly __rem(*this); auto __p = __rem._M_impl._M_start + __rem._divmod_naive(__x); poly __quot(__p, __rem._M_impl._M_finish); __rem._M_erase_at_end(__p); return {__quot, __rem}; } /** * @brief Differentiate. * * @return Derivative. */ poly deriv() const noexcept { if (auto __s = vec::_M_impl._M_start, __f = vec::_M_impl._M_finish; __s != __f) { poly __der(++__s, __f); __s = __der._M_impl._M_start, __f = __der._M_impl._M_finish; for (_Tp __i(1); __s != __f; ++__s, __i += 1) *__s *= __i; __der._erase_leading_zeros(); return __der; } return {}; } /** * @brief Differentiate at given point. * * @return Derivative coefficient. */ _Tp deriv(const _Tp& __a) const noexcept { _Tp __der(0); if (auto __s = vec::_M_impl._M_start, __f = vec::_M_impl._M_finish; __s != __f) for (_Tp __i(1), __p(1); ++__s != __f; __i += 1, __p *= __a) __der += *__s * __i * __p; return __der; } /** * @brief Integrate. * * @return Integral indefinite at the degrees divisible by the characteristic * of `_Tp`. Coefficients are set as 0 there. */ poly integ() const noexcept { if (auto __s = vec::_M_impl._M_start, __f = vec::_M_impl._M_finish; __s != __f) { poly __int(__f - __s + 1); __f = std::copy(__s, __f, __int._M_impl._M_start + 1); __s = __int._M_impl._M_start + 1; for (_Tp __i(1); __s != __f; ++__s, __i += 1) __i == _Tp(0) ? assert(*__s == _Tp(0)) : void(*__s /= __i); return __int; } return {}; } /** * @brief Integrate in given range. * * @return Definite integral over [0, __a]. */ _Tp integ(const _Tp& __a) const noexcept { _Tp __int(0); auto __s = vec::_M_impl._M_start, __f = vec::_M_impl._M_finish; for (_Tp __p(__a), __i(1); __s != __f; ++__s, __p *= __a, __i += 1) __int += *__s / __i * __p; return __int; } /** * @brief Integrate in given range. * * @return Definite integral over [__a, __b]. */ _Tp integ(const _Tp& __a, const _Tp& __b) const noexcept { _Tp __int(0); auto __s = vec::_M_impl._M_start, __f = vec::_M_impl._M_finish; for (_Tp __pa(__a), __pb(__b), __i(1); __s != __f; ++__s, __pa *= __a, __pb *= __b, __i += 1) __int += *__s / __i * (__pb - __pa); return __int; } }; } // namespace workspace #line 9 "Library/src/algebra/fps.hpp" namespace workspace { namespace fps_impl { template <class _Tp> decltype(auto) _exp(_Tp __x) { return exp(__x); } template <class _Tp> decltype(auto) _log(_Tp __x) { return log(__x); } template <class _Tp, class _Tp2> decltype(auto) _pow(_Tp __x, _Tp2 __y) { return pow(__x, __y); } template <class _T, typename = void> struct has_exp : std::false_type {}; template <class _T> struct has_exp<_T, std::enable_if_t<std::is_same< _T, decltype(exp(std::declval<_T>()))>::value>> : std::true_type {}; template <class _T, typename = void> struct has_log : std::false_type {}; template <class _T> struct has_log<_T, std::enable_if_t<std::is_same< _T, decltype(log(std::declval<_T>()))>::value>> : std::true_type {}; template <class _T1, class _T2, typename = void> struct has_pow : std::false_type {}; template <class _T1, class _T2> struct has_pow< _T1, _T2, std::enable_if_t<std::is_same< _T1, decltype(pow(std::declval<_T1>(), std::declval<_T2>()))>::value>> : std::true_type {}; } // namespace fps_impl template <class _Tp> class formal_power_series : public polynomial<_Tp> { using poly = polynomial<_Tp>; using fps = formal_power_series; public: using poly::poly; using typename poly::size_type; formal_power_series(const poly& __x) : poly(__x) {} formal_power_series(poly&& __x) : poly(std::move(__x)) {} formal_power_series& operator=(const poly& __x) { poly::operator=(__x); return *this; } formal_power_series& operator=(poly&& __x) { poly::operator=(std::move(__x)); return *this; } // private: using poly::_dft; using poly::_idft; fps _exp(size_type __n) const noexcept { if (__n == 1) return {1}; size_type __len = 2; while (__len < __n) __len <<= 1; fps __w(__len); auto __wp = __w.data(), __xp = new _Tp[__len], __yp = new _Tp[__len]{}, __zp = new _Tp[__len]; *__wp = 1, *(__wp + (__len >> 1)) = 1; for (size_type __i = 1;; __i <<= 1) { std::copy_n(__wp, __i, __xp); std::copy_n(__wp + (__len >> 1), __i, __yp); std::fill(std::copy_n(poly::_M_impl._M_start + 1, std::min(__i << 1, poly::size()) - 1, __zp), __zp + (__i << 1), _Tp(0)); for (size_type __j = 1; __j != (__i << 1); ++__j) __zp[__j - 1] *= __j; _dft(__xp, __xp + (__i << 1)); _dft(__yp, __yp + (__i << 1)); _dft(__zp, __zp + (__i << 1)); for (size_type __j = 0; __j != (__i << 1); ++__j) __zp[__j] *= __xp[__j]; _idft(__zp, __zp + (__i << 1)); std::fill(std::move(__zp + __i - 1, __zp + (__i << 1) - 1, __zp), __zp + (__i << 1), _Tp(0)); _dft(__zp, __zp + (__i << 1)); for (size_type __j = 0; __j != (__i << 1); ++__j) __zp[__j] *= __yp[__j]; _idft(__zp, __zp + (__i << 1)); for (size_type __j = 0; __j != __i; ++__j) __zp[__j] /= __j + __i, __zp[__j + __i] = 0; _dft(__zp, __zp + (__i << 1)); for (size_type __j = 0; __j != (__i << 1); ++__j) __zp[__j] *= static_cast<_Tp&&>(__xp[__j]); _idft(__zp, __zp + (__i << 1)); std::move(__zp, __zp + __i, __wp + __i); if (__i << 1 == __len) break; std::copy_n(__wp, __i << 1, __zp); _dft(__zp, __zp + (__i << 1)); for (size_type __j = 0; __j != (__i << 1); ++__j) __zp[__j] *= -__yp[__j]; _idft(__zp, __zp + (__i << 1)); std::fill(std::move(__zp + __i, __zp + (__i << 1), __zp), __zp + (__i << 1), _Tp(0)); _dft(__zp, __zp + (__i << 1)); for (size_type __j = 0; __j != (__i << 1); ++__j) __zp[__j] *= static_cast<_Tp&&>(__yp[__j]); _idft(__zp, __zp + (__i << 1)); std::move(__zp, __zp + __i, __wp + (__len >> 1) + __i); } delete[] __xp; delete[] __yp; delete[] __zp; __w.resize(__n); return __w; } public: fps sqrt(size_type __n) const noexcept {} /** * @brief Exponential. */ fps exp() const noexcept { return exp(poly::size()); } /** * @brief Exponential modulo x^n. * * @param __n Degree of modulus * @return */ fps exp(size_type __n) const noexcept { if (!__n) return {}; if (poly::empty()) return {1}; if constexpr (fps_impl::has_exp<_Tp>::value) { auto __x = _exp(__n); __x *= fps_impl::_exp(*poly::_M_impl._M_start); return __x; } assert(*poly::_M_impl._M_start == _Tp(0)); return _exp(__n); } /** * @brief Logarithm. */ fps log() const noexcept { return log(poly::size()); } /** * @brief Logarithm modulo x^n. * * @param __n Degree of modulus * @return */ fps log(size_type __n) const noexcept { assert(!poly::empty()); if constexpr (fps_impl::has_log<_Tp>::value) { assert(*poly::_M_impl._M_start != _Tp(0)); if (!__n) return {}; if (__n == 1) return {0}; auto __x = *this; auto __c = *poly::_M_impl._M_start; for (auto&& __a : __x) __a /= __c; __x = __x.inv(__n - 1) * __x.deriv(); __x.resize(__n); std::move_backward(__x.begin(), __x.begin() + __n - 1, __x.begin() + __n); for (size_type __i = 1; __i != __n; ++__i) __x[__i] /= __i; __x[0] = fps_impl::_log(__c); return __x; } assert(*poly::_M_impl._M_start == _Tp(1)); if (!__n) return {}; if (__n == 1) return {0}; auto __x = poly::inv(__n - 1).operator*(poly::deriv()); __x.resize(__n); std::move_backward(__x.begin(), __x.begin() + __n - 1, __x.begin() + __n); for (size_type __i = 1; __i != __n; ++__i) __x[__i] /= __i; __x[0] = 0; return __x; } template <class _Tp2> fps pow(_Tp2 __e) const noexcept { return pow(__e, poly::size()); } template <class _Tp2> fps pow(_Tp2 __e, size_type __n) const noexcept { if (!__n) return {}; if (__e == _Tp2(0)) { fps __x(__n, 0); __x[0] = 1; return __x; } auto __first = poly::_M_impl._M_start; if constexpr (is_integral_ext<_Tp2>::value) { if (__e > 0) { size_type __o = 0; while (__first != poly::_M_impl._M_finish && *__first == _Tp(0)) ++__first, ++__o; if ((__n - 1) / __e < __o || __first == poly::_M_impl._M_finish) return fps(__n, 0); __o *= __e; fps __x(__first, poly::_M_impl._M_finish); if constexpr (fps_impl::has_pow<_Tp, _Tp2>::value) { assert(*__first != _Tp(0)); auto __c = *__first; for (auto&& __a : __x) __a /= __c; __x = __x.log(__n - __o); __x *= __e; __x = __x.exp(__n - __o); __c = fps_impl::_pow(__c, __e); for (auto&& __a : __x) __a *= __c; __x.insert(__x.begin(), __o, _Tp(0)); return __x; } assert(*__first == _Tp(1)); __x = __x.log(__n - __o); __x *= __e; __x = __x.exp(__n - __o); __x.insert(__x.begin(), __o, _Tp(0)); return __x; } } if constexpr (fps_impl::has_pow<_Tp, _Tp2>::value) { assert(*__first != _Tp(0)); fps __x(__first, poly::_M_impl._M_finish); auto __c = *__first; for (auto&& __a : __x) __a /= __c; __x = __x.log(__n); __x *= __e; __x = __x.exp(__n); __c = fps_impl::_pow(__c, __e); for (auto&& __a : __x) __a *= __c; return __x; } assert(*__first == _Tp(1)); auto __x = log(__n); __x *= __e; return __x.exp(__n); } }; } // namespace workspace #line 29 "atcoder-workspace/2.cc" namespace workspace { // Capacity: 256 template <class T> struct iterable { T data[256] = {}; T *be, *en; iterable() : be(data), en(data) {} template <class It> iterable(It fi, It se) : be(data), en(data) { copy(fi, se); } template <class It> void copy(It fi, It se) { en = be; while (fi != se) *(en++) = *(fi++); } void pop() { --en; } void top() { return *(en - 1); } void clear() { en = be; } void insert(T x) { *(en++) = x; } void erase(T x) { T *p = en; while (*--p != x) ; // pass *p = *--en; } bool empty() const { return be == en; } size_t size() const { return en - be; } T *begin() { return be; } T *end() { return en; } }; struct A { iterable<A *> iter; void insert(A *x) { iter.insert(x); } }; void main() { // start here! using mint = modint<998244353>; int n; cin >> n; formal_power_series<mint> f(n); cin >> f; cout << f.log() << "\n"; } } // namespace workspace