// #define _GLIBCXX_DEBUG

#ifndef LIST__HPP_
#define LIST__HPP_

#include <algorithm>
#include <iostream>
#include <cassert>
#include <functional>
#include <utility>
using uint = unsigned int;
using ll = long long;
using namespace std;


template<typename T> concept has_identity = requires (){T();};
template<typename T> concept has_plus_operator = requires (T a, T b){a+b;};
template<typename T> concept has_compare_operator = requires (T a, T b){a < b;};
template<typename T> concept has_equal_operator = requires (T a, T b){a == b; a != b;};


template<typename T>
class list{
    public:

    using value_type = T;

    constexpr list(): sz_(0), max_sz_(0), A(nullptr){}
    constexpr list(uint sz): sz_(sz), max_sz_(sz){
        if (sz == 0){A = nullptr; return;}
        if (sz_ > 2147483647){cerr << "ValueError: too large or negative size" << endl; assert(false);}
        A = new T[sz_]{};
        if (A == nullptr){cerr << "MomoryError: failed to allocate memory" << endl; assert(false);}
    }
    constexpr list(uint sz, T elem): sz_(sz), max_sz_(sz){
        if (sz == 0){A = nullptr; return;}
        if (sz_ > 2147483647){cerr << "ValueError: too large or negative size" << endl; assert(false);}
        A = new T[sz_]{};
        if (A == nullptr){cerr << "MemoryError: failed to allocate memory" << endl; assert(false);}
        for (int i = 0; i < sz_; i++){
            A[i] = elem;
        }
    }
    constexpr list(initializer_list<T> init): sz_(init.size()), max_sz_(init.size()){
        A = new T[sz_]{};
        if (A == nullptr){cerr << "MemoryError: failed to allocate memory" << endl; assert(false);}
        for (int i = 0; i < sz_; i++){
            A[i] = *(init.begin()+i);
        }
    };
    // constexpr list<ll>(const range& r){
    //     if (r.step_ == 0){cerr << "ValueError: step argument must not be zero" << endl; assert(false);}
    //     if (r.step_ > 0){
    //         int imax = r.begin_ >= r.end_ ? 0 : (r.end_-r.begin_+r.step_-1)/r.step_;
    //         sz_ = imax;
    //         max_sz_ = imax;
    //         A = new T[sz_]{};
    //         if (A == nullptr){cerr << "MemoryError: failed to allocate memory" << endl; assert(false);}
    //         for (int i = 0; i < imax; i++){
    //             int target = r.begin_+r.step_*i;
    //             A[i] = target;
    //         }    
    //     }
    //     else{
    //         int imax = r.end_ >= r.begin_ ? 0 : (r.begin_-r.end_-r.step_-1)/(-r.step_);
    //         sz_ = imax;
    //         max_sz_ = imax;
    //         A = new T[sz_]{};
    //         if (A == nullptr){cerr << "MemoryError: failed to allocate memory" << endl; assert(false);}
    //         for (int i = 0; i < imax; i++){
    //             int target = r.begin_+r.step_*i;
    //             A[i] = target;
    //         }
    //     }
    // }
    constexpr list(const list<T>& other): sz_(other.sz_), max_sz_(other.max_sz_){
        A = new T[max_sz_]{};
        for (int i = 0; i < sz_; i++){
            A[i] = other.A[i];
        }
    }
    constexpr list(list<T>&& other) noexcept : sz_(other.sz_), max_sz_(other.max_sz_){
        A = other.A;
        other.sz_ = 0;
        other.max_sz_ = 0;
        other.A = nullptr;
    }
    constexpr ~list() {
        delete[] A;
    }

    class iterator{
        public:
        using difference_type = ll;
        using value_type = T;
        using pointer = T*;
        using reference = T&;
        using iterator_category = std::random_access_iterator_tag;
        
        iterator() = delete;
        constexpr iterator(const list* listptr_init, const int &idx_init): listptr(listptr_init), idx(idx_init){}

        inline constexpr value_type& operator* () const {return listptr->A[this->check_idx()];}
        inline constexpr value_type* operator->() const {return &listptr->A[this->check_idx()];}

        inline constexpr bool operator==(const iterator &other) const {return idx == other.idx;}
        inline constexpr bool operator!=(const iterator &other) const {return idx != other.idx;}
        inline constexpr bool operator<(const iterator &other) const {return idx < other.idx;}
        inline constexpr bool operator>(const iterator &other) const {return idx > other.idx;}
        inline constexpr bool operator<=(const iterator &other) const {return idx <= other.idx;}
        inline constexpr bool operator>=(const iterator &other) const {return idx >= other.idx;}

        inline constexpr iterator operator++(){
            this->idx++;
            return *this;
        }
        inline constexpr iterator operator++(int){
            iterator temp = *this;
            this->idx++;
            return temp;
        }
        inline constexpr iterator operator--(){
            this->idx--;
            return *this;
        }
        inline constexpr iterator operator--(int){
            iterator temp = *this;
            this->idx--;
            return temp;
        }

        inline constexpr iterator& operator+=(int amount){
            this->idx += amount;
            return *this;
        }
        inline constexpr iterator operator+(int amount) const {
            iterator temp = *this;
            temp.idx += amount;
            return temp;
        }
        inline constexpr iterator& operator-=(const int &amount){
            this->idx -= amount;
            return *this;
        }
        inline constexpr iterator operator-(const int &amount) const{
            iterator temp = *this;
            temp.idx -= amount;
            return temp;
        }

        inline constexpr difference_type operator-(const iterator &other) const{
            return idx - other.idx;
        }

        
        inline constexpr int check_idx() const {
            if (idx < 0 || idx >= listptr->sz_){
                cerr << "IteratorIndexError: index out of range. current size: " << listptr->sz_ << " target index: " << idx << endl; assert(false);
            }
            return idx;
        }

        private:
        
        const list* listptr;
        int idx;
    };

    inline constexpr iterator begin(){return iterator(this, 0);}
    inline constexpr iterator begin() const {return iterator(this, 0);}
    inline constexpr iterator end(){return iterator(this, sz_);}
    inline constexpr iterator end() const {return iterator(this, sz_);}

    inline constexpr list& operator=(const list& other){
        if (this == &other){return *this;}
        sz_ = other.sz_;
        max_sz_ = other.max_sz_;
        delete[] A;
        A = new T[max_sz_]{};
        for (int i = 0; i < sz_; i++){
            A[i] = other.A[i];
        }
        return *this;
    }
    inline constexpr list& operator=(list&& other) noexcept {
        if (this == &other){return *this;}
        delete[] A;
        sz_ = other.sz_;
        max_sz_ = other.max_sz_;
        A = other.A;
        other.sz_ = 0;
        other.max_sz_ = 0;
        other.A = nullptr;
        return *this;
    }
    inline constexpr bool operator==(const list& other) const {
        if (sz_ != other.sz_){return false;}
        bool res = true;
        for (int i = 0; i < sz_; i++){
            if (A[i] != other.A[i]){
                res = false;
                break;
            }
        }
        return res;
    }
    inline constexpr bool operator!=(const list& other) const {
        return !this->operator==(other);
    }
    inline constexpr bool operator<(const list& other) const {
        int lcp_back = 0;
        for (int i = 0, sz2 = min(sz_, other.sz_); i < sz2; i++){
            if (A[i] != other.A[i]){break;}
            lcp_back++;
        }
        if (lcp_back == sz_ && lcp_back == other.sz_){return false;}
        else if (lcp_back == sz_){return true;}
        else if (lcp_back == other.sz_){return false;}
        else{return A[lcp_back] < other.A[lcp_back];}
    }
    inline constexpr bool operator<=(const list& other) const {
        int lcp_back = 0;
        for (int i = 0, sz2 = min(sz_, other.sz_); i < sz2; i++){
            if (A[i] != other.A[i]){break;}
            lcp_back++;
        }
        if (lcp_back == sz_ && lcp_back == other.sz_){return true;}
        else if (lcp_back == sz_){return true;}
        else if (lcp_back == other.sz_){return false;}
        else{return A[lcp_back] < other.A[lcp_back];}
    }
    inline constexpr bool operator>(const list& other) const {return other < *this;}
    inline constexpr bool operator>=(const list& other) const {return other <= *this;}
    inline constexpr list& operator+=(const list& other){
        for (int i = 0, sz2 = other.size(); i < sz2; i++){
            append(other.A[i]);
        }
        return *this;
    }
    inline constexpr list operator+(const list& other) const {
        return list(*this) += other;
    }
    inline constexpr list& operator+=(const T& val){
        append(val);
        return *this;
    }
    inline constexpr list operator+(const T& val) const {
        return list(*this) += val;
    }
    
    

    inline constexpr void append(const T& val){
        if (sz_ < max_sz_){
            A[sz_] = val;
            sz_++;
        }
        else{
            if (max_sz_ < 2){max_sz_ = 2;}
            else{max_sz_ += max_sz_>>1;}
            auto B = new T[max_sz_]{};
            B[sz_] = val;
            for (int i = 0; i < sz_; i++){B[i] = std::move(A[i]);}
            sz_++;
            delete[] A;
            A = B;
        }
    }
    inline constexpr void append(T&& val){
        if (sz_ < max_sz_){
            A[sz_] = std::move(val);
            sz_++;
        }
        else{
            if (max_sz_ < 2){max_sz_ = 2;}
            else{max_sz_ += max_sz_>>1;}
            auto B = new T[max_sz_]{};
            B[sz_] = std::move(val);
            for (int i = 0; i < sz_; i++){B[i] = std::move(A[i]);}
            sz_++;
            delete[] A;
            A = B;
        }
    }
    inline constexpr T pop(){
        if (sz_ == 0){cerr << "IndexError: pop from empty list" << endl; assert(false);}
        sz_--;
        return std::move(A[sz_]);
    }
    inline constexpr void reserve(uint new_sz){
        if (sz_ > 2147483647){cerr << "ValueError: too large or negative size" << endl; assert(false);}
        if ((int)new_sz <= max_sz_){return;}
        max_sz_ = new_sz;
        auto B = new T[new_sz]{};
        if (B == nullptr){cerr << "MemoryError: failed to reallocate memory" << endl; assert(false);}
        for (int i = 0; i < sz_; i++){B[i] = std::move(A[i]);}
        delete[] A;
        A = B;
    }
    inline constexpr void resize(uint new_sz, T elem = T()){
        if (sz_ > 2147483647){cerr << "ValueError: too large or negative size" << endl; assert(false);}
        if ((int)new_sz <= max_sz_){
            while (sz_ > (int)new_sz){pop();}
            for (int i = sz_, new_sz_int = new_sz; i < new_sz_int; i++){A[i] = elem;}
            sz_ = new_sz;
            return;
        }
        max_sz_ = new_sz + (new_sz >> 1);
        auto B = new T[max_sz_]{};
        if (B == nullptr){cerr << "MemoryError: failed to reallocate memory" << endl; assert(false);}
        for (int i = 0; i < sz_; i++){B[i] = std::move(A[i]);}
        for (int i = sz_, new_sz_int = new_sz; i < new_sz_int; i++){B[i] = elem;}
        sz_ = new_sz;
        delete[] A;
        A = B;
    }
    template <class InputIterator> inline constexpr void assign(InputIterator first, InputIterator last){
        reserve(last-first);
        while (first != last){append(*first); first++;}
    }
    inline constexpr void assign(uint n, const T& val){
        reserve(n);
        sz_ = n;
        for (uint i = 0; i < n; i++){
            A[i] = val;
        }
    }
    inline constexpr iterator insert(iterator pos, const T& val){
        int idx = pos.check_idx();
        append(val);
        for (int i = idx; i < sz_-1; i++){swap(A[i], A[sz_-1]);}
        return pos;
    }
    // inline constexpr void insert(iterator pos, iterator begin, iterator end)
    inline constexpr iterator erase(iterator pos){
        int idx = pos.check_idx();
        for (int i = idx; i < sz_-1; i++){
            swap(A[i], A[i+1]);
        }
        sz_ -= 1;
        return pos;
    }
    inline constexpr iterator erase(iterator start, iterator end){
        int startidx = start.check_idx(), endidx = end.check_idx();
        if (startidx > endidx){cerr << "IndexError: invalid interval: [start, end) = []" << startidx << ", " << endidx << ")" << endl; assert(false);}
        for (int i = startidx; i < endidx; i++){
            A[i].~T();
        }
        if (startidx < endidx){
            for (int i = endidx; i < sz_; i++){
                swap(A[i+startidx-endidx], A[i]);
            }
        }
        sz_ -= endidx-startidx;
        return end;
    }
    inline constexpr void clear(){sz_ = 0; max_sz_ = 0; delete[] A;}
    inline constexpr const T& operator[](int idx) const {
        check_and_adjust_idx(idx);
        return A[idx];
    }
    inline constexpr T& operator[](int idx){
        check_and_adjust_idx(idx);
        return A[idx];
    }
    inline constexpr T& front(){return (*this)[0];}
    inline constexpr const T& front() const {return (*this)[0];}
    inline constexpr T& back(){return (*this)[-1];}
    inline constexpr const T& back() const {return (*this)[-1];}
    inline constexpr list<T> operator[](int left, int right, int step = 1) const {
        if (step == 0){cerr << "ValueError: step argument must not be zero" << endl; assert(false);}
        if (step > 0){
            int imax = left >= right ? 0 : (right-left+step-1)/step;
            list<T> B(imax);
            for (int i = 0; i < imax; i++){
                int target = left+step*i;
                if (target < 0 || target >= sz_){cerr << "IndexError: slice index out of range. current size: " << sz_ << " target index: " << target << endl; assert(false);}
                B[i] = A[target];
            }
            return std::move(B);
        }
        else{
            int imax = right >= left ? 0 : (left-right-step-1)/(-step);
            list<T> B(imax);
            for (int i = 0; i < imax; i++){
                int target = left+step*i;
                if (target < 0 || target >= sz_){cerr << "IndexError: slice index out of range. current size: " << sz_ << " target index: " << target << endl; assert(false);}
                B[i] = A[target];
            }
            return std::move(B);
        }
    }
    inline constexpr T* data() {return empty() ? nullptr : &A[0];}
    inline constexpr const T* data() const {return empty() ? nullptr : &A[0];}
    inline constexpr T sum() const requires has_plus_operator<T> && has_identity<T>{
        auto res = T();
        for (int i = 0; i < sz_; i++){
            res += A[i];
        }
        return res;
    }
    template<typename U> inline constexpr U sum(U init) const requires has_plus_operator<T>{
        for (int i = 0; i < sz_; i++){
            init += A[i];
        }
        return init;
    }
    inline constexpr void sort(bool reverse = false) requires has_compare_operator<T> {
        if (reverse){
            std::sort(begin(), end(), std::move([](const T& a, const T& b){return b < a;}));
        }
        else{
            std::sort(begin(), end());
        }
    }
    inline constexpr void sort(const function<bool(T,T)>& f) requires has_compare_operator<T> {
        std::sort(begin(), end(), f);
    }
    inline constexpr void reverse(){
        std::reverse(begin(), end());
    }
    inline constexpr ll count(const T& val){
        ll res = 0;
        for (int i = 0; i < sz_; i++){res += A[i] == val;}
        return res;
    }
    inline constexpr void apply(const function<T(T)>& f){for (int i = 0; i < sz_; i++){A[i] = f(A[i]);}}
    inline constexpr list<pair<int, T>> RLE(){
        list<pair<int, T>> res;
        for (int l = 0, r = 0; l < sz_;){
            while (r < sz_ && A[l] == A[r]){
                r++;
            }
            res.append(make_pair(r-l, A[l]));
            l = r;
        }
        return res;
    }
    inline constexpr iterator lower_bound(const T& val) const requires has_compare_operator<T> {return std::lower_bound(begin(), end(), val);}
    inline constexpr iterator upper_bound(const T& val) const requires has_compare_operator<T> {return std::upper_bound(begin(), end(), val);}
    inline constexpr ll size() const {return sz_;}
    inline constexpr ll maxsize() const {return max_sz_;}
    inline constexpr bool empty() const {return sz_ == 0;}

    private:

    inline constexpr void check_and_adjust_idx(int& idx) const {
        if (idx >= sz_){cerr << "IndexError: index out of range. current size: " << sz_ << " target index: " << idx << endl; assert(false);}
        if (idx < 0){
            if (idx < -sz_){cerr << "IndexError: index out of range. current size: " << sz_ << " target index: " << idx << endl; assert(false);}
            idx += sz_;
        }
    }

    int sz_;
    int max_sz_;
    T* A;
};


#endif /* LIST__HPP_ */
#include <iostream>
#include <iterator>
#include <type_traits>
#include <utility>
#include <string>
#include <cmath>
#include <unordered_set>
#include <unordered_map>
#include <queue>
#include <algorithm>
#include <iomanip>
#include <cassert>
#include <functional>
#include <random>
#include <bitset>

// #include <boost/multiprecision/cpp_int.hpp>



using namespace std;
using ll = long long;
using lll = __int128_t;
using ull = unsigned long long;
#ifdef BOOST_VERSION
using bll = boost::multiprecision::cpp_int;
#endif
using ld = long double;
using pii = array<int,2>;
using pll = array<ll,2>;
using plll = array<lll,2>;
using pld = array<ld, 2>;
using tll = array<ll,3>;

#define vall(A) A.begin(), A.end()
template<typename T> inline bool btest(T K, int i){return K&(1ull<<i);}
#ifndef LEN_DEFINED
#define LEN_DEFINED
template<typename T> inline constexpr ll len(const T& A){return ssize(A);}
#endif
constexpr ll pow2ll[63] = {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,32768,65536,131072,262144,524288,1048576,2097152,4194304,8388608,16777216,33554432,67108864,134217728,268435456,536870912,1073741824,2147483648,4294967296,8589934592,17179869184,34359738368,68719476736,137438953472,274877906944,549755813888,1099511627776,2199023255552,4398046511104,8796093022208,17592186044416,35184372088832,70368744177664,140737488355328,281474976710656,562949953421312,1125899906842624,2251799813685248,4503599627370496,9007199254740992,18014398509481984,36028797018963968,72057594037927936,144115188075855872,288230376151711744,576460752303423488,1152921504606846976,2305843009213693952,4611686018427387904};
constexpr ll pow10ll[19] = {1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000,10000000000,100000000000,1000000000000,10000000000000,100000000000000,1000000000000000,10000000000000000,100000000000000000,1000000000000000000};
constexpr ll di[4] = {0,1,0,-1};
constexpr ll di8[8] = {0,1,1,1,0,-1,-1,-1};
constexpr ll dj[4] = {1,0,-1,0};
constexpr ll dj8[8] = {1,1,0,-1,-1,-1,0,1};


#ifndef FASTIO__HPP_
#define FASTIO__HPP_

#include <array>
#include <bit>
#include <charconv>
#include <cstddef>
#include <cstdint>
#include <cstdio>
#include <cstring>
#include <iterator>
#include <string>
#include <string_view>
#include <tuple>
#include <type_traits>
#include <utility>

// C++20. Integer input must fit the destination type. For modint input, the
// magnitude must fit uint64_t; negative values are reduced modulo the modulus.
// FILE* is borrowed: it must outlive the reader/writer. Do not mix buffered
// readers on one stream. The default input buffer is intended for batch judges;
// interactive programs need BasicInput<1> and explicit output flush().
namespace fastio {
namespace detail {
template<class T> inline constexpr bool integer = std::is_integral_v<T>;
template<class T> struct unsigned_type { using type = std::make_unsigned_t<T>; };
template<> struct unsigned_type<bool> { using type = unsigned char; };
#if defined(__SIZEOF_INT128__)
template<> inline constexpr bool integer<__int128_t> = true;
template<> inline constexpr bool integer<__uint128_t> = true;
template<> struct unsigned_type<__int128_t> { using type = __uint128_t; };
template<> struct unsigned_type<__uint128_t> { using type = __uint128_t; };
#endif
template<class T> using unsigned_t = typename unsigned_type<T>::type;
template<class T> inline constexpr bool signed_integer = integer<T> && T(-1) < T(0);
template<class T> concept modular = requires(T x) { x.val; } &&
    (requires { T::mod_value; } || requires { T::mod(); });
template<class T> constexpr auto modulus() {
    if constexpr (requires { T::mod_value; }) return T::mod_value;
    else return T::mod();
}
template<class T> concept indexed_range = requires(T& x) { std::size(x); x[0]; };
template<class T> concept tuple_like = requires { std::tuple_size<T>::value; };
template<class> inline constexpr bool unsupported = false;
inline constexpr auto digit_pairs = [] {
    std::array<char, 200> result{};
    for (int i = 0; i < 100; ++i) {
        result[2*i] = char('0'+i/10);
        result[2*i+1] = char('0'+i%10);
    }
    return result;
}();
inline constexpr auto digit_quads = [] {
    std::array<std::array<char,4>,10000> result{};
    for (unsigned i = 0; i < 10000; ++i) {
        result[i] = {char('0'+i/1000),char('0'+i/100%10),
                     char('0'+i/10%10),char('0'+i%10)};
    }
    return result;
}();
inline char* leading_quad(char* p, unsigned value) {
    if (value >= 1000) { std::memcpy(p,digit_quads[value].data(),4); return p+4; }
    if (value >= 100) { *p++ = char('0'+value/100); value %= 100; }
    else if (value < 10) { *p++ = char('0'+value); return p; }
    std::memcpy(p,digit_pairs.data()+2*value,2);
    return p+2;
}
inline char* decimal32(char* p, std::uint32_t value) {
    if (value >= 100000000) {
        p = leading_quad(p,value/100000000);
        value %= 100000000;
        std::memcpy(p,digit_quads[value/10000].data(),4);
        std::memcpy(p+4,digit_quads[value%10000].data(),4);
        return p+8;
    }
    if (value >= 10000) {
        p = leading_quad(p,value/10000);
        std::memcpy(p,digit_quads[value%10000].data(),4);
        return p+4;
    }
    return leading_quad(p,value);
}
template<class T> inline constexpr bool output32 = [] {
    if constexpr (integer<T>) return sizeof(T) <= 4 && !std::is_same_v<T,char>;
    else if constexpr (modular<T>) return integer<decltype(T{}.val)> && sizeof(T{}.val) <= 4;
    else return false;
}();
template<class T> inline char* decimal_small(char* p, const T& value) {
    if constexpr (modular<T>) return decimal_small(p,value.val);
    else {
        std::uint32_t u = static_cast<std::uint32_t>(value);
        if constexpr (signed_integer<T>) {
            if (value < 0) { *p++ = '-'; u = 0-u; }
        }
        return decimal32(p,u);
    }
}
} // namespace detail

template<std::size_t BufferSize = 1 << 18>
class BasicInput {
    static_assert(BufferSize > 0);
    std::FILE* file_;
    char buffer_[BufferSize+1];
    std::size_t position_ = 0, size_ = 0;
    bool good_ = true;

    int peek() {
        if (position_ == size_) {
            size_ = std::fread(buffer_, 1, BufferSize, file_);
            buffer_[size_] = '\0';
            position_ = 0;
            if (size_ == 0) return EOF;
        }
        return static_cast<unsigned char>(buffer_[position_]);
    }
    int skip_space() {
        int c;
        while ((c = peek()) != EOF && c <= ' ') ++position_;
        return c;
    }
    template<class U, bool SingleBlock = false> bool magnitude(U& value, bool& negative) {
        int c = skip_space();
        if (c == EOF) return false;
        negative = c == '-';
        if (c == '-' || c == '+') { ++position_; c = peek(); }
        if (static_cast<unsigned>(c-'0') > 9) return false;
        value = 0;
        do {
            const char* p = buffer_+position_;
            if constexpr (sizeof(U) >= 4 && std::endian::native == std::endian::little) {
                while (buffer_+size_-p >= 8) {
                    std::uint64_t block;
                    std::memcpy(&block,p,8);
                    // All eight bytes must be ASCII decimal digits.
                    if (((block+0x4646464646464646ULL) | (block-0x3030303030303030ULL)) & 0x8080808080808080ULL) break;
                    block -= 0x3030303030303030ULL;
                    block = (block*10+(block>>8)) & 0x00ff00ff00ff00ffULL;
                    block = (block*100+(block>>16)) & 0x0000ffff0000ffffULL;
                    block = (block*10000+(block>>32)) & 0xffffffffULL;
                    value = value*100000000+static_cast<std::uint32_t>(block);
                    p += 8;
                    if constexpr (SingleBlock || sizeof(U) == 4) break;
                }
            }
            unsigned digit;
            // The sentinel stops this loop without an EOF/buffer check for
            // every digit. Only a token crossing a refill needs the outer loop.
            while ((digit = static_cast<unsigned>(static_cast<unsigned char>(*p)-'0')) <= 9) {
                value = value*10+digit;
                ++p;
            }
            position_ = static_cast<std::size_t>(p-buffer_);
            if (position_ != size_) break;
            c = peek();
        } while (static_cast<unsigned>(c-'0') <= 9);
        return true;
    }
    template<class T> bool read_one(T& value) {
        if constexpr (std::is_same_v<T, char>) {
            const int c = skip_space();
            if (c == EOF) return false;
            value = static_cast<char>(c);
            ++position_;
            return true;
        } else if constexpr (detail::integer<T>) {
            detail::unsigned_t<T> u;
            bool negative;
            if (!magnitude(u,negative)) return false;
            if constexpr (detail::signed_integer<T>) {
                // Never negate the minimum signed value, even for int128.
                value = negative && u != 0 ? -static_cast<T>(u-1)-1 : static_cast<T>(u);
            } else {
                if (negative || (std::is_same_v<T,bool> && u > 1)) return false;
                value = static_cast<T>(u);
            }
            return true;
        } else if constexpr (detail::modular<T>) {
            std::uint64_t u;
            bool negative;
            // Most modular coefficients have at most ten digits. After the
            // first eight, parse the short tail without another wide probe.
            if (!magnitude<std::uint64_t,true>(u,negative)) return false;
            const auto mod = detail::modulus<T>();
            if (u >= mod) u %= mod;
            value.val = negative && u != 0 ? mod-u : u;
            return true;
        } else if constexpr (std::is_same_v<T,std::string>) {
            int c = skip_space();
            if (c == EOF) return false;
            value.clear();
            do {
                value.push_back(static_cast<char>(c));
                ++position_;
                c = peek();
            } while (c != EOF && c > ' ');
            return true;
        } else if constexpr (std::is_floating_point_v<T>) {
            std::string token;
            if (!read_one(token)) return false;
            const char* first = token.data();
            const char* last = first+token.size();
            if (first != last && *first == '+') ++first;
            const auto result = std::from_chars(first,last,value);
            return result.ec == std::errc{} && result.ptr == last;
        } else if constexpr (detail::indexed_range<T>) {
            const auto n = std::size(value);
            if constexpr (requires { std::data(value); }) {
                auto* data = std::data(value);
                for (decltype(std::size(value)) i = 0; i < n; ++i)
                    if (!read_one(data[i])) return false;
            } else {
                for (decltype(std::size(value)) i = 0; i < n; ++i)
                    if (!read_one(value[i])) return false;
            }
            return true;
        } else if constexpr (detail::tuple_like<T>) {
            return std::apply([&](auto&... x) { return (read_one(x) && ...); },value);
        } else {
            static_assert(detail::unsupported<T>, "Unsupported fastio input type");
        }
    }
public:
    explicit BasicInput(std::FILE* file = stdin) : file_(file) {}
    BasicInput(const BasicInput&) = delete;
    BasicInput& operator=(const BasicInput&) = delete;

    template<class... T> bool read(T&... values) {
        if (!good_) return false;
        return good_ = (read_one(values) && ...);
    }
    template<class T> BasicInput& operator>>(T& value) { read(value); return *this; }
    explicit operator bool() const { return good_; }
};

template<std::size_t BufferSize = 1 << 18>
class BasicOutput {
    static_assert(BufferSize > 0);
    std::FILE* file_;
    char buffer_[BufferSize];
    std::size_t size_ = 0;
    bool good_ = true;

    bool drain() {
        if (!good_) return false;
        std::size_t done = 0;
        while (done < size_) {
            const auto n = std::fwrite(buffer_+done,1,size_-done,file_);
            if (n == 0) { size_ = 0; return good_ = false; }
            done += n;
        }
        size_ = 0;
        return true;
    }
    bool put(char c) {
        if (!good_ || (size_ == BufferSize && !drain())) return false;
        buffer_[size_++] = c;
        return true;
    }
    bool bytes(const char* data, std::size_t count) {
        if (!good_) return false;
        while (count != 0) {
            if (size_ == BufferSize && !drain()) return false;
            const auto room = BufferSize-size_;
            const auto n = count < room ? count : room;
            std::memcpy(buffer_+size_,data,n);
            size_ += n; data += n; count -= n;
        }
        return true;
    }
    template<class U> bool unsigned_number(U value, bool negative) {
        if constexpr (sizeof(U) <= 4 && BufferSize >= 12) {
            if (!good_ || (BufferSize-size_ < 12 && !drain())) return false;
            char* p = buffer_+size_;
            if (negative) *p++ = '-';
            p = detail::decimal32(p,value);
            size_ = static_cast<std::size_t>(p-buffer_);
            return true;
        }
        // 128-bit unsigned integers need at most 39 decimal digits.
        char digits[3*sizeof(U)+2];
        char* const end = digits+sizeof(digits);
        char* p = end;
        while (value >= 100) {
            const U q = value/100;
            const unsigned r = static_cast<unsigned>(value-q*100);
            p -= 2;
            std::memcpy(p,detail::digit_pairs.data()+2*r,2);
            value = q;
        }
        if (value < 10) *--p = char('0'+static_cast<unsigned>(value));
        else {
            p -= 2;
            std::memcpy(p,detail::digit_pairs.data()+2*static_cast<unsigned>(value),2);
        }
        if (negative) *--p = '-';
        if constexpr (BufferSize >= sizeof(digits)) {
            if (!good_ || (BufferSize-size_ < sizeof(digits) && !drain())) return false;
            auto count = static_cast<std::size_t>(end-p);
            char* dest = buffer_+size_;
            size_ += count;
            // Constant-size copies avoid a library call per integer. Each copy
            // stays inside the initialized digits and the reserved buffer area.
            if constexpr (sizeof(digits) >= 8)
                while (count >= 8) {
                    std::memcpy(dest,p,8);
                    dest += 8; p += 8; count -= 8;
                }
            if (count >= 4) { std::memcpy(dest,p,4); dest += 4; p += 4; count -= 4; }
            if (count >= 2) { std::memcpy(dest,p,2); dest += 2; p += 2; count -= 2; }
            if (count) *dest = *p;
            return true;
        }
        return bytes(p,static_cast<std::size_t>(end-p));
    }
    template<class T> bool write_one(const T& value) {
        if constexpr (std::is_same_v<T,char>) {
            return put(value);
        } else if constexpr (detail::integer<T>) {
            using U = detail::unsigned_t<T>;
            U u = static_cast<U>(value);
            bool negative = false;
            if constexpr (detail::signed_integer<T>) {
                negative = value < 0;
                if (negative) u = U(0)-u;
            }
            return unsigned_number(u,negative);
        } else if constexpr (detail::modular<T>) {
            return write_one(value.val);
        } else if constexpr (std::is_convertible_v<const T&,std::string_view>) {
            const std::string_view text(value);
            return bytes(text.data(),text.size());
        } else if constexpr (std::is_floating_point_v<T>) {
            char text[128];
            const auto result = std::to_chars(text,text+sizeof(text),value);
            if (result.ec != std::errc{}) return false;
            return bytes(text,static_cast<std::size_t>(result.ptr-text));
        } else if constexpr (detail::indexed_range<const T>) {
            const auto n = std::size(value);
            using Element = std::remove_cvref_t<decltype(value[0])>;
            if constexpr (detail::indexed_range<Element> || detail::tuple_like<Element>){
                for (decltype(std::size(value)) i = 0; i < n; ++i) {
                    if (i != 0 && !put('\n')) return false;
                    if (!write_one(value[i])) return false;
                }
                return true;
            }
            if constexpr (BufferSize >= 32 && detail::output32<Element> && requires { std::data(value); }) {
                const auto* data = std::data(value);
                char* p = buffer_+size_;
                for (decltype(std::size(value)) i = 0; i < n; ++i) {
                    if (buffer_+BufferSize-p < 12) {
                        size_ = static_cast<std::size_t>(p-buffer_);
                        if (!drain()) return false;
                        p = buffer_;
                    }
                    if (i != 0) *p++ = ' ';
                    p = detail::decimal_small(p,data[i]);
                }
                size_ = static_cast<std::size_t>(p-buffer_);
                return true;
            }
            for (decltype(std::size(value)) i = 0; i < n; ++i) {
                if (i != 0 && !put(' ')) return false;
                if (!write_one(value[i])) return false;
            }
            return true;
        } else if constexpr (detail::tuple_like<T>) {
            [[maybe_unused]] bool first = true;
            return std::apply([&](const auto&... x) {
                return (((std::exchange(first,false) || put(' ')) && write_one(x)) && ...);
            },value);
        } else {
            static_assert(detail::unsupported<T>, "Unsupported fastio output type");
        }
    }
public:
    explicit BasicOutput(std::FILE* file = stdout) : file_(file) {}
    BasicOutput(const BasicOutput&) = delete;
    BasicOutput& operator=(const BasicOutput&) = delete;
    ~BasicOutput() { flush(); }

    // write concatenates its arguments; print separates them with spaces and
    // appends a newline. Ranges and tuples use spaces between their elements.
    template<class... T> bool write(const T&... values) {
        if (!good_) return false;
        return good_ = (write_one(values) && ...);
    }
    template<class... T> bool print(const T&... values) {
        if (!good_) return false;
        [[maybe_unused]] bool first = true;
        good_ = (((std::exchange(first,false) || put(' ')) && write_one(values)) && ...);
        return good_ && put('\n');
    }
    template<class... T> bool printflush(const T&... values) {
        if (!good_) return false;
        [[maybe_unused]] bool first = true;
        good_ = (((std::exchange(first,false) || put(' ')) && write_one(values)) && ...);
        return good_ && put('\n') && flush();
    }
    bool flush() { return drain() && (good_ = std::fflush(file_) == 0); }
    template<class T> BasicOutput& operator<<(const T& value) { write(value); return *this; }
    explicit operator bool() const { return good_; }
};


#ifdef LOCAL_a1073741824
using FastInput = BasicInput<1>;
#else
using FastInput = BasicInput<>;
#endif
using FastOutput = BasicOutput<>;
// template<class... T> bool read(T&... values) { return in.read(values...); }
// template<class... T> bool write(const T&... values) { return out.write(values...); }
// template<class... T> bool print(const T&... values) { return out.print(values...); }
// template<class... T> bool printflush(const T&... values) { return out.print(values...) && out.flush(); }
// inline bool flush() { return out.flush(); }
} // namespace fastio

#endif /* FASTIO__HPP_ */

#ifndef FPS_OPERATION_HPP_
#define FPS_OPERATION_HPP_

#include <array>

#include <algorithm>
#include <functional>
#if defined(__x86_64__) || defined(__i386__)
#include <immintrin.h>
#endif
#ifndef DEQUE__HPP_
#define DEQUE__HPP_

#include <deque>
#include <iostream>
#include <iterator>
using namespace std;
using ll = long long;
using uint = unsigned int;

template<typename T>
class delist{
    public:

    inline constexpr delist(){A = deque<T>();}
    inline constexpr delist(uint sz, const T& val = T()){A = deque<T>(sz, val);}
    template <class InputIterator> inline constexpr delist(InputIterator first, InputIterator last){A = deque<T>(first, last);}

    inline constexpr void append(const T& val){A.push_back(val);}
    inline constexpr void appendleft(const T& val){A.push_front(val);}
    inline constexpr T pop(){
        if (A.empty()){cerr << "IndexError: pop from empty list" << endl;}
        auto temp = std::move(A.back());
        A.pop_back();
        return std::move(temp);
    }
    inline constexpr T popleft(){
        if (A.empty()){cerr << "IndexError: pop from empty list" << endl;}
        auto temp = std::move(A.front());
        A.pop_front();
        return std::move(temp);
    }
    inline constexpr T& front(){return A.front();}
    inline constexpr const T& front() const {return A.front();}
    inline constexpr T& back(){return A.back();}
    inline constexpr const T& back() const {return A.back();}

    inline constexpr _Deque_iterator<T, T&, T*> begin(){return A.begin();}
    inline constexpr _Deque_iterator<T, T&, T*> begin() const {return A.begin();}
    inline constexpr _Deque_iterator<T, T&, T*> end(){return A.end();}
    inline constexpr _Deque_iterator<T, T&, T*> end() const {return A.end();}


    inline constexpr ll size(){return A.size();}
    inline constexpr bool empty(){return A.empty();}

    private:


    deque<T> A;
};

#endif /* DEQUE__HPP_ */
#ifndef FPS__HPP_
#define FPS__HPP_

#include <iostream>
#include <iterator>

#include <cassert>
#ifndef MODINT__HPP_
#define MODINT__HPP_

#include <iostream>
#ifndef MATH_FUNCTION_HPP_
#define MATH_FUNCTION_HPP_

#include <array>
#include <cmath>
#include <cassert>
using namespace std;
using ll = long long;
using ld = long double;
using ull = unsigned long long;
using ulll = __uint128_t;
#ifdef BOOST_VERSION
using bll = boost::multiprecision::cpp_int;
#endif


constexpr ll getdigit(ull x){
    if (x < 1000000000){
        if (x < 100000){
            return x == 0 ? 0 : x < 10 ? 1 : x < 100 ? 2 : x < 1000 ? 3 : x < 10000 ? 4 : 5;
        }
        else{
            return x < 1000000 ? 6 : x < 10000000 ? 7 : x < 100000000 ? 8 : 9;
        }
    }
    else{
        if (x < 100000000000000){
            return x < 10000000000 ? 10 : x < 100000000000 ? 11 : x < 1000000000000 ? 12 : x < 10000000000000 ? 13 : 14;
        }
        else{
            return x < 1000000000000000 ? 15 : x < 10000000000000000 ? 16 : x < 100000000000000000 ? 17 : x < 1000000000000000000 ? 18 : x < 10000000000000000000u ? 19 : 20;
        }
    }
}

/// @brief a^bをmで割った余りを返す。bに関して対数時間で計算できる
constexpr ll modpow(ll a, ull b, const ll m){
    ll t = a%m;
    ll ans = (m == 1 ? 0 : 1);
    while (b > 0){
        if (b&1){
            ans = (ans*t)%m;
        }
        b >>= 1;
        t = (t*t)%m;
    }
    return ans;
}

/// @brief a^bをmで割った余りを返す。bに関して対数時間で計算できる。mはコンパイル時に決定している必要がある
template<ll m> constexpr ll modpow(ll a, ull b){
    ll t = a%m;
    ll ans = (m == 1 ? 0 : 1);
    while (b > 0){
        if (b&1){
            ans = (ans*t)%m;
        }
        b >>= 1;
        t = (t*t)%m;
    }
    return ans;
}

/// @brief a^nを返す。bに関して対数時間で計算できる。
template<typename T> constexpr T powll(T a, T n){
    T t = a;
    T ans = 1;
    while (n > 0){
        if (n%2){
            ans *= t;
        }
        n >>= 1;
        t *= t;
    }
    return ans;
}
/// @brief a^nを返す。bに関して対数時間で計算できる
template<> constexpr ll powll(ll a, ll n){
    ll t = a;
    ll ans = 1;
    while (n > 0){
        if (n&1){
            ans *= t;
        }
        n >>= 1;
        t *= t;
    }
    return ans;
}


/// @brief floor(sqrt(N))を返す。1.5×10^19まで対応
constexpr ll isqrt(ull N){
    assert(N <= 15000000000000000000ull);
    ull ret = sqrt(N);
    while (ret*ret > N){
        ret--;
    }
    while ((ret+1)*(ret+1) <= N){
        ret++;
    }
    return ret;
}
#ifdef BOOST_VERSION
/// @brief floor(sqrt(N))を返す。多倍長整数に対応
constexpr bll isqrt_large(bll N){
    bll ret = sqrt(N);
    while (ret*ret > N){
        ret--;
    }
    while ((ret+1)*(ret+1) <= N){
        ret++;
    }
    return ret;
}
#endif

/// @brief floor(log_a(L))を返す
constexpr ll ilog(ll a, ll L){
    __int128_t t = 1;
    ll ans = 0;
    while (t <= L){
        ans++;
        t *= a;
    }
    return ans-1;
}

/// @brief 有理数のfloorを求める。 floor(y/x)
template<typename T> constexpr inline T floor2(T y, T x){
    if ((x^y) > 0){
        x = x > 0 ? x : -x;
        y = y > 0 ? y : -y;
        return y/x;
    }
    else if ((x^y) < 0){
        x = x > 0 ? x : -x;
        y = y > 0 ? y : -y;
        return -((y+x-1)/x);
    }
    else{
        return y/x;
    }
}
/// @brief 有理数のceilを求める。 ceil(y/x)
template<typename T> constexpr inline T ceil2(T y, T x){
    if ((x^y) > 0){
        x = x > 0 ? x : -x;
        y = y > 0 ? y : -y;
        return (y+x-1)/x;
    }
    else if ((x^y) < 0){
        x = x > 0 ? x : -x;
        y = y > 0 ? y : -y;
        return -(y/x);
    }
    else{
        return y/x;
    }
}

/// @brief 一次不定方程式ax+by=gcd(a,b)の解を1つ見つける
/// @param a `a>=0`である必要がある
/// @param b `b>=0`である必要がある
/// @return {x,y,gcd(a,b)}
template<typename T>
constexpr array<T,3> axby1(T a, T b){
    T x = 1, y = 0;
    T z = 0, w = 1;
    T tmp = 0;
    while (b){
        T p = a/b, q = a%b;
        tmp = x - y * p; x = y; y = tmp;
        tmp = z - w * p; z = w; w = tmp;
        a = b; b = q;
    }
    return {x, z, a};
}

/// @brief 1/a mod Mを求める
template<typename T, typename U>
constexpr T inverse_mod(T a, U M){
    auto temp = axby1(a,(T)M);
    assert(temp[2] == 1);
    return (M+temp[0])%M;
}

/// @brief sqrt(a) mod Mを求める。ないなら-1が返される。
template<ll M>
inline constexpr ll cipolla(ll a){
    a %= M;
    if (M == 2) return a;
	if (a == 0) return 0;
    ll z = (M-1)/2;
    if (modpow<M>(a, z) != 1){return -1;}
    int b = 0;
    while (modpow<M>((b*b+M-a)%M, z) == 1){
        b++;
    }
    array<ll,2> x{1,0};
    array<ll,2> y{b, 1};
    ll w = (b*b+M-a)%M;
    z++;
    while (z){
        if (z&1){
            ll temp = x[0];
            x[0] = x[0]*y[0]%M+x[1]*y[1]%M*w%M;
            if (x[0] >= M){x[0] -= M;}
            x[1] = temp*y[1]%M+x[1]*y[0]%M;
            if (x[1] >= M){x[1] -= M;}
        }
        ll temp = y[0];
        y[0] = y[0]*y[0]%M+y[1]*y[1]%M*w%M;
        if (y[0] >= M){y[0] -= M;}
        y[1] = 2*temp*y[1]%M;
        z >>= 1;
    }
    return x[0];
}
inline constexpr ll cipolla(ll a, const ll M){
    a %= M;
    if (M == 2) return a;
	if (a == 0) return 0;
    ll z = (M-1)/2;
    if (modpow(a, z, M) != 1){return -1;}
    int b = 0;
    while (modpow((b*b+M-a)%M, z, M) == 1){
        b++;
    }
    array<ll,2> x{1,0};
    array<ll,2> y{b, 1};
    ll w = (b*b+M-a)%M;
    z++;
    while (z){
        if (z&1){
            ll temp = x[0];
            x[0] = x[0]*y[0]%M+x[1]*y[1]%M*w%M;
            if (x[0] >= M){x[0] -= M;}
            x[1] = temp*y[1]%M+x[1]*y[0]%M;
            if (x[1] >= M){x[1] -= M;}
        }
        ll temp = y[0];
        y[0] = y[0]*y[0]%M+y[1]*y[1]%M*w%M;
        if (y[0] >= M){y[0] -= M;}
        y[1] = 2*temp*y[1]%M;
        z >>= 1;
    }
    return x[0];
}

/// @brief x以下の最大の2冪を返す。0は0が返る。 
inline constexpr ull lowerpow2(ull x){
    if (x == 0){return 0;}
    return 1ull<<(63-__builtin_clzll(x));
}
/// @brief x以上の最小の2冪を返す。0は0が返る。 
inline constexpr ull upperpow2(ull x){
    if (x == 0){return 0;}
    if (x == 1){return 1;}
    return 1ull<<(64-__builtin_clzll(x-1));
}
/// @brief xのpopcountを求める 
inline constexpr int popcount(ull x){
    return __builtin_popcountll(x);
}
/// @brief xのbit lengthを求める。
inline constexpr int bit_length(ull x){
    if (x == 0){return 0;}
    return 64-__builtin_clzll(x);
}




#endif /* MATH_FUNCTION_HPP_ */

using ll = long long;
using uint = unsigned int;
using ull = unsigned long long;
using namespace std;

template <uint M>
struct constant_modint {
    static constexpr uint mod_value = M;
    uint val;
    constant_modint() : val(0) {}
    template <class T> constant_modint(T x) {
        ll y = (ll)(x % (ll)M);
        if (y < 0) y += M;
        val = y;
    }

    constant_modint& operator+=(const constant_modint& rhs) {
        val += rhs.val;
        if (val >= M) val -= M;
        return *this;
    }
    constant_modint& operator-=(const constant_modint& rhs) {
        if (val < rhs.val) val += M;
        val -= rhs.val;
        return *this;
    }
    constant_modint& operator*=(const constant_modint& rhs) {
        val = (ull)val * rhs.val % M;
        return *this;
    }
    constant_modint& operator/=(const constant_modint& rhs) {
        return *this *= rhs.inv();
    }

    constant_modint operator+() const { return *this; }
    constant_modint operator-() const { return constant_modint(0) - *this; }

    friend constant_modint operator+(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) += rhs; }
    friend constant_modint operator-(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) -= rhs; }
    friend constant_modint operator*(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) *= rhs; }
    friend constant_modint operator/(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) /= rhs; }

    friend bool operator==(const constant_modint& lhs, const constant_modint& rhs) { return lhs.val == rhs.val; }
    friend bool operator!=(const constant_modint& lhs, const constant_modint& rhs) { return lhs.val != rhs.val; }

    constant_modint pow(ull n) const {
        constant_modint res = 1, a = *this;
        while (n) {
            if (n&1) res *= a;
            a *= a;
            n >>= 1;
        }
        return res;
    }
    constant_modint inv() const {
        return inverse_mod<ll,ll>(val, M);
    }

    friend std::ostream& operator<<(std::ostream& os, const constant_modint& m) {
        return os << m.val;
    }
    friend std::istream& operator>>(std::istream& is, constant_modint& m) {
        ll x;
        is >> x;
        m = constant_modint(x);
        return is;
    }
};

template <int id = -1>
struct dynamic_modint {
    uint val;
    static uint& mod() {
        static uint M = 998244353;
        return M;
    }
    static void set_mod(uint m) {
        mod() = m;
    }

    dynamic_modint() : val(0) {}
    template <class T> dynamic_modint(T x) {
        ll y = x % (ll)mod();
        if (y < 0) y += mod();
        val = (uint)y;
    }

    dynamic_modint& operator+=(const dynamic_modint& rhs) {
        val += rhs.val;
        if (val >= mod()) val -= mod();
        return *this;
    }
    dynamic_modint& operator-=(const dynamic_modint& rhs) {
        if (val < rhs.val) val += mod();
        val -= rhs.val;
        return *this;
    }
    dynamic_modint& operator*=(const dynamic_modint& rhs) {
        val = (ull)val * rhs.val % mod();
        return *this;
    }
    dynamic_modint& operator/=(const dynamic_modint& rhs) {
        return *this *= rhs.inv();
    }

    dynamic_modint operator+() const { return *this; }
    dynamic_modint operator-() const { return dynamic_modint(0) - *this; }

    friend dynamic_modint operator+(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) += rhs; }
    friend dynamic_modint operator-(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) -= rhs; }
    friend dynamic_modint operator*(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) *= rhs; }
    friend dynamic_modint operator/(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) /= rhs; }

    friend bool operator==(const dynamic_modint& lhs, const dynamic_modint& rhs) { return lhs.val == rhs.val; }
    friend bool operator!=(const dynamic_modint& lhs, const dynamic_modint& rhs) { return lhs.val != rhs.val; }

    dynamic_modint pow(ull n) const {
        dynamic_modint res = 1, a = *this;
        while (n) {
            if (n & 1) res *= a;
            a *= a;
            n >>= 1;
        }
        return res;
    }
    dynamic_modint inv() const {
        ll a = val, b = mod(), u = 1, v = 0;
        while (b) {
            ll t = a / b;
            a -= t * b; swap(a, b);
            u -= t * v; swap(u, v);
        }
        return dynamic_modint(u);
    }

    friend std::ostream& operator<<(std::ostream& os, const dynamic_modint& m) {
        return os << m.val;
    }
    friend std::istream& operator>>(std::istream& is, dynamic_modint& m) {
        ll x;
        is >> x;
        m = dynamic_modint(x);
        return is;
    }
};

#endif /* MODINT__HPP_ */
#ifndef INV_TABLE_HPP_
#define INV_TABLE_HPP_


#include <cassert>

using uint = unsigned;
using ull = unsigned long long;
using ll = long long;
using namespace std;

/// @brief mod M上での階乗,逆元テーブルを保持する構造体
/// @tparam M 
template<uint M>
class mod_table{
    using mint = constant_modint<M>;
    list<mint> invmodlist;
    list<mint> factorialmodlist;
    list<mint> factorialmodinvlist;
    uint N_MAX;

    public:
    constexpr mod_table(const uint N_MAX__){
        N_MAX = max(1u, N_MAX__);
        invmodlist = list<mint>(N_MAX+1);
        factorialmodlist = list<mint>(N_MAX+1);
        factorialmodinvlist = list<mint>(N_MAX+1);
        invmodlist[1] = 1;
        for (uint i = 2; i <= N_MAX; i++){
            invmodlist[i] = (M-M/i)*invmodlist[M%i];
        }
        factorialmodinvlist[0] = 1;
        factorialmodlist[0] = 1;
        for (ull i = 1; i <= N_MAX; i++){
            factorialmodinvlist[i] = invmodlist[i]*factorialmodinvlist[i-1];
            factorialmodlist[i] = factorialmodlist[i-1]*i;
        }
    }
    inline constexpr mint inverse_mod(ll x) const {
        assert(0 <= x and x <= N_MAX);
        return invmodlist[x];
    }
    inline constexpr mint factorialmod(ll x) const {
        assert(0 <= x and x <= N_MAX);
        return factorialmodlist[x];
    }
    inline constexpr mint factorialmodinv(ll x) const {
        assert(0 <= x and x <= N_MAX);
        return factorialmodinvlist[x];
    }
    inline constexpr mint inverse_mod(mint x) const {return inverse_mod(x.val);}
    inline constexpr mint factorialmod(mint x) const {return factorialmod(x.val);}
    inline constexpr mint factorialmodinv(mint x) const {return factorialmodinv(x.val);}
};

#endif /* INV_TABLE_HPP_ */
using namespace std;
using ll = long long;
using uint = unsigned;

#ifndef LEN_DEFINED
#define LEN_DEFINED
template<typename T> inline constexpr ll len(const T& A){return ssize(A);}
#endif


template<uint T>
struct FormalPowerSeries{
    using mint = constant_modint<T>;
    list<mint> f;
    ll sz;
    constexpr FormalPowerSeries(){
        sz = 0;
    }
    /// @brief 指定したサイズのfpsを作成する。
    constexpr FormalPowerSeries(uint _init_sz) : f(_init_sz), sz(_init_sz) {}
    constexpr FormalPowerSeries(uint _init_sz, mint _init_val) : f(_init_sz,_init_val), sz(_init_sz) {}
    template<typename U> constexpr FormalPowerSeries(const list<U>& _init){
        f.assign(_init.begin(), _init.end());
        sz = f.size();
    }
    template<typename U> inline constexpr FormalPowerSeries(initializer_list<U> _init) : FormalPowerSeries(list<U>(_init)){}
    inline constexpr uint size() const {return sz;}
    inline constexpr bool empty() const {return size() == 0;}
    inline constexpr mint* data() {return empty() ? nullptr : &f[0];}
    inline constexpr const mint* data() const {return empty() ? nullptr : &f[0];}
    constexpr void resize(uint _new_size){
        if (sz == _new_size){return;}
        f.resize(_new_size, 0);
        sz = _new_size;
    }
    inline constexpr mint& operator[](uint deg) {return f[deg];}
    inline constexpr const mint& operator[](uint deg) const {return f[deg];}
    template<typename U> constexpr operator list<U>() const {
        list<U> res(sz);
        for (int i = 0; i < sz; i++){
            res[i] = f[i].val;
        }
        return res;
    }

    constexpr void differential(){
        if (sz == 0){return;}
        for (int i = 0; i < sz-1; i++){
            f[i] = mint(i+1)*f[i+1];
        }
        f.pop();
        sz--;
    }
    constexpr void integral(const mod_table<T>& mtable){
        if (sz == 0){return;}
        sz++;
        f.append(0);
        for (int i = sz-1; i > 0; i--){
            f[i] = mtable.inverse_mod(i)*f[i-1];
        }
        f[0] = 0;
    }
    constexpr mint evaluate(const mint& x){
        mint res = 0;
        for (int i = 0; i < sz; i++){
            res *= x;
            res += f[sz-1-i];
        }
        return res;
    }
};

#endif /* FPS__HPP_ */



#ifndef PRIME_AND_DIVISORS_HPP_
#define PRIME_AND_DIVISORS_HPP_

#include <array>


#include <algorithm>
#include <cassert>
#include <numeric>
#include <iostream>
using namespace std;
using ll = long long;
using ulll = __uint128_t;
using pll = array<ll,2>;
using pull = array<ull,2>;
using pii = array<int,2>;



/// @brief 試し割り法で正の整数Nを素因数分解する
/// @return list<array<ll,2>>{{素因数1,個数}, {素因数2,個数}, {素因数3,個数}...}
template<typename T> constexpr list<array<T,2>> trial_division(T N){
    assert(N > 0);
    if (N == 1){
        return list<array<T,2>> {{1,0}};
    }
    list<array<T,2>> R;//戻り値用リスト

    const T M = isqrt(static_cast<ull>(N));
    for (T i = 2; i <= M; i++){
        if (N % i == 0){
            T divide_count = 0;
            while (N % i == 0){
                divide_count++;
                N /= i;
            }
            R.append({i,divide_count});
        }
    }
    if (N != 1){
        R.append({N,1});
    }
    return R;
}

/// @brief 素因数分解リストを受け取って約数関数の値を求める
/// @return 約数のK乗和
template<typename T> constexpr T divisor_function(const list<array<T,2>>& vv, T K, const T MOD = -1){
    if (vv[0][0] == 1){
        return 1;
    }
    T R = 1;
    if (K == 0){
        for (auto x : vv){
            R *= x[1]+1;
            if (MOD > 0){
                R %= MOD;
            }
        }
    }
    else{
        for (auto x : vv){
            ll r = powll(x[0],K);
            R *= (powll(r,x[1]+1) - 1)/(r - 1);
            if (MOD > 0){
                R %= MOD;
            }
        }
    }
    return R;
}

/// @brief 素因数分解の結果pを受け取って、約数リストを生成する。
template<typename T> constexpr list<T> enumerate_divisor(const list<array<T,2>>& p){
    list<T> d{1};
    if (p[0][0] == 1){
        return d;
    }
    for (auto &v : p){
        T t = d.size();
        T temp = 1;
        for (T w = 0; w < v[1]; w++){
            temp *= v[0];
            for (T i = 0; i < t; i++){
                d.append(d[i]*temp);
            }
        }
    }
    sort(vall(d));
    return d;
}

/// @brief 線形篩
/// @attention コンストラクタに整数Nを渡すことでN以下の整数を扱うことができる。
struct LinearSieve{
    list<int> p_list;
    list<int> lpf;
    //Nを渡すことで1以上N以下の整数を扱うことができる
    LinearSieve(int N): lpf(N+1,-1){
        lpf[1] = 1;
        int p_list_size = 0;
        for (int i = 2; i <= N; i++){
            if (lpf[i] < 0){
                p_list.append(i);
                p_list_size++;
                lpf[i] = i;
            }
            for (int j = 0; j < p_list_size && p_list[j] <= lpf[i] && p_list[j]*i <= N; j++){
                lpf[p_list[j]*i] = p_list[j];
            }
        }
    }
    list<pii> factorize(int x){
        if (x == 1){return {{1,0}};}
        list<pii> r;
        do{
            if (r.empty() || lpf[x] != r.back()[0]){
                r.append({lpf[x], 1});
            }
            else{
                r.back()[1]++;
            }
            x /= lpf[x];
        }while(x > 1);
        return r;
    }
    list<list<pii>> factorize_all(int N){
        list<list<pii>> r(N+1);
        r[1].append({1,0});
        for (int i = 2; i <= N; i++){
            r[i] = factorize(i);
        }
        return r;
    }
    /// @brief N以下の整数に対するメビウス関数の値を列挙する。
    list<int> enumerate_mobius(int N){
        list<int> ret(N+1);
        ret[1] = 1;
        for (int i = 2; i <= N; i++){
            if (lpf[i] == i){ret[i] = -1; continue;}
            int temp = i/lpf[i];
            if (lpf[i] == lpf[temp]){
                ret[i] = 0;
            }
            else{
                ret[i] = ret[lpf[i]]*ret[temp];
            }
        }
        return ret;
    }
};

constexpr ull __MillerRabin_small[3] = {2,7,61};
constexpr ull __MillerRabin_large[7] = {2,325,9375,28178,450775,9780504,1795265022};


struct Montgomery64 {
    ull n, ni, r2;
    constexpr Montgomery64(ull n) : n(n), ni(n), r2(-ulll(n) % n){
        for (int i = 0; i < 5; ++i) ni *= 2ull - n * ni;
        ni *= -1;
    }
    constexpr ull reduce(ulll x) const {
        ull m = (ull)x * ni;
        ull res = (x + ulll(m) * n) >> 64;
        return res >= n ? res - n : res;
    }
    constexpr ull mul(ull x, ull y) const {
        return reduce(ulll(x) * y);
    }
    constexpr ull pow_not_reduced(ull a, ull b)const{
        ull res = reduce(r2);
        a = reduce(ulll(a) * r2);
        while (b) {
            if (b & 1) res = mul(res, a);
            a = mul(a, a);
            b >>= 1;
        }
        return res;
    }
};


///@brief ミラーラビン素数判定法
template<typename T> constexpr bool MillerRabin(T N){
    if (N <= 1) return false;
    if (N == 2 || N == 3) return true;
    if (N % 2 == 0) return false;

    T d = N - 1;
    T s = 0;
    while (d % 2 == 0){
        d >>= 1;
        s++;
    }
    auto modpow_T = [&](T a,T b,T m){
        T t = a%m;
        T ans = (m == 1 ? 0 : 1);
        while (b > 0){
            if (b&1){
                ans = (ans*t)%m;
            }
            b >>= 1;
            t = (t*t)%m;
        }
        return ans;
    };
    auto check = [&](T a){
        if (a >= N) a %= N;
        if (a == 0) return true;
        T x = modpow_T(a, d, N);
        if (x == 1 || x == N-1) return true;
        for (int r = 1; r < s; r++){
            x = x*x%N;
            if (x == N - 1) return true;
        }
        return false;
    };

    if (N < 4759123141){
        for (auto a : __MillerRabin_small){
            if (!check(static_cast<T>(a))) return false;
        }
    }
    else if (N < 1000000000000000000){
        for (auto a : __MillerRabin_large){
            if (!check(static_cast<T>(a))) return false;
        }
    }
    else{
        for (auto a : __MillerRabin_large){
            if (!check(static_cast<T>(a))) return false;
        }
        //50個の整数でテスト
        for (T a = 638245987265732141; a < 638245987265732191; a++){
            if (!check(a)) return false;
        }
    }
    return true;
}

/// @brief ミラーラビン素数判定法のull範囲限定版
template<> constexpr bool MillerRabin(ull N){
    if (N <= 1) return false;
    if (N == 2 || N == 3) return true;
    if (N % 2 == 0) return false;

    ull d = N - 1;
    int s = 0;
    while (d % 2 == 0){
        d >>= 1;
        s++;
    }
    Montgomery64 mg(N);
    auto check = [&](ull a){
        if (a >= N) a %= N;
        if (a == 0) return true;
        ull x = mg.pow_not_reduced(a, d);
        ull y = mg.reduce(x);
        if (y == 1 || y == N-1) return true;
        for (int r = 1; r < s; r++){
            x = mg.mul(x, x);
            if (mg.reduce(x) == N - 1) return true;
        }
        return false;
    };

    if (N < 4759123141ull){
        for (auto a : __MillerRabin_small){
            if (!check(a)) return false;
        }
    }
    else{
        for (auto a : __MillerRabin_large){
            if (!check(a)) return false;
        }
    }
    return true;
}



/// @brief ポラード・ロー法による素因数分解
template<typename T> constexpr list<array<T, 2>> Pollard_rho(T N){
    assert(N > 0);
    if (N == 1){return {{1,0}};}
    list<T> res;
    list<T> dq;
    dq.append(N);
    while (!dq.empty()){
        T n = dq.back();
        dq.pop();
        bool next_loop = false;
        while (!MillerRabin<T>(n)){
            if ((n&1) == 0){
                while ((n&1) == 0){res.append(2); n>>=1;}
                if (n == 1){next_loop = true; break;}
                continue;
            }
            T blocksize = sqrt(sqrt(sqrt(n)));
            for (T c = 1; c <= n; c++){
                T y = (0b1010110101101011110011ull^c^n)%n;
                T r = 1;
                T k = 0;
                T k_mod_blocksize = 0;
                bool next_loop_inner = false;
                for (int _ = 0;;){
                    T q = 1;
                    T q_old = q;
                    T y_old = y;
                    T x = y;
                    while (k < r){
                        k++;
                        k_mod_blocksize++;
                        y = y*y%n+c; if (y >= n){ y -= n;}
                        q = q*(y >= x ? (y-x) : (x-y))%n;
                        if (k_mod_blocksize == blocksize){
                            k_mod_blocksize -= blocksize;
                            T g = gcd(q,n);
                            if (g > 1){
                                if (g != n){
                                    dq.append(g);
                                    dq.append(n/g);
                                    next_loop = true; break;
                                }
                                else{
                                    y = y_old;
                                    q = q_old;
                                    for (int i = 0; i < blocksize; i++){
                                        y = y*y%n+c; if (y >= n){ y -= n;}
                                        q = q*(y >= x ? (y-x) : (x-y))%n;
                                        g = gcd(q,n);
                                        if (g > 1){
                                            if (g != n){
                                                dq.append(g);
                                                dq.append(n/g);
                                                next_loop = true; break;
                                            }
                                            next_loop_inner = true; break;
                                        }
                                    }
                                    if (next_loop || next_loop_inner) break;
                                }
                            }
                            q_old = q;
                            y_old = y;
                        }
                    }
                    if (next_loop || next_loop_inner) break;
                    q_old = gcd(q,n);
                    if (q_old > 1){
                        if (q_old != n){
                            dq.append(q_old);
                            dq.append(n/q_old);
                            next_loop = true; break;
                        }
                        next_loop_inner = true; break;
                    }
                    r <<= 1;
                }
                if (next_loop) break;
            }
            if (next_loop) break;
        }
        if (!next_loop) {
            res.append(n);
        }
    }
    sort(res.begin(), res.end());
    list<array<T,2>> res2;
    for (ull l = 0, r = 0, sz = res.size(); l < sz;){
        while (r < sz && res[l] == res[r]){
            r++;
        }
        res2.append({res[l], r-l});
        l = r;
    }
    return res2;
}

/// @brief ポラード・ロー法による素因数分解のull範囲限定版
template<> constexpr list<pull> Pollard_rho(ull N){
    assert(N > 0);
    if (N == 1){return {{1,0}};}
    list<ull> res;
    list<ull> dq;
    dq.append(N);
    while (!dq.empty()){
        ull n = dq.back();
        dq.pop();
        bool next_loop = false;
        while (!MillerRabin<ull>(n)){
            if ((n&1) == 0){
                while ((n&1) == 0){res.append(2); n>>=1;}
                if (n == 1){next_loop = true; break;}
                continue;
            }
            Montgomery64 mg(n);
            int blocksize = max<int>(1, pow(n, 0.125));
            for (ull c = 1; c <= n; c++){
                ull reduced_c = mg.reduce(c*(ulll)mg.r2);
                ull y = (0b1010110101101011110011ull^c^n)%n; y = mg.reduce(y*(ulll)mg.r2);
                ull r = 1;
                ull k = 0;
                int k_mod_blocksize = 0;
                bool next_loop_inner = false;
                for (int _ = 0; _ < 63; _++){
                    ull q = mg.reduce(mg.r2);
                    ull q_old = q;
                    ull y_old = y;
                    ull x = y;
                    while (k < r){
                        k++;
                        k_mod_blocksize++;
                        y = mg.reduce(y*(ulll)y)+reduced_c; if (y >= n){ y -= n;}
                        q = mg.reduce(q*(y >= x ? (ulll)(y-x) : (ulll)(x-y)));
                        if (k_mod_blocksize == blocksize){
                            k_mod_blocksize -= blocksize;
                            ull g = gcd(mg.reduce(q),n);
                            if (g > 1){
                                if (g != n){
                                    dq.append(g);
                                    dq.append(n/g);
                                    next_loop = true; break;
                                }
                                else{
                                    y = y_old;
                                    q = q_old;
                                    for (int i = 0; i < blocksize; i++){
                                        y = mg.reduce(y*(ulll)y)+reduced_c; if (y >= n){ y -= n;}
                                        q = mg.reduce(q*(y >= x ? (ulll)(y-x) : (ulll)(x-y)));
                                        g = gcd(mg.reduce(q),n);
                                        if (g > 1){
                                            if (g != n){
                                                dq.append(g);
                                                dq.append(n/g);
                                                next_loop = true; break;
                                            }
                                            next_loop_inner = true; break;
                                        }
                                    }
                                    if (next_loop || next_loop_inner) break;
                                }
                            }
                            q_old = q;
                            y_old = y;
                        }
                    }
                    if (next_loop || next_loop_inner) break;
                    q_old = gcd(mg.reduce(q),n);
                    if (q_old > 1){
                        if (q_old != n){
                            dq.append(q_old);
                            dq.append(n/q_old);
                            next_loop = true; break;
                        }
                        next_loop_inner = true; break;
                    }
                    r <<= 1;
                }
                if (next_loop) break;
            }
            if (next_loop) break;
        }
        if (!next_loop) {
            res.append(n);
        }
    }
    sort(res.begin(), res.end());
    list<pull> res2;
    for (int l = 0, r = 0, sz = res.size(); l < sz;){
        while (r < sz && res[l] == res[r]){
            r++;
        }
        res2.append({res[l], static_cast<ull>(r-l)});
        l = r;
    }
    return res2;
}


/// @brief 素因数分解を行う。`N`の大きさによって、試し割り法とポラード・ロー法を使い分けてくれる。
template<typename T> constexpr list<array<T, 2>> factorize(T N){
    if (N < 640000){
        return trial_division<T>(N);
    }
    if (N < 1000000000000000000){
        return Pollard_rho<T>(N);
    }
    return Pollard_rho<T>(N);
}
template<> constexpr list<pull> factorize(ull N){
    if (N < 640000ull){
        return trial_division<ull>(N);
    }
    return Pollard_rho<ull>(N);
}
template<> constexpr list<pll> factorize(ll N){
    auto res = factorize<ull>(N);
    list<pll> res2(res.size());
    for (int i = 0, sz = res.size(); i < sz; i++){
        res2[i][0] = res[i][0];
        res2[i][1] = res[i][1];
    }
    return res2;
}
template<> constexpr list<pii> factorize(int N){
    auto res = factorize<ull>(N);
    list<pii> res2(res.size());
    for (int i = 0, sz = res.size(); i < sz; i++){
        res2[i][0] = res[i][0];
        res2[i][1] = res[i][1];
    }
    return res2;
}

struct Xorshift64{
    ull state;
    constexpr Xorshift64(ull seed = 881726454633252252ull) : state(seed) {}
    constexpr ull next() {
        state ^= state << 13;
        state ^= state >> 7;
        state ^= state << 17;
        return state;
    }
};

/// @brief mod Pにおける原始根を1つ探す。Pはull範囲内の素数でなければならない。
constexpr ull primitive_root(ull P){
    assert(MillerRabin<ull>(P));
    if (P == 2){
        return 1;
    }
    Montgomery64 mg(P);
    Xorshift64 rng(0x1234567E89ABCDEF);
    auto factorized = factorize<ull>(P-1);
    while (true) {
        ull r = rng.next()%P;
        if (r == 0){continue;}
        bool ok = true;
        for (auto& q : factorized){
            if (mg.reduce(mg.pow_not_reduced(r, (P-1)/q[0])) == 1){
                ok = false;
                break;
            }
        }
        if (ok) return r;
    }
}

#endif /* PRIME_AND_DIVISORS_HPP_ */




using namespace std;
using ll = long long;
using uint = unsigned;
using ull = unsigned long long;
#define vall(A) A.begin(), A.end()

template<typename T> inline bool btest_for_fps(T K, int i){return K&(1ull<<i);}


/// @brief mod M上での形式的冪級数の計算を行う構造体
template<uint M>
struct fps_operator{
    template<uint T> using fps = FormalPowerSeries<T>;
    template<uint T> using mint = constant_modint<T>;
    uint sum_e[30];
    uint sum_ie[30];
    uint log_max_length;
    uint last_powroot;

    // --- モンゴメリ乗算用の定数と関数 ---
    static constexpr uint get_r() {
        uint res = M;
        for (int i = 0; i < 4; ++i) res *= 2 - M * res;
        return res;
    }
    static constexpr uint R = get_r();
    static constexpr uint R2 = -ull(M) % M;
    static constexpr uint NEG_INV = 0 - R;

    static inline constexpr uint reduce(ull x) {
        uint res = (x + ull(uint(x) * NEG_INV) * M) >> 32;
        return res >= M ? res - M : res;
    }
    static inline constexpr uint to_montgomery(uint x) {
        return reduce(ull(x) * R2);
    }
    static inline constexpr uint from_montgomery(uint x) {
        return reduce(x);
    }
    static inline constexpr uint montgomery_mul(uint x, uint y) {
        return reduce(ull(x) * y);
    }
    static inline constexpr uint mod_add(uint x, uint y) {
        return x + y >= M ? x + y - M : x + y;
    }
    static inline constexpr uint mod_sub(uint x, uint y) {
        return x < y ? x + M - y : x - y;
    }

    // --- AVX-512 / AVX2 SIMD Intrinsics ---
#if defined(__AVX512F__) && defined(__AVX512DQ__)
    static inline __m512i montgomery_mul_simd(__m512i a, __m512i b, __m512i mod, __m512i neg_inv) {
        __m512i mul0 = _mm512_mul_epu32(a, b);
        __m512i mul1 = _mm512_mul_epu32(_mm512_srli_epi64(a, 32), _mm512_srli_epi64(b, 32));
        __m512i q0 = _mm512_mul_epu32(_mm512_mullo_epi32(mul0, neg_inv), mod);
        __m512i q1 = _mm512_mul_epu32(_mm512_mullo_epi32(mul1, neg_inv), mod);
        __m512i res0 = _mm512_srli_epi64(_mm512_add_epi64(mul0, q0), 32);
        __m512i res1 = _mm512_srli_epi64(_mm512_add_epi64(mul1, q1), 32);
        __m512i res = _mm512_or_si512(res0, _mm512_slli_epi64(res1, 32));
        __mmask16 cmp = _mm512_cmpge_epu32_mask(res, mod);
        return _mm512_mask_sub_epi32(res, cmp, res, mod);
    }
    static inline __m512i mod_add_simd(__m512i a, __m512i b, __m512i mod) {
        __m512i res = _mm512_add_epi32(a, b);
        __mmask16 cmp = _mm512_cmpge_epu32_mask(res, mod);
        return _mm512_mask_sub_epi32(res, cmp, res, mod);
    }
    static inline __m512i mod_sub_simd(__m512i a, __m512i b, __m512i mod) {
        __mmask16 cmp = _mm512_cmplt_epu32_mask(a, b);
        __m512i res = _mm512_sub_epi32(a, b);
        return _mm512_mask_add_epi32(res, cmp, res, mod);
    }
#elif defined(__AVX2__)
    // AVX2用 モンゴメリ乗算・剰余加減算
    static inline __m256i montgomery_mul_simd(__m256i a, __m256i b, __m256i mod, __m256i neg_inv) {
        __m256i mul0 = _mm256_mul_epu32(a, b);
        __m256i mul1 = _mm256_mul_epu32(_mm256_srli_epi64(a, 32), _mm256_srli_epi64(b, 32));
        __m256i q0 = _mm256_mul_epu32(_mm256_mullo_epi32(mul0, neg_inv), mod);
        __m256i q1 = _mm256_mul_epu32(_mm256_mullo_epi32(mul1, neg_inv), mod);
        __m256i res0 = _mm256_srli_epi64(_mm256_add_epi64(mul0, q0), 32);
        __m256i res1 = _mm256_srli_epi64(_mm256_add_epi64(mul1, q1), 32);
        
        __m256i res = _mm256_blend_epi32(res0, _mm256_slli_epi64(res1, 32), 0xAA);
        __m256i diff = _mm256_sub_epi32(res, mod);
        // diff < 0 なら -1 (全ビット1)、そうでないなら 0 となるマスク
        __m256i mask = _mm256_srai_epi32(diff, 31);
        return _mm256_add_epi32(diff, _mm256_and_si256(mask, mod));
    }
    static inline __m256i mod_add_simd(__m256i a, __m256i b, __m256i mod) {
        __m256i diff = _mm256_sub_epi32(_mm256_add_epi32(a, b), mod);
        __m256i mask = _mm256_srai_epi32(diff, 31);
        return _mm256_add_epi32(diff, _mm256_and_si256(mask, mod));
    }
    static inline __m256i mod_sub_simd(__m256i a, __m256i b, __m256i mod) {
        __m256i diff = _mm256_sub_epi32(a, b);
        __m256i mask = _mm256_srai_epi32(diff, 31);
        return _mm256_add_epi32(diff, _mm256_and_si256(mask, mod));
    }
#endif


    constexpr fps_operator(){
        list<ll> powroot{1};
        list<ll> powrootinv;
        while (powroot.back() >= 0){
            powroot.append(cipolla<M>(powroot.back()));
            if (powroot.back() == powroot[powroot.size()-2]){
                powroot.back() = M-powroot.back();
            }
        }
        powroot.pop();
        log_max_length = powroot.size()-1;
        last_powroot = powroot.back();
        for (auto v : powroot){
            powrootinv.append(inverse_mod<ll,ll>(v,M));
        }
        int cnt2 = powroot.size()-1;
        uint now = reduce(R2);
        for (int i = 0; i <= cnt2-2; i++){
            sum_e[i] = montgomery_mul(to_montgomery(powroot[i+2]), now);
            now = montgomery_mul(now, to_montgomery(powrootinv[i+2]));
        }
        uint inow = reduce(R2);
        for (int i = 0; i <= cnt2-2; i++){
            sum_ie[i] = montgomery_mul(to_montgomery(powrootinv[i+2]), inow);
            inow = montgomery_mul(inow, to_montgomery(powroot[i+2]));
        }
        for (int i = cnt2-1; i < 30; i++){
            sum_e[i] = 0;
            sum_ie[i] = 0;
        }
    }

    // Public transforms return canonical residues in the original bit-reversed
    // order. IDFT remains unnormalized. Transposes share the same kernels.
private:
    template<bool InverseButterfly, bool InverseRoots>
    void transform_scalar(mint<M>* a, int n) const {
        const uint* rate = InverseRoots ? sum_ie : sum_e;
        const int h = __builtin_ctz(unsigned(n));
        for (int stage = 0; stage < h; ++stage) {
            const int p = InverseButterfly ? 1 << stage : n >> (stage + 1);
            uint w = to_montgomery(1);
            for (int offset = 0, s = 0; offset < n; offset += 2*p, ++s) {
                for (int i = 0; i < p; ++i) {
                    const uint l = a[offset+i].val, r = a[offset+i+p].val;
                    if constexpr (InverseButterfly) {
                        a[offset+i].val = mod_add(l, r);
                        a[offset+i+p].val = montgomery_mul(mod_sub(l, r), w);
                    } else {
                        const uint t = montgomery_mul(r, w);
                        a[offset+i].val = mod_add(l, t);
                        a[offset+i+p].val = mod_sub(l, t);
                    }
                }
                w = montgomery_mul(w, rate[__builtin_ctz(~unsigned(s))]);
            }
        }
    }

#if defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
    // Internal SIMD residues are in [0, 2M). M < 2^30 guarantees that
    // additions fit uint32_t and Montgomery products reduce to [0, 2M).
    __attribute__((target("avx2")))
    static inline __m256i mul_lazy256(__m256i a, __m256i b) {
        const __m256i mod = _mm256_set1_epi32(M);
        const __m256i ni = _mm256_set1_epi32(NEG_INV);
        const __m256i x = _mm256_mul_epu32(a, b);
        const __m256i y = _mm256_mul_epu32(_mm256_srli_epi64(a, 32), _mm256_srli_epi64(b, 32));
        // Only the low 32 bits of each 64-bit lane are consumed below; a
        // widening multiply avoids computing the unused odd 32-bit products.
        const __m256i qx = _mm256_mul_epu32(_mm256_mul_epu32(x, ni), mod);
        const __m256i qy = _mm256_mul_epu32(_mm256_mul_epu32(y, ni), mod);
        const __m256i rx = _mm256_srli_epi64(_mm256_add_epi64(x, qx), 32);
        const __m256i ry = _mm256_add_epi64(y, qy);
        const __m256i r = _mm256_blend_epi32(rx, ry, 0xaa);
        return r;
    }
    __attribute__((target("avx2")))
    static inline __m256i mul256(__m256i a, __m256i b) {
        const __m256i r = mul_lazy256(a,b);
        return _mm256_min_epu32(r,_mm256_sub_epi32(r,_mm256_set1_epi32(M)));
    }
    __attribute__((target("avx2")))
    static inline __m256i add_lazy256(__m256i a, __m256i b) {
        const __m256i r = _mm256_add_epi32(a, b);
        return _mm256_min_epu32(r, _mm256_sub_epi32(r, _mm256_set1_epi32(2*M)));
    }
    __attribute__((target("avx2")))
    static inline __m256i sub_lazy256(__m256i a, __m256i b) {
        const __m256i r = _mm256_sub_epi32(a, b);
        return _mm256_min_epu32(r, _mm256_add_epi32(r, _mm256_set1_epi32(2*M)));
    }

    __attribute__((target("avx2")))
    static inline __m256i shoup256(__m256i a, __m256i b, __m256i bq) {
        // b is an ordinary residue; bq = floor(b*2^32/M). The result is
        // congruent to a*b and lies in [0, 2M), without a final correction.
        const __m256i q0 = _mm256_srli_epi64(_mm256_mul_epu32(a,bq),32);
        const __m256i q1 = _mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(bq,32));
        const __m256i q = _mm256_blend_epi32(q0,q1,0xaa);
        const __m256i mod = _mm256_set1_epi32(M);
        const __m256i r = _mm256_sub_epi32(_mm256_mullo_epi32(a,b),_mm256_mullo_epi32(q,mod));
        return r;
    }

    template<bool InverseButterfly, bool InverseRoots, bool Canonical = true, bool ZeroUpperHalf = false>
    __attribute__((target("avx2")))
    void transform_avx2(mint<M>* a, int n) const {
        const uint* rate = InverseRoots ? sum_ie : sum_e;
        const int h = __builtin_ctz(unsigned(n));
        uint rate4[29];
        for (int i = 0; i+1 < h; ++i) rate4[i] = montgomery_mul(rate[0],rate[i+1]);
        const __m256i imag = _mm256_set1_epi32(rate[0]);
        const uint im = from_montgomery(rate[0]);
        const __m256i imv = _mm256_set1_epi32(im);
        const __m256i imq = _mm256_set1_epi32((ull(im)<<32)/M);
        int first_stage = 0;
        if constexpr (ZeroUpperHalf && !InverseButterfly) {
            // The first radix-4 block has no twiddle and its upper two inputs
            // are zero. Do not load padding or execute the redundant butterflies.
            const int p = n/4;
            for (int i = 0; i < p; i += 8) {
                const __m256i x = _mm256_loadu_si256((const __m256i*)(a+i));
                const __m256i y = _mm256_loadu_si256((const __m256i*)(a+i+p));
                const __m256i t = shoup256(y,imv,imq);
                _mm256_storeu_si256((__m256i*)(a+i),add_lazy256(x,y));
                _mm256_storeu_si256((__m256i*)(a+i+p),sub_lazy256(x,y));
                _mm256_storeu_si256((__m256i*)(a+i+2*p),add_lazy256(x,t));
                _mm256_storeu_si256((__m256i*)(a+i+3*p),sub_lazy256(x,t));
            }
            first_stage = 2;
        }
        for (int stage = first_stage; stage < h; ++stage) {
            const int p4 = InverseButterfly ? 1 << stage : n >> (stage + 2);
            if (stage + 1 < h && p4 >= 8) {
                uint w = to_montgomery(1);
                for (int offset = 0, s = 0; offset < n; offset += 4*p4, ++s) {
                    const uint w2 = montgomery_mul(w,w), w3 = montgomery_mul(w2,w);
                    const uint b1 = from_montgomery(w), b2 = from_montgomery(w2), b3 = from_montgomery(w3);
                    const __m256i v1 = _mm256_set1_epi32(b1), q1 = _mm256_set1_epi32((ull(b1)<<32)/M);
                    const __m256i v2 = _mm256_set1_epi32(b2), q2 = _mm256_set1_epi32((ull(b2)<<32)/M);
                    const __m256i v3 = _mm256_set1_epi32(b3), q3 = _mm256_set1_epi32((ull(b3)<<32)/M);
                    for (int i = 0; i < p4; i += 8) {
                        __m256i a0 = _mm256_loadu_si256((const __m256i*)(a+offset+i));
                        __m256i a1 = _mm256_loadu_si256((const __m256i*)(a+offset+i+p4));
                        __m256i a2 = _mm256_loadu_si256((const __m256i*)(a+offset+i+2*p4));
                        __m256i a3 = _mm256_loadu_si256((const __m256i*)(a+offset+i+3*p4));
                        __m256i x,y,z,t;
                        if constexpr (InverseButterfly) {
                            x = add_lazy256(a0,a1); y = sub_lazy256(a0,a1);
                            z = add_lazy256(a2,a3); t = shoup256(sub_lazy256(a2,a3),imv,imq);
                            a0 = add_lazy256(x,z); a1 = add_lazy256(y,t);
                            a2 = sub_lazy256(x,z); a3 = sub_lazy256(y,t);
                            if (s != 0) {
                                a1 = shoup256(a1,v1,q1); a2 = shoup256(a2,v2,q2); a3 = shoup256(a3,v3,q3);
                            }
                        } else {
                            if (s != 0) {
                                a1 = shoup256(a1,v1,q1); a2 = shoup256(a2,v2,q2); a3 = shoup256(a3,v3,q3);
                            }
                            x = add_lazy256(a0,a2); y = add_lazy256(a1,a3);
                            z = sub_lazy256(a0,a2); t = shoup256(sub_lazy256(a1,a3),imv,imq);
                            a0 = add_lazy256(x,y); a1 = sub_lazy256(x,y);
                            a2 = add_lazy256(z,t); a3 = sub_lazy256(z,t);
                        }
                        _mm256_storeu_si256((__m256i*)(a+offset+i),a0);
                        _mm256_storeu_si256((__m256i*)(a+offset+i+p4),a1);
                        _mm256_storeu_si256((__m256i*)(a+offset+i+2*p4),a2);
                        _mm256_storeu_si256((__m256i*)(a+offset+i+3*p4),a3);
                    }
                    w = montgomery_mul(w,rate4[__builtin_ctz(~unsigned(s))]);
                }
                ++stage;
                continue;
            }
            if (stage + 1 < h && n >= 32) {
                const int p = p4, groups = 8/p;
                alignas(32) uint lane1[8], lane2[8], lane3[8];
                uint w = to_montgomery(1), last = w;
                for (int s = 0; s < groups; ++s) {
                    const uint w2 = montgomery_mul(w,w), w3 = montgomery_mul(w2,w);
                    for (int i = 0; i < p; ++i) {
                        const int j = p == 1 ? (s/2 + (s%2)*4) : s*p+i;
                        lane1[j] = w; lane2[j] = w2; lane3[j] = w3;
                    }
                    last = w;
                    w = montgomery_mul(w,rate4[__builtin_ctz(~unsigned(s))]);
                }
                __m256i v1 = _mm256_load_si256((const __m256i*)lane1);
                __m256i v2 = _mm256_load_si256((const __m256i*)lane2);
                __m256i v3 = _mm256_load_si256((const __m256i*)lane3);
                const int shift = __builtin_ctz(unsigned(groups));
                uint advance1[30], advance2[30], advance3[30];
                for (int i = 0; i+shift+1 < h; ++i) {
                    const uint d = montgomery_mul(last,rate4[i+shift]);
                    advance1[i] = d;
                    advance2[i] = montgomery_mul(d,d);
                    advance3[i] = montgomery_mul(advance2[i],d);
                }
                for (int offset = 0, s = 0; offset < n; offset += 32, ++s) {
                    __m256i x0 = _mm256_loadu_si256((const __m256i*)(a+offset));
                    __m256i x1 = _mm256_loadu_si256((const __m256i*)(a+offset+8));
                    __m256i x2 = _mm256_loadu_si256((const __m256i*)(a+offset+16));
                    __m256i x3 = _mm256_loadu_si256((const __m256i*)(a+offset+24));
                    __m256i a0,a1,a2,a3;
                    if (p == 1) {
                        const __m256i t0 = _mm256_unpacklo_epi32(x0,x1), t1 = _mm256_unpackhi_epi32(x0,x1);
                        const __m256i t2 = _mm256_unpacklo_epi32(x2,x3), t3 = _mm256_unpackhi_epi32(x2,x3);
                        a0 = _mm256_unpacklo_epi64(t0,t2); a1 = _mm256_unpackhi_epi64(t0,t2);
                        a2 = _mm256_unpacklo_epi64(t1,t3); a3 = _mm256_unpackhi_epi64(t1,t3);
                    } else if (p == 2) {
                        const __m256i t0 = _mm256_unpacklo_epi64(x0,x1), t1 = _mm256_unpackhi_epi64(x0,x1);
                        const __m256i t2 = _mm256_unpacklo_epi64(x2,x3), t3 = _mm256_unpackhi_epi64(x2,x3);
                        a0 = _mm256_permute2x128_si256(t0,t2,0x20); a1 = _mm256_permute2x128_si256(t1,t3,0x20);
                        a2 = _mm256_permute2x128_si256(t0,t2,0x31); a3 = _mm256_permute2x128_si256(t1,t3,0x31);
                    } else {
                        a0 = _mm256_permute2x128_si256(x0,x2,0x20); a1 = _mm256_permute2x128_si256(x0,x2,0x31);
                        a2 = _mm256_permute2x128_si256(x1,x3,0x20); a3 = _mm256_permute2x128_si256(x1,x3,0x31);
                    }
                    __m256i x,y,z,t;
                    if constexpr (InverseButterfly) {
                        x = add_lazy256(a0,a1); y = sub_lazy256(a0,a1);
                        z = add_lazy256(a2,a3); t = mul_lazy256(sub_lazy256(a2,a3),imag);
                        a0 = add_lazy256(x,z); a1 = mul_lazy256(add_lazy256(y,t),v1);
                        a2 = mul_lazy256(sub_lazy256(x,z),v2); a3 = mul_lazy256(sub_lazy256(y,t),v3);
                    } else {
                        a1 = mul_lazy256(a1,v1); a2 = mul_lazy256(a2,v2); a3 = mul_lazy256(a3,v3);
                        x = add_lazy256(a0,a2); y = add_lazy256(a1,a3);
                        z = sub_lazy256(a0,a2); t = mul_lazy256(sub_lazy256(a1,a3),imag);
                        a0 = add_lazy256(x,y); a1 = sub_lazy256(x,y);
                        a2 = add_lazy256(z,t); a3 = sub_lazy256(z,t);
                    }
                    if (p == 1) {
                        const __m256i t0 = _mm256_unpacklo_epi32(a0,a1), t1 = _mm256_unpackhi_epi32(a0,a1);
                        const __m256i t2 = _mm256_unpacklo_epi32(a2,a3), t3 = _mm256_unpackhi_epi32(a2,a3);
                        x0 = _mm256_unpacklo_epi64(t0,t2); x1 = _mm256_unpackhi_epi64(t0,t2);
                        x2 = _mm256_unpacklo_epi64(t1,t3); x3 = _mm256_unpackhi_epi64(t1,t3);
                    } else if (p == 2) {
                        const __m256i t0 = _mm256_permute2x128_si256(a0,a2,0x20), t1 = _mm256_permute2x128_si256(a1,a3,0x20);
                        const __m256i t2 = _mm256_permute2x128_si256(a0,a2,0x31), t3 = _mm256_permute2x128_si256(a1,a3,0x31);
                        x0 = _mm256_unpacklo_epi64(t0,t1); x1 = _mm256_unpackhi_epi64(t0,t1);
                        x2 = _mm256_unpacklo_epi64(t2,t3); x3 = _mm256_unpackhi_epi64(t2,t3);
                    } else {
                        x0 = _mm256_permute2x128_si256(a0,a1,0x20); x1 = _mm256_permute2x128_si256(a2,a3,0x20);
                        x2 = _mm256_permute2x128_si256(a0,a1,0x31); x3 = _mm256_permute2x128_si256(a2,a3,0x31);
                    }
                    _mm256_storeu_si256((__m256i*)(a+offset),x0);
                    _mm256_storeu_si256((__m256i*)(a+offset+8),x1);
                    _mm256_storeu_si256((__m256i*)(a+offset+16),x2);
                    _mm256_storeu_si256((__m256i*)(a+offset+24),x3);
                    const int k = __builtin_ctz(~unsigned(s));
                    v1 = mul256(v1,_mm256_set1_epi32(advance1[k]));
                    v2 = mul256(v2,_mm256_set1_epi32(advance2[k]));
                    v3 = mul256(v3,_mm256_set1_epi32(advance3[k]));
                }
                ++stage;
                continue;
            }
            const int p = InverseButterfly ? 1 << stage : n >> (stage + 1);
            uint w = to_montgomery(1);
            if (p >= 8) {
                for (int offset = 0, s = 0; offset < n; offset += 2*p, ++s) {
                    const __m256i vw = _mm256_set1_epi32(w);
                    for (int i = 0; i < p; i += 8) {
                        __m256i l = _mm256_loadu_si256((const __m256i*)(a+offset+i));
                        __m256i r = _mm256_loadu_si256((const __m256i*)(a+offset+i+p));
                        __m256i x, y;
                        if constexpr (InverseButterfly) {
                            x = add_lazy256(l, r); y = mul_lazy256(sub_lazy256(l, r), vw);
                        } else {
                            r = mul_lazy256(r, vw); x = add_lazy256(l, r); y = sub_lazy256(l, r);
                        }
                        _mm256_storeu_si256((__m256i*)(a+offset+i), x);
                        _mm256_storeu_si256((__m256i*)(a+offset+i+p), y);
                    }
                    w = montgomery_mul(w, rate[__builtin_ctz(~unsigned(s))]);
                }
            } else {
                // Pack 8 butterflies even when each block has only 1, 2 or 4 pairs.
                alignas(32) uint lane[8]{};
                uint t = to_montgomery(1), last = t;
                const int groups = 8/p;
                for (int s = 0; s < groups; ++s) {
                    for (int i = 0; i < p; ++i) lane[s*p+i] = t;
                    last = t;
                    t = montgomery_mul(t, rate[__builtin_ctz(~unsigned(s))]);
                }
                const __m256i base = _mm256_load_si256((const __m256i*)lane);
                const __m256i unpack = p == 1 ? _mm256_setr_epi32(0,2,4,6,1,3,5,7)
                                             : _mm256_setr_epi32(0,1,4,5,2,3,6,7);
                const __m256i pack = p == 1 ? _mm256_setr_epi32(0,4,1,5,2,6,3,7) : unpack;
                const int shift = __builtin_ctz(unsigned(groups));
                uint advance[30];
                for (int i = 0; i + shift < h; ++i) advance[i] = montgomery_mul(last, rate[i+shift]);
                for (int offset = 0, s = 0; offset < n; offset += 16, ++s) {
                    __m256i x = _mm256_loadu_si256((const __m256i*)(a+offset));
                    __m256i y = _mm256_loadu_si256((const __m256i*)(a+offset+8));
                    if (p != 4) {
                        x = _mm256_permutevar8x32_epi32(x, unpack);
                        y = _mm256_permutevar8x32_epi32(y, unpack);
                    }
                    __m256i l = _mm256_permute2x128_si256(x,y,0x20);
                    __m256i r = _mm256_permute2x128_si256(x,y,0x31);
                    const __m256i vw = mul256(base, _mm256_set1_epi32(w));
                    if constexpr (InverseButterfly) {
                        x = add_lazy256(l,r); y = mul_lazy256(sub_lazy256(l,r),vw);
                    } else {
                        r = mul_lazy256(r,vw); x = add_lazy256(l,r); y = sub_lazy256(l,r);
                    }
                    l = _mm256_permute2x128_si256(x,y,0x20);
                    r = _mm256_permute2x128_si256(x,y,0x31);
                    if (p != 4) {
                        l = _mm256_permutevar8x32_epi32(l,pack);
                        r = _mm256_permutevar8x32_epi32(r,pack);
                    }
                    _mm256_storeu_si256((__m256i*)(a+offset),l);
                    _mm256_storeu_si256((__m256i*)(a+offset+8),r);
                    w = montgomery_mul(w,advance[__builtin_ctz(~unsigned(s))]);
                }
            }
        }
        const __m256i mod = _mm256_set1_epi32(M);
        if constexpr (Canonical) for (int i = 0; i < n; i += 8) {
            __m256i x = _mm256_loadu_si256((const __m256i*)(a+i));
            x = _mm256_min_epu32(x,_mm256_sub_epi32(x,mod));
            _mm256_storeu_si256((__m256i*)(a+i),x);
        }
    }

    template<bool Canonical = true>
    __attribute__((target("avx2")))
    static void multiply_avx2(mint<M>* dst, const mint<M>* a, const mint<M>* b, int n, uint scale) {
        const __m256i vs = _mm256_set1_epi32(scale);
        for (int i = 0; i < n; i += 8) {
            __m256i x = _mm256_loadu_si256((const __m256i*)(a+i));
            const __m256i y = _mm256_loadu_si256((const __m256i*)(b+i));
            x = mul_lazy256(mul_lazy256(x,y),vs);
            if constexpr (Canonical)
                x = _mm256_min_epu32(x,_mm256_sub_epi32(x,_mm256_set1_epi32(M)));
            _mm256_storeu_si256((__m256i*)(dst+i),x);
        }
    }
#endif

    static bool use_avx2() {
#if defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
        // This check also makes ordinary -O2/-O3 builds use SIMD on supported CPUs.
#ifdef FPS_DISABLE_SIMD
        return false;
#else
        return M < (1u << 30) && __builtin_cpu_supports("avx2");
#endif
#else
        return false;
#endif
    }

    // Pointwise product with optional normalization, all operands in ordinary
    // representation. dst may alias either input.
    void pointwise_product(fps<M>& dst, const fps<M>& lhs, const fps<M>& rhs,
                           int n, mint<M> factor = 1) const {
        if (n == 0) return;
        auto* out = &dst[0];
        const auto* a = &lhs[0];
        const auto* b = &rhs[0];
        int i = 0;
#if defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
        if (use_avx2()) {
            i = n & ~7;
            multiply_avx2(out,a,b,i,ull(R2)*factor.val%M);
        }
#endif
        for (; i < n; ++i) out[i] = a[i]*b[i]*factor;
    }

    template<bool InverseButterfly, bool InverseRoots>
    void transform(fps<M>& F) const {
        if (F.empty()) return;
        F.f.reserve(upperpow2(F.sz));
        F.resize(upperpow2(F.sz));
        const int n = F.size();
        assert(unsigned(n) <= (1u << log_max_length));
        auto* a = &F[0];
#if defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
        if (n >= 16 && use_avx2()) {
            transform_avx2<InverseButterfly, InverseRoots>(a,n);
            return;
        }
#endif
        transform_scalar<InverseButterfly, InverseRoots>(a,n);
    }
public:
    void inplaceDFT(fps<M>& F) const { transform<false,false>(F); }
    void inplaceIDFT(fps<M>& F) const { transform<true,true>(F); }
    void inplaceDFT_T(fps<M>& F) const { transform<true,false>(F); }
    void inplaceIDFT_T(fps<M>& F) const { transform<false,true>(F); }

private:
    static fps<M> convolution_copy(const fps<M>& src, int capacity) {
        assert(capacity > 0 && src.size() <= unsigned(capacity));
        fps<M> result(capacity);
        result.resize(src.size());
        if (!src.empty()) std::copy_n(src.data(),src.size(),result.data());
        return result;
    }
    static int convolution_capacity(const fps<M>& a, const fps<M>& b) {
        if (a.empty() || b.empty() || std::min(a.size(),b.size()) <= 60) return 0;
        if (std::max(a.size(),b.size())/std::min(a.size(),b.size()) >= 16) return 0;
        return upperpow2(a.size()+b.size()-1);
    }
    fps<M> convolution_owned(fps<M> F1, fps<M> F2) const {
        int n = F1.size(), m = F2.size();
        if (n == 0 || m == 0) return fps<M>(0);
        if (std::min(n,m) <= 60) {
            if (n < m) { std::swap(n,m); std::swap(F1,F2); }
            fps<M> ans(n+m-1);
            auto* c = &ans[0];
            const auto* a = &F1[0];
            const auto* b = &F2[0];
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < m; ++j) c[i+j] += a[i]*b[j];
            return ans;
        }
        if (std::max(n,m)/std::min(n,m) >= 16) {
            if (n < m) { std::swap(n,m); std::swap(F1,F2); }
            // Overlap-add: reuse the short operand's transform for each block
            // instead of padding it to the length of the entire long operand.
            const int length = upperpow2(2*m-1), step = length-m+1;
            F2.f.reserve(length); F2.resize(length);
            inplaceDFT(F2);
            fps<M> block(length), ans(n+m-1);
            const mint<M> scale = mint<M>(length).inv();
            for (int start = 0; start < n; start += step) {
                const int count = std::min(step,n-start);
                std::copy_n(F1.data()+start,count,block.data());
                std::fill(block.data()+count,block.data()+length,mint<M>(0));
                inplaceDFT(block);
                pointwise_product(block,block,F2,length,scale);
                inplaceIDFT(block);
                for (int i = 0; i < count+m-1; ++i) ans.data()[start+i] += block.data()[i];
            }
            return ans;
        }
        const int size = upperpow2(n+m-1);
        assert(unsigned(size) <= (1u << log_max_length));
        F1.f.reserve(size); F2.f.reserve(size);
        F1.resize(size); F2.resize(size);
        auto* a = &F1[0];
        const auto* b = &F2[0];
        // Two Montgomery reductions cancel R^2. Include 1/size here, so the
        // inverse transform directly returns normalized, ordinary residues.
        const uint scale = ull(R2)*modpow<M>(size,M-2)%M;
#if defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__))
        if (use_avx2()) {
            if (n <= size/2) transform_avx2<false,false,false,true>(a,size);
            else transform_avx2<false,false,false>(a,size);
            if (m <= size/2) transform_avx2<false,false,false,true>(&F2[0],size);
            else transform_avx2<false,false,false>(&F2[0],size);
            multiply_avx2<false>(a,a,b,size,scale);
            transform_avx2<true,true>(a,size);
        } else
#endif
        {
            inplaceDFT(F1); inplaceDFT(F2);
            for (int i = 0; i < size; ++i)
                a[i].val = montgomery_mul(montgomery_mul(a[i].val,b[i].val),scale);
            inplaceIDFT(F1);
        }
        F1.resize(n+m-1);
        return F1;
    }
public:
    // Copy lvalues directly into the final NTT allocation. Rvalues retain their
    // storage, so std::move still avoids copying when the caller is finished.
    fps<M> convolution(const fps<M>& a, const fps<M>& b) const {
        const int capacity = convolution_capacity(a,b);
        if (capacity == 0) return convolution_owned(a,b);
        return convolution_owned(convolution_copy(a,capacity),convolution_copy(b,capacity));
    }
    fps<M> convolution(fps<M>&& a, fps<M>&& b) const {
        // Moving the same series twice must retain the original square semantics.
        if (&a == &b) return convolution_owned(a,b);
        return convolution_owned(std::move(a),std::move(b));
    }
    fps<M> convolution(fps<M>&& a, const fps<M>& b) const {
        const int capacity = convolution_capacity(a,b);
        auto copy = capacity ? convolution_copy(b,capacity) : b;
        return convolution_owned(std::move(a),std::move(copy));
    }
    fps<M> convolution(const fps<M>& a, fps<M>&& b) const {
        const int capacity = convolution_capacity(a,b);
        auto copy = capacity ? convolution_copy(a,capacity) : a;
        return convolution_owned(std::move(copy),std::move(b));
    }
    fps<M> add(const fps<M>& F1, const fps<M>& F2) const {
        fps<M> ret(max(F1.size(),F2.size()));
        for (uint i = 0; i < F1.size(); ++i) ret[i] = F1[i];
        for (uint i = 0; i < F2.size(); ++i) ret[i] += F2[i];
        return ret;
    }
    fps<M> subtract(const fps<M>& F1, const fps<M>& F2) const {
        fps<M> ret(max(F1.size(),F2.size()));
        for (uint i = 0; i < F1.size(); ++i) ret[i] = F1[i];
        for (uint i = 0; i < F2.size(); ++i) ret[i] -= F2[i];
        return ret;
    }
    /// @brief F*G == 1 mod x^n となるGを求める F[x^0] != 0 が必要 
    fps<M> inv(const fps<M>& F, const int n) const {
        assert(n >= 0);
        if (n == 0) return fps<M>(0);
        assert(!F.empty() && F[0] != 0);
        fps<M> res{F[0].inv()};
        for (int d = 1; d < n; d <<= 1) {
            fps<M> f(2*d), g(2*d);
            for (int j = 0, limit = min<int>(F.size(),2*d); j < limit; ++j) f[j] = F[j];
            for (int j = 0; j < d; ++j) g[j] = res[j];
            inplaceDFT(f); inplaceDFT(g);
            const mint<M> iz = mint<M>(2*d).inv();
            pointwise_product(f,f,g,2*d,iz);
            inplaceIDFT(f);
            for (int j = 0; j < d; ++j) { f[j] = 0; f[j+d] = -f[j+d]; }
            inplaceDFT(f);
            pointwise_product(f,f,g,2*d,iz);
            inplaceIDFT(f);
            for (int j = 0; j < d; ++j) f[j] = res[j];
            res = std::move(f);
        }
        res.resize(n);
        return res;
    }
    /// @brief F÷G の商と余りを求める。{商, 余り}として返される。
    array<fps<M>,2> div(fps<M> F, fps<M>G) const {
        while (!F.f.empty() && F.f.back() == 0){
            F.f.pop();
            F.sz--;
        }
        while (!G.f.empty() && G.f.back() == 0){
            G.f.pop();
            G.sz--;
        }
        int n = F.sz, m = G.sz;

        if (n < m){return {fps<M>{}, F};}

        if (min(m, n-m+1) <= 100){
            fps<M> res(n-m+1);
            mint<M> t = G.f.back().inv();
            for (int i = n-1; i >= m-1; i--){
                auto u = t*F[i];
                res[i-m+1] = u;
                for (int j = 0; j < m; j++){
                    F[i-j] -= G[m-1-j]*u;
                }
            }
            while (!F.f.empty() && F.f.back() == 0){
                F.f.pop();
                F.sz--;
            }
            while (!G.f.empty() && G.f.back() == 0){
                G.f.pop();
                G.sz--;
            }
            return {res, F};
        }

        auto F2 = F;
        auto G2 = G;

        reverse(vall(F.f));
        reverse(vall(G.f));
        while (len(F.f) > n-m+1){
            F.f.pop();
            F.sz--;
        }
        while (len(G.f) > n-m+1){
            G.f.pop();
            G.sz--;
        }
        auto res = convolution(F, inv(G, n-m+1));
        while (len(res.f) > n-m+1){
            res.f.pop();
            res.sz--;
        }
        reverse(vall(res.f));
        auto res2 = subtract(F2, convolution(res, G2));
        while (!res2.f.empty() && res2.f.back() == 0){
            res2.f.pop();
            res2.sz--;
        }
        return {res, res2};
    }
    /// @brief F を G で割った余りを求める。
    fps<M> remainder(fps<M> F, const fps<M>& G) const {
        while (!F.f.empty() && F.f.back() == 0){
            F.f.pop();
            F.sz--;
        }
        assert(!G.f.empty() && G.f.back() != 0);

        int n = F.sz, m = G.sz;
        if (n < m) return F;

        if (min(m, n-m+1) <= 100){
            const bool monic = (G[m-1] == mint<M>(1));
            const mint<M> inverse_leading_coefficient = monic ? mint<M>(1) : G[m-1].inv();
            for (int i = n-1; i >= m-1; i--){
                if (F[i] == 0) continue;
                mint<M> q = monic ? F[i] : inverse_leading_coefficient*F[i];
                for (int j = 0; j < m; j++){
                    F[i-j] -= G[m-1-j]*q;
                }
            }
            F.resize(m-1);
            while (!F.f.empty() && F.f.back() == 0){
                F.f.pop();
                F.sz--;
            }
            return F;
        }

        auto qr = div(std::move(F), G);
        return std::move(qr[1]);
    }
    /// @brief exp(F) mod x^n を求める F[x^0] == 0 が必要
    fps<M> exp(const fps<M>& F, const int n, const mod_table<M>& mtable) const {
        assert(n >= 0);
        if (n == 0) return fps<M>(0);
        assert(F.empty() || F[0] == 0);
        ull sz_f = F.sz;
        fps<M> F0{1},f0{1},g0{1};
        for(int d = 1; d < n; d<<=1){
            auto G0 = g0;
            inplaceDFT(G0);
            const mint<M> iz = mint<M>(d).inv(), iz2 = mint<M>(2*d).inv();
            fps<M> Delta(d);
            pointwise_product(Delta,F0,G0,d,iz);
            inplaceIDFT(Delta);
            Delta[0] -= 1;
            fps<M> delta(2*d);
            for(int j = 0; j < d; j++) delta[d+j] = Delta[j];
            fps<M> epsilon(2*d);

            fps<M> DF0(d-1);
            for (ull j = 0, ulld = d; j < ulld-1; j++) DF0[j] = f0[j+1]*mint<M>(j+1);
            DF0.f.append(0);
            DF0.sz++;
            inplaceDFT(DF0);
            pointwise_product(DF0,DF0,G0,d,iz);
            inplaceIDFT(DF0);
            for (ull j = 0, limit = d-1; j < limit; ++j) {
                const auto df = j+1 < sz_f ? F[j+1]*mint<M>(j+1) : mint<M>(0);
                epsilon[j] = df;
                epsilon[j+d] = DF0[j]-df;
            }
            epsilon[d-1] = DF0[d-1];
            Delta = std::move(delta);
            inplaceDFT(Delta);
            fps<M> DH0(d-1);
            for (ull j = 0, limit = d-1; j < limit; j++) DH0[j] = j+1 < sz_f ? F[j+1]*mint<M>(j+1) : mint<M>(0);
            DH0.resize(2*d);
            inplaceDFT(DH0);
            pointwise_product(Delta,Delta,DH0,2*d,iz2);
            inplaceIDFT(Delta);
            for(int j = 0; j < d; j++){
                epsilon[j+d] -= Delta[j+d];
            }
        
            for (ull i = 2*d-1; i >= 1; i--){
                epsilon[i] = epsilon[i-1]*mtable.inverse_mod(i)-(i < sz_f ? F[i] : mint<M>(0));
            }
            epsilon[0] = 0;
            
            auto Epsilon = std::move(epsilon);
            inplaceDFT(Epsilon);
            for (int j = 0; j < d; j++){DH0[j] = f0[j]; DH0[j+d].val = 0;}
            inplaceDFT(DH0);
            pointwise_product(Epsilon,Epsilon,DH0,2*d,iz2);
            inplaceIDFT(Epsilon);
            f0.resize(2*d);
            for (int j = 0; j < d; j++){
                f0[j+d] -= Epsilon[j+d];
            }
            if(2*d >= n) break;
            
            G0.resize(2*d);
            for (int j = 0; j < d; j++) G0[j] = g0[j];
            inplaceDFT(G0);
            F0 = f0;
            inplaceDFT(F0);
            fps<M> T(2*d);
            pointwise_product(T,F0,G0,2*d,iz2);
            inplaceIDFT(T);
            for (int j = 0; j < d; j++){
                T[j].val = 0;
                T[j+d] = -T[j+d];
            }
            inplaceDFT(T);
            pointwise_product(T,T,G0,2*d,iz2);
            inplaceIDFT(T);
            for (int j = 0; j < d; j++) T[j] = g0[j];
            g0 = std::move(T);
        }
        f0.resize(n);
        return f0;
    }
    /// @brief log(F) mod x^n を求める F[x^0] == 1 が必要 
    fps<M> log(const fps<M>& F, const int n, const mod_table<M>& mtable) const {
        assert(n >= 0);
        if (n == 0) return fps<M>(0);
        assert(!F.empty() && F[0] == 1);
        const int degree = min<int>(F.size(),n)-1;
        fps<M> derivative(degree);
        for (int i = 0; i < degree; ++i) derivative[i] = F[i+1]*mint<M>(i+1);
        auto res = convolution(std::move(derivative),inv(F,n));
        res.resize(n);
        for (int i = n-1; i >= 1; --i) res[i] = res[i-1]*mtable.inverse_mod(i);
        res[0] = 0;
        return res;
    }
    /// @brief F^k mod x^n を求める
    fps<M> pow(const fps<M>& F, const ull k, const int n, const mod_table<M>& mtable) const {
        assert(n >= 0);
        if (n == 0) return fps<M>(0);
        fps<M> ret(n);
        if (k == 0) { ret[0] = 1; return ret; }
        int lowest = 0;
        while (lowest < min<int>(F.size(),n) && F[lowest] == 0) ++lowest;
        if (lowest == min<int>(F.size(),n)) return ret;
        // Check before multiplying, including k = ULLONG_MAX.
        if (lowest > 0 && k > ull((n-1)/lowest)) return ret;
        const int offset = ull(lowest)*k;
        const int precision = n-offset;
        const auto leading = F[lowest], inverse_leading = leading.inv();
        fps<M> g(precision);
        for (int i = 0, limit = min<int>(precision,F.size()-lowest); i < limit; ++i)
            g[i] = F[lowest+i]*inverse_leading;
        auto logarithm = log(g,precision,mtable);
        const mint<M> exponent = k;
        for (auto& x : logarithm.f) x *= exponent;
        auto result = exp(logarithm,precision,mtable);
        const auto scale = leading.pow(k);
        for (int i = 0; i < precision; ++i) ret[offset+i] = result[i]*scale;
        return ret;
    }
    // f(g(x)) mod x^len(f) を求める
    fps<M> composition(fps<M> f, fps<M> g, const mod_table<M>& mtable) const {
        f.resize(max(f.size(),g.size()));
        int N = f.size();
        if (N == 0) return fps<M>(0);

        int n = 1;
        while (n < (int)f.size()) n *= 2;
        f.resize(n), g.resize(n);

        fps<M> W(2*n);
        {
            // bit reverse order
            list<int> btr(2*n);
            int log = 31-__builtin_clz(2*n);
            for (int i = 0; i < 2*n; i++){btr[i] = (btr[i >> 1] >> 1) + ((i & 1) << (log - 1));}
            int t = log_max_length;
            uint r = last_powroot;
            ull dw = modpow<M>(inverse_mod<ll,ll>(r, M), (1<<t)/(4*n));
            uint w = 1;
            for (auto i: btr) { W[i].val = w, w = w*dw%M; }
        }

        auto rec = [&](auto &rec, int n, int k, fps<M> &Q) -> fps<M> {
            if (n == 1) {
                reverse(f.f.begin(), f.f.end());
                fps<M> p(2*k);
                for (int i = 0; i < k; i++) p[2 * i] = f[i];
                return p;
            }
            Q.resize(4*n*k);
            Q[2*n*k].val = 1;
            inplaceDFT(Q);
            fps<M> nxt_Q(2*n*k);
            for (int i = 0; i < 2*n*k; i++) nxt_Q[i] = Q[2*i]*Q[2*i+1];
            
            inplaceIDFT(nxt_Q);
            ull iz = inverse_mod<ll,ll>(nxt_Q.size(), M);
            for (uint i = 0; i < nxt_Q.size(); i++) nxt_Q[i].val = nxt_Q[i].val*iz%M;

            for (int j = 0; j < 2*k; j++) for (int i = n/2; i < n; i++) nxt_Q[n*j+i].val = 0;
            nxt_Q[0].val = 0;
            fps<M> p = rec(rec, n/2, 2*k, nxt_Q);
            for (int j = 0; j < 2*k; j++) for (int i = n/2; i < n; i++) p[n * j + i].val = 0;
            
            iz = mtable.inverse_mod(p.size()).val;
            for (uint i = 0; i < p.size(); i++) p[i].val = p[i].val*iz%M;
            inplaceIDFT_T(p);
            
            p.resize(4*n*k);
            for (int i = 2*n*k-1; i >= 0; i--){
                const auto t = mtable.inverse_mod(2)*W[i]*p[i];
                p[2*i+1] = -t*Q[2*i];
                p[2*i] = t*Q[2*i+1];
            }
            
            inplaceDFT_T(p);
            
            p.resize(2*n*k);
            return p;
        };

        fps<M> Q(2*n);
        for(int i = 0; i < n; i++) Q[i].val = g[i].val == 0 ? 0 : M-g[i].val;
        fps<M> p = rec(rec, n, 1, Q);
        p.resize(n);
        reverse(p.f.begin(), p.f.end());
        p.resize(N);
        return p;
    }
    fps<M> power_projection(fps<M> f, fps<M> wt, int m, const mod_table<M>& mtable) const {
        assert(len(f) == len(wt));
        if (f.empty()) { return fps<M>(m + 1, 0); }
        if (f[0] != 0) {
            auto c = f[0];
            f[0] = 0;
            fps<M> A = power_projection(f, wt, m, mtable);
            for (int p = 0; p < m + 1; p++) A[p] *= mtable.factorialmodinv(p);
            fps<M> B(m + 1);
            mint<M> pow = 1;
            for (int q = 0; q < m + 1; q++) B[q] = pow * mtable.factorialmodinv(q), pow *= c;
            A = convolution(A, B);
            A.resize(m + 1);
            for (int i = 0; i < m + 1; i++) A[i] *= mtable.factorialmod(i);
            return A;
        }

        int n = 1;
        while (n < len(f)) n *= 2;

        for (auto& x: f.f) x = -x;
        f.resize(n), wt.resize(n);
        reverse(vall(wt.f));
        fps<M>&P = wt, &Q = f;
        P.resize(4 * n), Q.resize(4 * n);

        fps<M> W(n);
        {
            // bit reverse order
            list<int> btr(n);
            int log = 31-__builtin_clz(n);
            for (int i = 1; i < n; i++) { btr[i] = (btr[i >> 1] >> 1) + ((i & 1) << (log - 1)); }
            int t = log_max_length;
            mint<M> r = last_powroot;
            mint<M> dw = r.inv().pow((1 << t) / (2 * n));
            mint<M> w = 1;
            for (auto i: btr) { W[i] = w, w *= dw; }
        }

        int k = 1;
        while (n > 1) {
            auto normalized_IDFT = [&](fps<M>& A) {
                inplaceIDFT(A);
                mint<M> iz = mint<M>(A.size()).inv();
                for (auto& x: A.f) x *= iz;
            };
            auto doubling_y = [&](fps<M>& A, int l, int r) {
                mint<M> z = W[k / 2].inv();
                fps<M> t(k);
                for (int i = l; i < r; i++) {
                    for (int j = 0; j < k; j++) t[j] = A[2 * n * j + i];
                    normalized_IDFT(t);
                    mint<M> x = 1;
                    for (int j = 0; j < k; j++) t[j] *= x, x *= z;
                    inplaceDFT(t);
                    for (int j = 0; j < k; j++) A[2 * n * (k + j) + i] = t[j];
                }
            };
            auto FFT_x = [&](fps<M>& A, int l, int r) {
                fps<M> t(2 * n);
                for (int j = l; j < r; j++) {
                    copy(A.f.begin() + 2 * n * j, A.f.begin() + 2 * n * (j + 1), t.f.begin());
                    inplaceDFT(t);
                    copy(vall(t.f), A.f.begin() + 2 * n * j);
                }
            };

            if (n <= k) {
                doubling_y(P, 0, n), doubling_y(Q, 1, n);
                FFT_x(P, 0, 2 * k), FFT_x(Q, 0, 2 * k);
            } else {
                FFT_x(P, 0, k), FFT_x(Q, 0, k);
                doubling_y(P, 0, 2 * n), doubling_y(Q, 0, 2 * n);
            }
            for (int i = 0; i < 2 * n * k; i++) Q[i] += 1;
            for (int i = 2 * n * k; i < 4 * n * k; i++) Q[i] -= 1;

            fps<M> F(2 * n), G(2 * n), p(n), q(n);
            for (int j = 0; j < 2 * k; j++) {
                copy(P.f.begin() + 2 * n * j, P.f.begin() + 2 * n * (j + 1), F.f.begin());
                copy(Q.f.begin() + 2 * n * j, Q.f.begin() + 2 * n * (j + 1), G.f.begin());
                for (int i = 0; i < n; i++) {
                    p[i] = W[i] * (F[2 * i] * G[2 * i + 1] - F[2 * i + 1] * G[2 * i]);
                    q[i] = G[2 * i] * G[2 * i + 1];
                }
                normalized_IDFT(p), normalized_IDFT(q);
                fill(p.f.begin() + n / 2, p.f.end(), 0);
                fill(q.f.begin() + n / 2, q.f.end(), 0);
                copy(vall(p.f), P.f.begin() + n * j);
                copy(vall(q.f), Q.f.begin() + n * j);
            }
            fill(P.f.begin() + 2 * n * k, P.f.end(), 0);
            fill(Q.f.begin() + 2 * n * k, Q.f.end(), 0);
            for (int j = 0; j < 4 * k; j++) Q[n * j] = 0;
            n /= 2, k *= 2;
        }
        for (int i = 0; i < k; i++) P[i] = P[2 * i];
        P.resize(k);
        mint<M> ik = mint<M>(k).inv();
        for (auto& x: P.f) x *= ik;
        inplaceIDFT(P);
        for (auto& x: P.f) x *= ik;
        reverse(vall(P.f));
        P.resize(m + 1);
        return P;
    }
    fps<M> compositional_inverse(fps<M> f, mod_table<M>& mtable) const {
        /*
        仮定：[x^1]f = c は 0 ではない.
        f を改めて cf と書く. f の逆関数を g とすれば
        g(cf(x)/c)=x なので g(x/c) が求めるもの．よって [x^1]f=1 に帰着.

        f, g は n 次とする. Lagrange inversion formula
        n[x^n]f(x)^i = [x^{n-i}] i(g(x)/x)^{-n}
        左辺を i=1,...,n で求めると (g(x)/x)^{-n} が n-1 次まで求まる.
        [x^1]g=1 に帰着していたので -1/n 乗根がとれて, g(x)/x が n-1 次まで求まる.
        */
        const int n = len(f) - 1;
        if (n == -1) return {};
        assert(f[0] == 0);
        if (n == 0) return f;
        assert(f[1] != 0);
        mint<M> c = f[1];
        mint<M> ic = c.inv();
        for (auto& x: f.f) x *= ic;
        fps<M> wt(n + 1);
        wt[n] = 1;

        fps<M> A = power_projection(f, wt, n, mtable);
        fps<M> g(n);
        for (int i = 1; i < n + 1; i++) g[n - i] = mint<M>(n) * A[i] * mint<M>(i).inv();
        g = pow(g, M - mint<M>(n).inv().val, len(g), mtable);
        g.f.insert(g.f.begin(), 0);
        g.sz++;

        mint<M> pow = 1;
        for (int i = 0; i < len(g); i++) g[i] *= pow, pow *= ic;
        return g;
    }
    // f(x+a)を求める
    fps<M> polynominal_taylor_shift(const fps<M>& F, ll a, const mod_table<M>& mtable){
        const int n = F.size();
        const mint<M> shift = a;
        if (n == 0 || shift == 0) return F;
        fps<M> g(n), h(n);
        mint<M> power = 1;
        for (int i = 0; i < n; ++i) {
            g[n-1-i] = F[i]*mtable.factorialmod(i);
            h[i] = power*mtable.factorialmodinv(i);
            power *= shift;
        }
        g = convolution(std::move(g),std::move(h));
        fps<M> result(n);
        for (int i = 0; i < n; ++i) result[i] = g[n-1-i]*mtable.factorialmodinv(i);
        return result;
    }
    // fの微分を求める。
    fps<M> differential(const fps<M>& F) const {
        fps<M> G = F;
        G.differential();
        return G;
    }
    // fの積分を求める。積分定数は0
    fps<M> integral(const fps<M>& F, const mod_table<M>& mtable) const {
        fps<M> G = F;
        G.integral(mtable);
        return G;
    }
};



/// @brief `[x^N].val(P(x)/Q(x))` をmod Mで求める
template<uint M>
constant_modint<M> Bostan_Mori(const ll N, FormalPowerSeries<M> P, FormalPowerSeries<M> Q, const fps_operator<M>& op){
    assert(N >= 0);
    if (N == 0){
        return P[0].val*inverse_mod<ll,ll>(Q[0].val, M)%M;
    }
    FormalPowerSeries<M> Q_minus;
    const int maxloop = (N == 1 ? 1 : 65-__builtin_clzll(N-1));
    for (int _i_ = 0; _i_ < maxloop; _i_++){
        Q_minus.resize(Q.size());
        for (size_t i = 0; i < Q_minus.size(); i += 2){
            Q_minus[i] = Q[i];
        }
        for (size_t i = 1; i < Q_minus.size(); i += 2){
            Q_minus[i].val = (M-1)*Q[i].val%M;
        }
        auto A = op.convolution(P,Q_minus);
        auto B = op.convolution(Q,Q_minus);
        Q.resize((B.size()+1)/2);
        for (size_t i = 0; i < B.size(); i += 2){
            Q[i/2] = B[i];
        }

        P.resize((A.size()+!btest_for_fps(N,_i_))/2);
        for (size_t i = btest_for_fps(N,_i_); i < A.size(); i += 2){
            P[i/2] = A[i];
        }
    }
    return P[0].val*inverse_mod<ll,ll>(Q[0].val, M)%M;
}

template<uint M>
list<FormalPowerSeries<M>> multipoint_product_tree(const list<constant_modint<M>>& samples, const fps_operator<M>& op){
    const int n = len(samples);
    if (n == 0) return {};

    list<FormalPowerSeries<M>> tree(2*n-1);
    int next_id = 0;
    auto build = [&](auto&& self, int l, int r) -> int {
        const int node = next_id++;
        if (r-l == 1){
            tree[node] = FormalPowerSeries<M>{-samples[l], constant_modint<M>(1)};
            return node;
        }

        const int mid = (l+r)/2;
        const int left = self(self, l, mid);
        const int right = self(self, mid, r);
        if (r-l == 2){
            const auto x0 = samples[l], x1 = samples[l+1];
            tree[node] = FormalPowerSeries<M>{x0*x1, -(x0+x1), constant_modint<M>(1)};
        }
        else if (r-l == 3){
            const auto x0 = samples[l], x1 = samples[l+1], x2 = samples[l+2];
            tree[node] = FormalPowerSeries<M>{
                -x0*x1*x2,
                x0*x1+x0*x2+x1*x2,
                -(x0+x1+x2),
                constant_modint<M>(1)
            };
        }
        else{
            tree[node] = op.convolution(tree[left], tree[right]);
        }
        return node;
    };
    build(build, 0, n);
    return tree;
}

template<uint M>
list<constant_modint<M>> multipoint_evaluation_with_product_tree(const FormalPowerSeries<M>& F, const list<constant_modint<M>>& samples, const list<FormalPowerSeries<M>>& tree, const fps_operator<M>& op){
    const int n = len(samples);
    if (n == 0) return {};

    list<constant_modint<M>> ans(n);
    list<FormalPowerSeries<M>> remainders;
    remainders.reserve(n);
    list<array<int, 4>> stack;
    stack.reserve(64);

    remainders.append(F);
    stack.append({0,n,0,0});

    while (!stack.empty()){
        auto [l,r,node,remainder_index] = stack.pop();


        if (l >= r){continue;}
        if (r-l == 1){
            constant_modint<M> value = 0;
            for (int i = int(remainders[remainder_index].sz)-1; i >= 0; i--){
                value = value*samples[l]+remainders[remainder_index][i];
            }
            ans[l] = value;
            continue;
        }

        if (remainders[remainder_index].size() >= tree[node].size()){
            remainders.append(op.remainder(remainders[remainder_index], tree[node]));
            remainder_index = int(remainders.size())-1;
        }
        const int mid = (l+r)/2;
        const int left = node+1;
        const int right = node+2*(mid-l);
        stack.append({mid,r,right,remainder_index});
        stack.append({l,mid,left,remainder_index});
    }

    return ans;
}

template <uint M>
list<constant_modint<M>> multipoint_evaluation(const FormalPowerSeries<M>& F, const list<constant_modint<M>>& samples, const fps_operator<M>& op){
    auto tree = multipoint_product_tree(samples, op);
    return multipoint_evaluation_with_product_tree(F, samples, tree, op);
}

template<uint M>
FormalPowerSeries<M> product_of_polynomial_sequence(const list<FormalPowerSeries<M>>& A, const fps_operator<M>& op){
    delist<FormalPowerSeries<M>> dq;
    dq.append(FormalPowerSeries<998244353>{1});
    for (auto& f : A){
        dq.append(std::move(f));
    }
    while (dq.size() > 1){
        dq.append(op.convolution(std::move(*dq.begin()), std::move(*++dq.begin())));
        dq.popleft();
        dq.popleft();
    }
    return dq.front();
}

template<uint M>
FormalPowerSeries<M> polynomial_interpolation(const list<constant_modint<M>>& X, const list<constant_modint<M>>& Y, const fps_operator<M>& op){
    assert(len(X) == len(Y));
    const int N = len(X);
    if (N == 0) return {};

    const auto tree = multipoint_product_tree(X, op);
    auto product_derivative = op.differential(tree[0]);
    list<constant_modint<M>> Z =
        multipoint_evaluation_with_product_tree(product_derivative, X, tree, op);

    list<constant_modint<M>> prefix_product(N+1);
    prefix_product[0] = 1;
    for (int i = 0; i < N; i++){
        prefix_product[i+1] = prefix_product[i]*Z[i];
    }
    constant_modint<M> inverse_suffix_product = prefix_product[N].inv();
    for (int i = N-1; i >= 0; i--){
        const auto denominator = Z[i];
        Z[i] = Y[i]*prefix_product[i]*inverse_suffix_product;
        inverse_suffix_product *= denominator;
    }

    auto interpolate = [&](auto&& self, int l, int r, int node) -> FormalPowerSeries<M> {
        if (r-l == 1){
            return FormalPowerSeries<M>{Z[l]};
        }
        if (r-l == 2){
            const auto x0 = X[l], x1 = X[l+1];
            const auto z0 = Z[l], z1 = Z[l+1];
            return FormalPowerSeries<M>{-(z0*x1+z1*x0), z0+z1};
        }
        if (r-l == 3){
            const auto x0 = X[l], x1 = X[l+1], x2 = X[l+2];
            const auto z0 = Z[l], z1 = Z[l+1], z2 = Z[l+2];
            return FormalPowerSeries<M>{
                z0*x1*x2+z1*x0*x2+z2*x0*x1,
                -(z0*(x1+x2)+z1*(x0+x2)+z2*(x0+x1)),
                z0+z1+z2
            };
        }

        const int mid = (l+r)/2;
        const int left = node+1;
        const int right = node+2*(mid-l);
        auto left_interpolation = self(self, l, mid, left);
        auto right_interpolation = self(self, mid, r, right);
        return op.add(
            op.convolution(std::move(left_interpolation), tree[right]),
            op.convolution(tree[left], std::move(right_interpolation))
        );
    };
    return interpolate(interpolate, 0, N, 0);
}

#undef vall
#endif /* FPS_OPERATION_HPP_ */


using mint = constant_modint<998244353>;
fps_operator<998244353> op;
mod_table<998244353> mtable(1000000);
fastio::FastInput in;
fastio::FastOutput out;


void solve(){

    ll N,M,K;
    in.read(N,M,K);

    mint nCr = 1;

    FormalPowerSeries<998244353> F(K+1);

    for (ll i = 0; i <= min(M, K); i++){
        F[i] = nCr*mint(M-i).pow(N);
        nCr *= mint(M-i)*mtable.inverse_mod(i+1);
    }

    // out.print(F);

    FormalPowerSeries<998244353> G(K+1);
    G[0] = 0;
    for (ll i = 1; i <= K; i++){
        G[i] = mtable.factorialmodinv(i);
        if (i%2){G[i] *= -1;}
    }

    // out.print(G);


    F = op.composition(F, G, mtable);

    for (ll i = 0; i <= K; i++){
        G[i] = mint(M).pow(i)*mtable.factorialmodinv(i);
    }

    F = op.convolution(F,G);

    for (ll i = 1; i <= K; i++){
        out.print(F[i]*mtable.factorialmod(i));
    }

    
}





int main(){

    ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    ll T = 1;
    // cin >> T;
    while (T--){
        solve();
    }
}