#include #include #pragma GCC target("avx2") #pragma GCC optimize("O3") #pragma GCC optimize("unroll-loops") #define FOR(i,n) for(int i = 0; i < (n); i++) #define sz(c) ((int)(c).size()) #define ten(x) ((int)1e##x) #define all(v) (v).begin(), (v).end() using namespace std; using ll=long long; using P = pair; const long double PI=acos(-1); const ll INF=1e18; const int inf=1e9; struct Edge { ll to; ll cost; }; using Graph=vector>; template bool chmax(T &a,const T& b){ if (a bool chmin(T &a,const T& b){ if (a>b){ a=b; return true; } return false; } template struct Fp{ ll val; constexpr Fp(long long v = 0) noexcept : val(v % MOD) { if (val < 0) val += MOD; } static constexpr int getmod() { return MOD; } constexpr Fp operator - () const noexcept { return val ? MOD - val : 0; } constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; } constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; } constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; } constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; } constexpr Fp& operator += (const Fp& r) noexcept { val += r.val; if (val >= MOD) val -= MOD; return *this; } constexpr Fp& operator -= (const Fp& r) noexcept { val -= r.val; if (val < 0) val += MOD; return *this; } constexpr Fp& operator *= (const Fp& r) noexcept { val = val * r.val % MOD; return *this; } constexpr Fp& operator /= (const Fp& r) noexcept { ll a = r.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); } val = val * u % MOD; if (val < 0) val += MOD; return *this; } constexpr bool operator == (const Fp& r) const noexcept { return this->val == r.val; } constexpr bool operator != (const Fp& r) const noexcept { return this->val != r.val; } constexpr bool operator < (const Fp& r) const noexcept { return this->val < r.val; } friend constexpr istream& operator >> (istream& is, Fp& x) noexcept { is >> x.val; x.val %= MOD; if (x.val < 0) x.val += MOD; return is; } friend constexpr ostream& operator << (ostream& os, const Fp& x) noexcept { return os << x.val; } friend constexpr Fp modpow(const Fp& a, long long n) noexcept { Fp res=1,r=a; while(n){ if(n&1) res*=r; r*=r; n>>=1; } return res; } friend constexpr Fp modinv(const Fp& r) noexcept { long long a = r.val, b = MOD, u = 1, v = 0; while (b) { long long t = a / b; a -= t * b, swap(a, b); u -= t * v, swap(u, v); } return Fp(u); } ll get(){ return val; } explicit operator bool()const{ return val; } }; template< uint32_t mod, bool fast = false > struct MontgomeryModInt { using mint = MontgomeryModInt; using i32 = int32_t; using i64 = int64_t; using u32 = uint32_t; using u64 = uint64_t; static constexpr u32 get_r() { u32 ret = mod; for(i32 i = 0; i < 4; i++) ret *= 2 - mod * ret; return ret; } static constexpr u32 r = get_r(); static constexpr u32 n2 = -u64(mod) % mod; static_assert(r * mod == 1, "invalid, r * mod != 1"); static_assert(mod < (1 << 30), "invalid, mod >= 2 ^ 30"); static_assert((mod & 1) == 1, "invalid, mod % 2 == 0"); u32 a; MontgomeryModInt() : a{} {} MontgomeryModInt(const i64 &x) : a(reduce(u64(fast ? x : (x % mod + mod)) * n2)) {} static constexpr u32 reduce(const u64 &b) { return u32(b >> 32) + mod - u32((u64(u32(b) * r) * mod) >> 32); } constexpr mint& operator+=(const mint &p) { if(i32(a += p.a - 2 * mod) < 0) a += 2 * mod; return *this; } constexpr mint& operator-=(const mint &p) { if(i32(a -= p.a) < 0) a += 2 * mod; return *this; } constexpr mint& operator*=(const mint &p) { a = reduce(u64(a) * p.a); return *this; } constexpr mint& operator/=(const mint &p) { *this *= modinv(p); return *this; } constexpr mint operator-() const { return mint() - *this; } constexpr mint operator+(const mint &p) const { return mint(*this) += p; } constexpr mint operator-(const mint &p) const { return mint(*this) -= p; } constexpr mint operator*(const mint &p) const { return mint(*this) *= p; } constexpr mint operator/(const mint &p) const { return mint(*this) /= p; } constexpr bool operator==(const mint &p) const { return (a >= mod ? a - mod : a) == (p.a >= mod ? p.a - mod : p.a); } constexpr bool operator!=(const mint &p) const { return (a >= mod ? a - mod : a) != (p.a >= mod ? p.a - mod : p.a); } u32 get() const { u32 ret = reduce(a); return ret >= mod ? ret - mod : ret; } friend constexpr MontgomeryModInt modpow(const MontgomeryModInt &x,u64 n) noexcept { MontgomeryModInt ret(1), mul(x); while(n > 0) { if(n & 1) ret *= mul; mul *= mul; n >>= 1; } return ret; } friend constexpr MontgomeryModInt modinv(const MontgomeryModInt &r) noexcept { u64 a = r.get(), b = mod, u = 1, v = 0; while (b) { long long t = a / b; a -= t * b, swap(a, b); u -= t * v, swap(u, v); } return MontgomeryModInt(u); } friend ostream &operator<<(ostream &os, const mint &p) { return os << p.get(); } friend istream &operator>>(istream &is, mint &a) { i64 t; is >> t; a = mint(t); return is; } static constexpr u32 getmod() { return mod; } }; template struct SegmentTree{ int n; vector dat; SegmentTree(int N){ n=1; while(n>=1; dat[k]=op(dat[k*2],dat[k*2+1]); } } void apply(int k,T x){ k+=n; dat[k]=op(dat[k],x); while(k){ k>>=1; dat[k]=op(dat[k*2],dat[k*2+1]); } } void set(int k,T x){ k+=n; dat[k]=x; while(k){ k>>=1; dat[k]=op(dat[k*2],dat[k*2+1]); } } T query(int l,int r){ T prodl=e(),prodr=e(); l+=n; r+=n; while(l>=1; r>>=1; } return op(prodl,prodr); } }; struct FenwickTree{ int n; vector dat; FenwickTree(int N){ n=1; while(n struct LazySegTree{ private: int _n,size=1,idx=0; vectorseq; vectorlazy; void update(int k){seq[k]=op(seq[2*k],seq[2*k+1]);} void all_apply(int k,F f){ seq[k]=mapping(f,seq[k]); if(k(n,e())){} LazySegTree(const vector&v):_n(int(v.size())){ while(size<_n)size<<=1,idx++; seq=vector(2*size,e()); lazy=vector(2*size,id()); for(int i=0;i<_n;i++)seq[size+i]=v[i]; for(int i=size-1;i>=1;i--)update(i); } void set(int p,S x){ p+=size; for(int i=idx;i>=1;i--)eval(p>>i); seq[p]=x; for(int i=1;i<=idx;i++)update(p>>i); } void add(int p,ll x){ p+=size; for(int i=idx;i>=1;i--)eval(p>>i); seq[p].value+=x; for(int i=1;i<=idx;i++)update(p>>i); } S operator[](int p){ p+=size; for(int i=idx;i>=1;i--)eval(p>>i); return seq[p]; } S query(int l,int r){ if(l==r)return e(); S sml=e(),smr=e(); l+=size,r+=size; for(int i=idx;i>=1;i--){ if(((l>>i)<>i); if(((r>>i)<>i); } while(l>=1,r>>=1; } return op(sml,smr); } S all_query()const{return seq[1];} void apply(int p,F f){ p+=size; for(int i=idx;i>=1;i--)eval(p>>i); seq[p]=mapping(f,seq[p]); for(int i=1;i<=idx;i++)update(p>>i); } void apply(int l,int r,F f){ if(l==r)return ; l+=size; r+=size; for(int i=idx;i>=1;i--){ if(((l>>i)<>i); if(((r>>i)<>i); } int l2=l,r2=r; while(l>=1; r>>=1; } l=l2,r=r2; for(int i=1;i<=idx;i++){ if(((l>>i)<>i); if(((r>>i)<>i); } } }; ll mod(ll a,ll MOD){ if(a<0) a+=MOD; return a%MOD; } ll modpow(ll a,ll n,ll mod){ ll res=1; a%=mod; while (n>0){ if (n & 1) res*=a; a *= a; a%=mod; n >>= 1; res%=mod; } return res; } vector

prime_factorize(ll N) { vector

res; for (ll a = 2; a * a <= N; ++a) { if (N % a != 0) continue; ll ex = 0; while(N % a == 0){ ++ex; N /= a; } res.push_back({a, ex}); } if (N != 1) res.push_back({N, 1}); return res; } ll modinv(ll a, ll mod) { ll b = mod, u = 1, v = 0; while (b) { ll t = a/b; a -= t * b, swap(a, b); u -= t * v, swap(u, v); } u %= mod; if (u < 0) u += mod; return u; } ll extGcd(ll a, ll b, ll &p, ll &q) { if (b == 0) { p = 1; q = 0; return a; } ll d = extGcd(b, a%b, q, p); q -= a/b * p; return d; } P ChineseRem(const vector &b, const vector &m) { ll r = 0, M = 1; for (int i = 0; i < (int)b.size(); ++i) { ll p, q; ll d = extGcd(M, m[i], p, q); if ((b[i] - r) % d != 0) return make_pair(0, -1); ll tmp = (b[i] - r) / d * p % (m[i]/d); r += M * tmp; M *= m[i]/d; } return make_pair(mod(r, M), M); } namespace NTT { using i64 = int64_t; __attribute__((target("sse4.2"))) inline __m128i my128_mullo_epu32( const __m128i &a, const __m128i &b) { return _mm_mullo_epi32(a, b); } __attribute__((target("sse4.2"))) inline __m128i my128_mulhi_epu32( const __m128i &a, const __m128i &b) { __m128i a13 = _mm_shuffle_epi32(a, 0xF5); __m128i b13 = _mm_shuffle_epi32(b, 0xF5); __m128i prod02 = _mm_mul_epu32(a, b); __m128i prod13 = _mm_mul_epu32(a13, b13); __m128i prod = _mm_unpackhi_epi64(_mm_unpacklo_epi32(prod02, prod13), _mm_unpackhi_epi32(prod02, prod13)); return prod; } __attribute__((target("sse4.2"))) inline __m128i montgomery_mul_128( const __m128i &a, const __m128i &b, const __m128i &r, const __m128i &m1) { return _mm_sub_epi32( _mm_add_epi32(my128_mulhi_epu32(a, b), m1), my128_mulhi_epu32(my128_mullo_epu32(my128_mullo_epu32(a, b), r), m1)); } __attribute__((target("sse4.2"))) inline __m128i montgomery_add_128( const __m128i &a, const __m128i &b, const __m128i &m2, const __m128i &m0) { __m128i ret = _mm_sub_epi32(_mm_add_epi32(a, b), m2); return _mm_add_epi32(_mm_and_si128(_mm_cmpgt_epi32(m0, ret), m2), ret); } __attribute__((target("sse4.2"))) inline __m128i montgomery_sub_128( const __m128i &a, const __m128i &b, const __m128i &m2, const __m128i &m0) { __m128i ret = _mm_sub_epi32(a, b); return _mm_add_epi32(_mm_and_si128(_mm_cmpgt_epi32(m0, ret), m2), ret); } __attribute__((target("avx2"))) inline __m256i my256_mullo_epu32( const __m256i &a, const __m256i &b) { return _mm256_mullo_epi32(a, b); } __attribute__((target("avx2"))) inline __m256i my256_mulhi_epu32( const __m256i &a, const __m256i &b) { __m256i a13 = _mm256_shuffle_epi32(a, 0xF5); __m256i b13 = _mm256_shuffle_epi32(b, 0xF5); __m256i prod02 = _mm256_mul_epu32(a, b); __m256i prod13 = _mm256_mul_epu32(a13, b13); __m256i prod = _mm256_unpackhi_epi64(_mm256_unpacklo_epi32(prod02, prod13), _mm256_unpackhi_epi32(prod02, prod13)); return prod; } __attribute__((target("avx2"))) inline __m256i montgomery_mul_256( const __m256i &a, const __m256i &b, const __m256i &r, const __m256i &m1) { return _mm256_sub_epi32( _mm256_add_epi32(my256_mulhi_epu32(a, b), m1), my256_mulhi_epu32(my256_mullo_epu32(my256_mullo_epu32(a, b), r), m1)); } __attribute__((target("avx2"))) inline __m256i montgomery_add_256( const __m256i &a, const __m256i &b, const __m256i &m2, const __m256i &m0) { __m256i ret = _mm256_sub_epi32(_mm256_add_epi32(a, b), m2); return _mm256_add_epi32(_mm256_and_si256(_mm256_cmpgt_epi32(m0, ret), m2), ret); } __attribute__((target("avx2"))) inline __m256i montgomery_sub_256( const __m256i &a, const __m256i &b, const __m256i &m2, const __m256i &m0) { __m256i ret = _mm256_sub_epi32(a, b); return _mm256_add_epi32(_mm256_and_si256(_mm256_cmpgt_epi32(m0, ret), m2), ret); } int calc_primitive_root(int mod) { if (mod == 2) return 1; if (mod == 167772161) return 3; if (mod == 469762049) return 3; if (mod == 754974721) return 11; if (mod == 998244353) return 3; int divs[20] = {}; divs[0] = 2; int cnt = 1; long long x = (mod - 1) / 2; while (x % 2 == 0) x /= 2; for (long long i = 3; i * i <= x; i += 2) { if (x % i == 0) { divs[cnt++] = i; while (x % i == 0) x /= i; } } if (x > 1) divs[cnt++] = x; for (int g = 2;; g++) { bool ok = true; for (int i = 0; i < cnt; i++) { if (modpow(g, (mod - 1) / divs[i], mod) == 1) { ok = false; break; } } if (ok) return g; } } int get_fft_size(int N, int M) { int size_a = 1, size_b = 1; while (size_a < N) size_a <<= 1; while (size_b < M) size_b <<= 1; return max(size_a, size_b) << 1; } constexpr int bsf_constexpr(unsigned int n) { int x = 0; while (!(n & (1 << x))) x++; return x; } int bsf(unsigned int n) { #ifdef _MSC_VER unsigned long index; _BitScanForward(&index, n); return index; #else return __builtin_ctz(n); #endif } template struct fft_info{ static constexpr int rank2 = bsf_constexpr(mint::getmod() - 1); std::array root; // root[i]^(2^i) == 1 std::array iroot; // root[i] * iroot[i] == 1 std::array rate2; std::array irate2; std::array rate3; std::array irate3; int g; fft_info(){ int MOD=mint::getmod(); g=calc_primitive_root(MOD); root[rank2] = modpow(mint(g),(MOD - 1) >> rank2); iroot[rank2] = modinv(root[rank2]); for (int i = rank2 - 1; i >= 0; i--) { root[i] = root[i + 1] * root[i + 1]; iroot[i] = iroot[i + 1] * iroot[i + 1]; } { mint prod = 1, iprod = 1; for (int i = 0; i <= rank2 - 2; i++) { rate2[i] = root[i + 2] * prod; irate2[i] = iroot[i + 2] * iprod; prod *= iroot[i + 2]; iprod *= root[i + 2]; } } { mint prod = 1, iprod = 1; for (int i = 0; i <= rank2 - 3; i++) { rate3[i] = root[i + 3] * prod; irate3[i] = iroot[i + 3] * iprod; prod *= iroot[i + 3]; iprod *= root[i + 3]; } } } }; int ceil_pow2(int n) { int x = 0; while ((1U << x) < (unsigned int)(n)) x++; return x; } // number-theoretic transform template void trans(std::vector& a) { int n = int(a.size()); int h = ceil_pow2(n); int MOD=a[0].getmod(); static const fft_info info; int len = 0; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed while (len < h) { if (h - len == 1) { int p = 1 << (h - len - 1); mint rot = 1; for (int s = 0; s < (1 << len); s++) { int offset = s << (h - len); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p] * rot; a[i + offset] = l + r; a[i + offset + p] = l - r; } if (s + 1 != (1 << len)) rot *= info.rate2[bsf(~(unsigned int)(s))]; } len++; } else { // 4-base int p = 1 << (h - len - 2); mint rot = 1, imag = info.root[2]; for (int s = 0; s < (1 << len); s++) { mint rot2 = rot * rot; mint rot3 = rot2 * rot; int offset = s << (h - len); for (int i = 0; i < p; i++) { auto mod2 = 1ULL * MOD * MOD; auto a0 = 1ULL * a[i + offset].get(); auto a1 = 1ULL * a[i + offset + p].get() * rot.get(); auto a2 = 1ULL * a[i + offset + 2 * p].get() * rot2.get(); auto a3 = 1ULL * a[i + offset + 3 * p].get() * rot3.get(); auto a1na3imag = 1ULL * mint(a1 + mod2 - a3).get() * imag.get(); auto na2 = mod2 - a2; a[i + offset] = a0 + a2 + a1 + a3; a[i + offset + 1 * p] = a0 + a2 + (2 * mod2 - (a1 + a3)); a[i + offset + 2 * p] = a0 + na2 + a1na3imag; a[i + offset + 3 * p] = a0 + na2 + (mod2 - a1na3imag); } if (s + 1 != (1 << len)) rot *= info.rate3[bsf(~(unsigned int)(s))]; } len += 2; } } } template void trans_inv(std::vector& a) { int n = int(a.size()); int h = ceil_pow2(n); static const fft_info info; int MOD=a[0].getmod(); int len = h; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed while (len) { if (len == 1) { int p = 1 << (h - len); mint irot = 1; for (int s = 0; s < (1 << (len - 1)); s++) { int offset = s << (h - len + 1); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p]; a[i + offset] = l + r; a[i + offset + p] = (unsigned long long)(MOD + l.get() - r.get()) * irot.get(); ; } if (s + 1 != (1 << (len - 1))) irot *= info.irate2[bsf(~(unsigned int)(s))]; } len--; } else { // 4-base int p = 1 << (h - len); mint irot = 1, iimag = info.iroot[2]; for (int s = 0; s < (1 << (len - 2)); s++) { mint irot2 = irot * irot; mint irot3 = irot2 * irot; int offset = s << (h - len + 2); for (int i = 0; i < p; i++) { auto a0 = 1ULL * a[i + offset + 0 * p].get(); auto a1 = 1ULL * a[i + offset + 1 * p].get(); auto a2 = 1ULL * a[i + offset + 2 * p].get(); auto a3 = 1ULL * a[i + offset + 3 * p].get(); auto a2na3iimag = 1ULL * mint((MOD + a2 - a3) * iimag.get()).get(); a[i + offset] = a0 + a1 + a2 + a3; a[i + offset + 1 * p] = (a0 + (MOD - a1) + a2na3iimag) * irot.get(); a[i + offset + 2 * p] = (a0 + a1 + (MOD - a2) + (MOD - a3)) * irot2.get(); a[i + offset + 3 * p] = (a0 + (MOD - a1) + (MOD - a2na3iimag)) * irot3.get(); } if (s + 1 != (1 << (len - 2))) irot *= info.irate3[bsf(~(unsigned int)(s))]; } len -= 2; } } } namespace ntt_inner { using u64 = uint64_t; constexpr uint32_t get_pr(uint32_t mod) { if (mod == 2) return 1; u64 ds[32] = {}; int idx = 0; u64 m = mod - 1; for (u64 i = 2; i * i <= m; ++i) { if (m % i == 0) { ds[idx++] = i; while (m % i == 0) m /= i; } } if (m != 1) ds[idx++] = m; uint32_t pr = 2; while (1) { int flg = 1; for (int i = 0; i < idx; ++i) { u64 a = pr, b = (mod - 1) / ds[i], r = 1; while (b) { if (b & 1) r = r * a % mod; a = a * a % mod; b >>= 1; } if (r == 1) { flg = 0; break; } } if (flg == 1) break; ++pr; } return pr; } constexpr int SZ_FFT_BUF = 1 << 23; uint32_t _buf1[SZ_FFT_BUF] __attribute__((aligned(64))); uint32_t _buf2[SZ_FFT_BUF] __attribute__((aligned(64))); } // namespace ntt_inner template struct NumberTheoreticTransform { static constexpr uint32_t mod = mint::getmod(); static constexpr uint32_t pr = ntt_inner::get_pr(mint::getmod()); static constexpr int level = __builtin_ctzll(mod - 1); mint dw[level], dy[level]; mint *buf1, *buf2; constexpr NumberTheoreticTransform() { setwy(level); union raw_cast { mint dat; uint32_t _; }; buf1 = &(((raw_cast *)(ntt_inner::_buf1))->dat); buf2 = &(((raw_cast *)(ntt_inner::_buf2))->dat); } constexpr void setwy(int k) { mint w[level], y[level]; w[k - 1] = modpow(mint(pr),(mod - 1) / (1 << k)); y[k - 1] = modinv(w[k - 1]); for (int i = k - 2; i > 0; --i) w[i] = w[i + 1] * w[i + 1], y[i] = y[i + 1] * y[i + 1]; dw[0] = dy[0] = w[1] * w[1]; dw[1] = w[1], dy[1] = y[1], dw[2] = w[2], dy[2] = y[2]; for (int i = 3; i < k; ++i) { dw[i] = dw[i - 1] * y[i - 2] * w[i]; dy[i] = dy[i - 1] * w[i - 2] * y[i]; } } __attribute__((target("avx2"))) void ntt(mint *a, int n) { int k = n ? __builtin_ctz(n) : 0; if (k == 0) return; if (k == 1) { mint a1 = a[1]; a[1] = a[0] - a[1]; a[0] = a[0] + a1; return; } if (k & 1) { int v = 1 << (k - 1); if (v < 8) { for (int j = 0; j < v; ++j) { mint ajv = a[j + v]; a[j + v] = a[j] - ajv; a[j] += ajv; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m2 = _mm256_set1_epi32(mod + mod); int j0 = 0; int j1 = v; for (; j0 < v; j0 += 8, j1 += 8) { __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); __m256i naj = montgomery_add_256(T0, T1, m2, m0); __m256i najv = montgomery_sub_256(T0, T1, m2, m0); _mm256_storeu_si256((__m256i *)(a + j0), naj); _mm256_storeu_si256((__m256i *)(a + j1), najv); } } } int u = 1 << (2 + (k & 1)); int v = 1 << (k - 2 - (k & 1)); mint one = mint(1); mint imag = dw[1]; while (v) { if (v == 1) { mint ww = one, xx = one, wx = one; for (int jh = 0; jh < u;) { ww = xx * xx, wx = ww * xx; mint t0 = a[jh + 0], t1 = a[jh + 1] * xx; mint t2 = a[jh + 2] * ww, t3 = a[jh + 3] * wx; mint t0p2 = t0 + t2, t1p3 = t1 + t3; mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag; a[jh + 0] = t0p2 + t1p3, a[jh + 1] = t0p2 - t1p3; a[jh + 2] = t0m2 + t1m3, a[jh + 3] = t0m2 - t1m3; xx *= dw[__builtin_ctz((jh += 4))]; } } else if (v == 4) { const __m128i m0 = _mm_set1_epi32(0); const __m128i m1 = _mm_set1_epi32(mod); const __m128i m2 = _mm_set1_epi32(mod + mod); const __m128i r = _mm_set1_epi32(mint::r); const __m128i Imag = _mm_set1_epi32(imag.a); mint ww = one, xx = one, wx = one; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = j1 + v; int j3 = j2 + v; int je = v; for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i T0P2 = montgomery_add_128(T0, T2, m2, m0); const __m128i T1P3 = montgomery_add_128(T1, T3, m2, m0); const __m128i T0M2 = montgomery_sub_128(T0, T2, m2, m0); const __m128i T1M3 = montgomery_mul_128( montgomery_sub_128(T1, T3, m2, m0), Imag, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_sub_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j2), montgomery_add_128(T0M2, T1M3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j3), montgomery_sub_128(T0M2, T1M3, m2, m0)); } } else { ww = xx * xx, wx = ww * xx; const __m128i WW = _mm_set1_epi32(ww.a); const __m128i WX = _mm_set1_epi32(wx.a); const __m128i XX = _mm_set1_epi32(xx.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i MT1 = montgomery_mul_128(T1, XX, r, m1); const __m128i MT2 = montgomery_mul_128(T2, WW, r, m1); const __m128i MT3 = montgomery_mul_128(T3, WX, r, m1); const __m128i T0P2 = montgomery_add_128(T0, MT2, m2, m0); const __m128i T1P3 = montgomery_add_128(MT1, MT3, m2, m0); const __m128i T0M2 = montgomery_sub_128(T0, MT2, m2, m0); const __m128i T1M3 = montgomery_mul_128( montgomery_sub_128(MT1, MT3, m2, m0), Imag, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_sub_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j2), montgomery_add_128(T0M2, T1M3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j3), montgomery_sub_128(T0M2, T1M3, m2, m0)); } } xx *= dw[__builtin_ctz((jh += 4))]; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m1 = _mm256_set1_epi32(mod); const __m256i m2 = _mm256_set1_epi32(mod + mod); const __m256i r = _mm256_set1_epi32(mint::r); const __m256i Imag = _mm256_set1_epi32(imag.a); mint ww = one, xx = one, wx = one; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = j1 + v; int j3 = j2 + v; int je = v; for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i T0P2 = montgomery_add_256(T0, T2, m2, m0); const __m256i T1P3 = montgomery_add_256(T1, T3, m2, m0); const __m256i T0M2 = montgomery_sub_256(T0, T2, m2, m0); const __m256i T1M3 = montgomery_mul_256( montgomery_sub_256(T1, T3, m2, m0), Imag, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_sub_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j2), montgomery_add_256(T0M2, T1M3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j3), montgomery_sub_256(T0M2, T1M3, m2, m0)); } } else { ww = xx * xx, wx = ww * xx; const __m256i WW = _mm256_set1_epi32(ww.a); const __m256i WX = _mm256_set1_epi32(wx.a); const __m256i XX = _mm256_set1_epi32(xx.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i MT1 = montgomery_mul_256(T1, XX, r, m1); const __m256i MT2 = montgomery_mul_256(T2, WW, r, m1); const __m256i MT3 = montgomery_mul_256(T3, WX, r, m1); const __m256i T0P2 = montgomery_add_256(T0, MT2, m2, m0); const __m256i T1P3 = montgomery_add_256(MT1, MT3, m2, m0); const __m256i T0M2 = montgomery_sub_256(T0, MT2, m2, m0); const __m256i T1M3 = montgomery_mul_256( montgomery_sub_256(MT1, MT3, m2, m0), Imag, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_sub_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j2), montgomery_add_256(T0M2, T1M3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j3), montgomery_sub_256(T0M2, T1M3, m2, m0)); } } xx *= dw[__builtin_ctz((jh += 4))]; } } u <<= 2; v >>= 2; } } __attribute__((target("avx2"))) void intt(mint *a, int n, int normalize = true) { int k = n ? __builtin_ctz(n) : 0; if (k == 0) return; if (k == 1) { mint a1 = a[1]; a[1] = a[0] - a[1]; a[0] = a[0] + a1; if (normalize) { a[0] *= modinv(mint(2)); a[1] *= modinv(mint(2)); } return; } int u = 1 << (k - 2); int v = 1; mint one = mint(1); mint imag = dy[1]; while (u) { if (v == 1) { mint ww = one, xx = one, yy = one; u <<= 2; for (int jh = 0; jh < u;) { ww = xx * xx, yy = xx * imag; mint t0 = a[jh + 0], t1 = a[jh + 1]; mint t2 = a[jh + 2], t3 = a[jh + 3]; mint t0p1 = t0 + t1, t2p3 = t2 + t3; mint t0m1 = (t0 - t1) * xx, t2m3 = (t2 - t3) * yy; a[jh + 0] = t0p1 + t2p3, a[jh + 2] = (t0p1 - t2p3) * ww; a[jh + 1] = t0m1 + t2m3, a[jh + 3] = (t0m1 - t2m3) * ww; xx *= dy[__builtin_ctz(jh += 4)]; } } else if (v == 4) { const __m128i m0 = _mm_set1_epi32(0); const __m128i m1 = _mm_set1_epi32(mod); const __m128i m2 = _mm_set1_epi32(mod + mod); const __m128i r = _mm_set1_epi32(mint::r); const __m128i Imag = _mm_set1_epi32(imag.a); mint ww = one, xx = one, yy = one; u <<= 2; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = v + v; int j3 = j2 + v; for (; j0 < v; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i T0P1 = montgomery_add_128(T0, T1, m2, m0); const __m128i T2P3 = montgomery_add_128(T2, T3, m2, m0); const __m128i T0M1 = montgomery_sub_128(T0, T1, m2, m0); const __m128i T2M3 = montgomery_mul_128( montgomery_sub_128(T2, T3, m2, m0), Imag, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P1, T2P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j2), montgomery_sub_128(T0P1, T2P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_add_128(T0M1, T2M3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j3), montgomery_sub_128(T0M1, T2M3, m2, m0)); } } else { ww = xx * xx, yy = xx * imag; const __m128i WW = _mm_set1_epi32(ww.a); const __m128i XX = _mm_set1_epi32(xx.a); const __m128i YY = _mm_set1_epi32(yy.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i T0P1 = montgomery_add_128(T0, T1, m2, m0); const __m128i T2P3 = montgomery_add_128(T2, T3, m2, m0); const __m128i T0M1 = montgomery_mul_128( montgomery_sub_128(T0, T1, m2, m0), XX, r, m1); __m128i T2M3 = montgomery_mul_128( montgomery_sub_128(T2, T3, m2, m0), YY, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P1, T2P3, m2, m0)); _mm_storeu_si128( (__m128i *)(a + j2), montgomery_mul_128(montgomery_sub_128(T0P1, T2P3, m2, m0), WW, r, m1)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_add_128(T0M1, T2M3, m2, m0)); _mm_storeu_si128( (__m128i *)(a + j3), montgomery_mul_128(montgomery_sub_128(T0M1, T2M3, m2, m0), WW, r, m1)); } } xx *= dy[__builtin_ctz(jh += 4)]; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m1 = _mm256_set1_epi32(mod); const __m256i m2 = _mm256_set1_epi32(mod + mod); const __m256i r = _mm256_set1_epi32(mint::r); const __m256i Imag = _mm256_set1_epi32(imag.a); mint ww = one, xx = one, yy = one; u <<= 2; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = v + v; int j3 = j2 + v; for (; j0 < v; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i T0P1 = montgomery_add_256(T0, T1, m2, m0); const __m256i T2P3 = montgomery_add_256(T2, T3, m2, m0); const __m256i T0M1 = montgomery_sub_256(T0, T1, m2, m0); const __m256i T2M3 = montgomery_mul_256( montgomery_sub_256(T2, T3, m2, m0), Imag, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P1, T2P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j2), montgomery_sub_256(T0P1, T2P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_add_256(T0M1, T2M3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j3), montgomery_sub_256(T0M1, T2M3, m2, m0)); } } else { ww = xx * xx, yy = xx * imag; const __m256i WW = _mm256_set1_epi32(ww.a); const __m256i XX = _mm256_set1_epi32(xx.a); const __m256i YY = _mm256_set1_epi32(yy.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i T0P1 = montgomery_add_256(T0, T1, m2, m0); const __m256i T2P3 = montgomery_add_256(T2, T3, m2, m0); const __m256i T0M1 = montgomery_mul_256( montgomery_sub_256(T0, T1, m2, m0), XX, r, m1); const __m256i T2M3 = montgomery_mul_256( montgomery_sub_256(T2, T3, m2, m0), YY, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P1, T2P3, m2, m0)); _mm256_storeu_si256( (__m256i *)(a + j2), montgomery_mul_256(montgomery_sub_256(T0P1, T2P3, m2, m0), WW, r, m1)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_add_256(T0M1, T2M3, m2, m0)); _mm256_storeu_si256( (__m256i *)(a + j3), montgomery_mul_256(montgomery_sub_256(T0M1, T2M3, m2, m0), WW, r, m1)); } } xx *= dy[__builtin_ctz(jh += 4)]; } } u >>= 4; v <<= 2; } if (k & 1) { v = 1 << (k - 1); if (v < 8) { for (int j = 0; j < v; ++j) { mint ajv = a[j] - a[j + v]; a[j] += a[j + v]; a[j + v] = ajv; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m2 = _mm256_set1_epi32(mod + mod); int j0 = 0; int j1 = v; for (; j0 < v; j0 += 8, j1 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); __m256i naj = montgomery_add_256(T0, T1, m2, m0); __m256i najv = montgomery_sub_256(T0, T1, m2, m0); _mm256_storeu_si256((__m256i *)(a + j0), naj); _mm256_storeu_si256((__m256i *)(a + j1), najv); } } } if (normalize) { mint invn = modinv(mint(n)); for (int i = 0; i < n; i++) a[i] *= invn; } } __attribute__((target("avx2"))) void inplace_multiply( int l1, int l2, int zero_padding = true) { int l = l1 + l2 - 1; int M = 4; while (M < l) M <<= 1; if (zero_padding) { for (int i = l1; i < M; i++) ntt_inner::_buf1[i] = 0; for (int i = l2; i < M; i++) ntt_inner::_buf2[i] = 0; } const __m256i m0 = _mm256_set1_epi32(0); const __m256i m1 = _mm256_set1_epi32(mod); const __m256i r = _mm256_set1_epi32(mint::r); const __m256i N2 = _mm256_set1_epi32(mint::n2); for (int i = 0; i < l1; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf1 + i)); __m256i b = montgomery_mul_256(a, N2, r, m1); _mm256_storeu_si256((__m256i *)(ntt_inner::_buf1 + i), b); } for (int i = 0; i < l2; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf2 + i)); __m256i b = montgomery_mul_256(a, N2, r, m1); _mm256_storeu_si256((__m256i *)(ntt_inner::_buf2 + i), b); } ntt(buf1, M); ntt(buf2, M); for (int i = 0; i < M; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf1 + i)); __m256i b = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf2 + i)); __m256i c = montgomery_mul_256(a, b, r, m1); _mm256_storeu_si256((__m256i *)(ntt_inner::_buf1 + i), c); } intt(buf1, M, false); const __m256i INVM = _mm256_set1_epi32((mint(M).inverse()).a); for (int i = 0; i < l; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf1 + i)); __m256i b = montgomery_mul_256(a, INVM, r, m1); __m256i c = my256_mulhi_epu32(my256_mullo_epu32(b, r), m1); __m256i d = _mm256_and_si256(_mm256_cmpgt_epi32(c, m0), m1); __m256i e = _mm256_sub_epi32(d, c); _mm256_storeu_si256((__m256i *)(ntt_inner::_buf1 + i), e); } } void ntt(vector &a) { int M = (int)a.size(); for (int i = 0; i < M; i++) buf1[i].a = a[i].a; ntt(buf1, M); for (int i = 0; i < M; i++) a[i].a = buf1[i].a; } void intt(vector &a) { int M = (int)a.size(); for (int i = 0; i < M; i++) buf1[i].a = a[i].a; intt(buf1, M, true); for (int i = 0; i < M; i++) a[i].a = buf1[i].a; } vector multiply(const vector &a, const vector &b) { if (a.size() == 0 && b.size() == 0) return vector{}; int l = a.size() + b.size() - 1; if (min(a.size(), b.size()) <= 40) { vector s(l); for (int i = 0; i < (int)a.size(); ++i) for (int j = 0; j < (int)b.size(); ++j) s[i + j] += a[i] * b[j]; return s; } assert(l <= ntt_inner::SZ_FFT_BUF); int M = 4; while (M < l) M <<= 1; for (int i = 0; i < (int)a.size(); ++i) buf1[i].a = a[i].a; for (int i = (int)a.size(); i < M; ++i) buf1[i].a = 0; for (int i = 0; i < (int)b.size(); ++i) buf2[i].a = b[i].a; for (int i = (int)b.size(); i < M; ++i) buf2[i].a = 0; ntt(buf1, M); ntt(buf2, M); for (int i = 0; i < M; ++i) buf1[i].a = mint::reduce(uint64_t(buf1[i].a) * buf2[i].a); intt(buf1, M, false); vector s(l); mint invm = modinv(mint(M)); for (int i = 0; i < l; ++i) s[i] = buf1[i] * invm; return s; } void ntt_doubling(vector &a) { int M = (int)a.size(); for (int i = 0; i < M; i++) buf1[i].a = a[i].a; intt(buf1, M); mint r = 1, zeta = modpow(mint(pr),(mint::get_mod() - 1) / (M << 1)); for (int i = 0; i < M; i++) buf1[i] *= r, r *= zeta; ntt(buf1, M); a.resize(2 * M); for (int i = 0; i < M; i++) a[M + i].a = buf1[i].a; } }; // for garner static constexpr int m0 = 167772161; static constexpr int m1 = 469762049; static constexpr int m2 = 754974721; using mint0 = MontgomeryModInt; using mint1 = MontgomeryModInt; using mint2 = MontgomeryModInt; static constexpr int r01 = 104391568; static constexpr int r02 = 323560596; static constexpr int r12 = 399692502; static constexpr int r02r12 = 190329765; static constexpr i64 w1 = m0; static constexpr i64 w2 = i64(m0) * m1; using mint998 = MontgomeryModInt<998244353>; NumberTheoreticTransform ntt998; NumberTheoreticTransform ntt0; NumberTheoreticTransform ntt1; NumberTheoreticTransform ntt2; // small case (T = mint, long long) template vector naive_mul (const vector &A, const vector &B) { if (A.empty() || B.empty()) return {}; int N = (int)A.size(), M = (int)B.size(); vector res(N + M - 1); for (int i = 0; i < N; ++i) for (int j = 0; j < M; ++j) res[i + j] += A[i] * B[j]; return res; } // mint template vector mul(vector A,vector B) { if (A.empty() || B.empty()) return {}; int n = int(A.size()), m = int(B.size()); if (min(n, m) < 30) return naive_mul(A, B); int MOD = A[0].getmod(); if (MOD == 998244353) { vector a(n),b(m); for(int i=0;i c=ntt998.multiply(a,b); vector res(n+m-1); for(int i=0;i a0(n), b0(m); vector a1(n), b1(m); vector a2(n), b2(m); for (int i = 0; i < n; ++i) a0[i] = mint0(A[i].get()), a1[i] = mint1(A[i].get()), a2[i] = mint2(A[i].get()); for (int i = 0; i < m; ++i) b0[i] = mint0(B[i].get()), b1[i] = mint1(B[i].get()), b2[i] = mint2(B[i].get()); static const int W1 = w1%MOD, W2 = w2%MOD; vector c0=ntt0.multiply(a0,b0); vector c1=ntt1.multiply(a1,b1); vector c2=ntt2.multiply(a2,b2); vector res(n + m - 1); for (int i = 0; i < n + m - 1; ++i) { int n1 = c1[i].get(), n2 = c2[i].get(), a = c0[i].get(); int b = i64(n1 + m1 - a) * r01 % m1; int c = (i64(n2 + m2 - a) * r02r12 + i64(m2 - b) * r12) % m2; res[i] = mint(i64(a) + i64(b) * W1 + i64(c) * W2); } return res; } }; // Formal Power Series template struct FPS : vector { using vector::vector; /* template FPS(Args...args) : vector(args...){} */ // constructor FPS(const vector& r) : vector(r) {} // core operator inline FPS pre(int siz) const { return FPS(begin(*this), begin(*this) + min((int)this->size(), siz)); } inline FPS rev() const { FPS res = *this; reverse(begin(res), end(res)); return res; } inline FPS& normalize() { while (!this->empty() && this->back() == 0) this->pop_back(); return *this; } // basic operator inline FPS operator - () const noexcept { FPS res = (*this); for (int i = 0; i < (int)res.size(); ++i) res[i] = -res[i]; return res; } inline FPS operator + (const mint& v) const { return FPS(*this) += v; } inline FPS operator + (const FPS& r) const { return FPS(*this) += r; } inline FPS operator - (const mint& v) const { return FPS(*this) -= v; } inline FPS operator - (const FPS& r) const { return FPS(*this) -= r; } inline FPS operator * (const mint& v) const { return FPS(*this) *= v; } inline FPS operator * (const FPS& r) const { return FPS(*this) *= r; } inline FPS operator / (const mint& v) const { return FPS(*this) /= v; } inline FPS operator << (int x) const { return FPS(*this) <<= x; } inline FPS operator >> (int x) const { return FPS(*this) >>= x; } inline FPS& operator += (const mint& v) { if (this->empty()) this->resize(1); (*this)[0] += v; return *this; } inline FPS& operator += (const FPS& r) { if (r.size() > this->size()) this->resize(r.size()); for (int i = 0; i < (int)r.size(); ++i) (*this)[i] += r[i]; return this->normalize(); } inline FPS& operator -= (const mint& v) { if (this->empty()) this->resize(1); (*this)[0] -= v; return *this; } inline FPS& operator -= (const FPS& r) { if (r.size() > this->size()) this->resize(r.size()); for (int i = 0; i < (int)r.size(); ++i) (*this)[i] -= r[i]; return this->normalize(); } inline FPS& operator *= (const mint& v) { for (int i = 0; i < (int)this->size(); ++i) (*this)[i] *= v; return *this; } inline FPS& operator *= (const FPS& r) { return *this = NTT::ntt998.multiply((*this), r); } inline FPS& operator /= (const mint& v) { assert(v != 0); mint iv = modinv(v); for (int i = 0; i < (int)this->size(); ++i) (*this)[i] *= iv; return *this; } inline FPS& operator <<= (int x) { FPS res(x, 0); res.insert(res.end(), begin(*this), end(*this)); return *this = res; } inline FPS& operator >>= (int x) { FPS res; res.insert(res.end(), begin(*this) + x, end(*this)); return *this = res; } inline mint eval(const mint& v){ mint res = 0; for (int i = (int)this->size()-1; i >= 0; --i) { res *= v; res += (*this)[i]; } return res; } inline friend FPS gcd(const FPS& f, const FPS& g) { if (g.empty()) return f; return gcd(g, f % g); } // advanced operation // df/dx inline friend FPS diff(const FPS& f) { int n = (int)f.size(); FPS res(n-1); for (int i = 1; i < n; ++i) res[i-1] = f[i] * i; return res; } // \int f dx inline friend FPS integral(const FPS& f) { int n = (int)f.size(); FPS res(n+1, 0); for (int i = 0; i < n; ++i) res[i+1] = f[i] / (i+1); return res; } // inv(f), f[0] must not be 0 inline friend FPS inv(const FPS& f, int deg) { assert(f[0] != 0); if (deg < 0) deg = (int)f.size(); FPS res({mint(1) / f[0]}); for (int i = 1; i < deg; i <<= 1) { res = (res + res - res * res * f.pre(i << 1)).pre(i << 1); } res.resize(deg); return res; } inline friend FPS inv(const FPS& f) { return inv(f, f.size()); } // division, r must be normalized (r.back() must not be 0) inline FPS& operator /= (const FPS& r) { const int n=(*this).size(),m=r.size(); if(nnormalize(); if (this->size() < r.size()) { this->clear(); return *this; } int need = (int)this->size() - (int)r.size() + 1; *this = ((*this).rev().pre(need) * inv(r.rev(), need)).pre(need).rev(); return *this; } inline FPS& operator %= (const FPS &r) { const int n=(*this).size(),m=r.size(); if(nnormalize(); FPS q = (*this) / r; return *this -= q * r; } inline FPS operator / (const FPS& r) const { return FPS(*this) /= r; } inline FPS operator % (const FPS& r) const { return FPS(*this) %= r; } // log(f) = \int f'/f dx, f[0] must be 1 inline friend FPS log(const FPS& f, int deg) { assert(f[0] == 1); FPS res = integral((diff(f) * inv(f, deg)).pre(deg-1)); return res; } inline friend FPS log(const FPS& f) { return log(f, f.size()); } // exp(f), f[0] must be 0 inline friend FPS exp(const FPS& f, int deg) { assert(f[0] == 0); FPS res(1, 1); for (int i = 1; i < deg; i <<= 1) { res = res * (f.pre(i<<1) - log(res, i<<1) + 1).pre(i<<1); } res.resize(deg); return res; } inline friend FPS exp(const FPS& f) { return exp(f, f.size()); } // pow(f) = exp(e * log f) inline friend FPS pow(const FPS& f, long long e, int deg) { long long i = 0; while (i < (int)f.size() && f[i] == 0) ++i; if (i == (int)f.size()) return FPS(deg, 0); if (i * e >= deg) return FPS(deg, 0); mint k = f[i]; FPS res = exp(log((f >> i) / k, deg) * e, deg) * modpow(k, e) << (e * i); res.resize(deg); return res; } inline friend FPS pow(const FPS& f, long long e) { return pow(f, e, f.size()); } // sqrt(f), f[0] must be 1 inline friend FPS sqrt_base(const FPS& f, int deg) { assert(f[0] == 1); mint inv2 = mint(1) / 2; FPS res(1, 1); for (int i = 1; i < deg; i <<= 1) { res = (res + f.pre(i << 1) * inv(res, i << 1)).pre(i << 1); for (mint& x : res) x *= inv2; } res.resize(deg); return res; } inline friend FPS sqrt_base(const FPS& f) { return sqrt_base(f, f.size()); } FPS taylor_shift(mint c) const { int n = (int) this->size(); vector fact(n), rfact(n); fact[0] = rfact[0] = mint(1); for(int i = 1; i < n; i++) fact[i] = fact[i - 1] * mint(i); rfact[n - 1] = mint(1) / fact[n - 1]; for(int i = n - 1; i > 1; i--) rfact[i - 1] = rfact[i] * mint(i); FPS p(*this); for(int i = 0; i < n; i++) p[i] *= fact[i]; p = p.rev(); FPS bs(n, mint(1)); for(int i = 1; i < n; i++) bs[i] = bs[i - 1] * c * rfact[i] * fact[i - 1]; p = (p * bs).pre(n); p = p.rev(); for(int i = 0; i < n; i++) p[i] *= rfact[i]; return p; } }; template FPS product(vector> a){ int siz=1; while(siz> res(siz*2-1,{1}); for(int i=0;i=0;--i){ res[i]=res[2*i+1]*res[2*i+2]; } return res[0]; } template FPS inv_sum(vector> f){ int siz=1; while(siz> mol(siz*2-1),dem(siz*2-1,{1}); for(size_t i=0;i=0;--i){ dem[i]=dem[2*i+1]*dem[2*i+2]; mol[i]=mol[2*i+1]*dem[2*i+2]+mol[2*i+2]*dem[2*i+1]; } mol[0]*=inv(dem[0]); return mol[0]; } template FPS rev(FPS p) { reverse(p.begin(),p.end()); return p; } template FPS RSZ(FPS p, int x) { p.resize(x); return p; } //fast Input by yosupo #include #include #include #include #include #include #include #include #include #include namespace fastio{ /* quote from yosupo's submission in Library Checker */ int bsr(unsigned int n) { return 8 * (int)sizeof(unsigned int) - 1 - __builtin_clz(n); } // @param n `1 <= n` // @return maximum non-negative `x` s.t. `(n & (1 << x)) != 0` int bsr(unsigned long n) { return 8 * (int)sizeof(unsigned long) - 1 - __builtin_clzl(n); } // @param n `1 <= n` // @return maximum non-negative `x` s.t. `(n & (1 << x)) != 0` int bsr(unsigned long long n) { return 8 * (int)sizeof(unsigned long long) - 1 - __builtin_clzll(n); } // @param n `1 <= n` // @return minimum non-negative `x` s.t. `(n & (1 << x)) != 0` int bsr(unsigned __int128 n) { unsigned long long low = (unsigned long long)(n); unsigned long long high = (unsigned long long)(n >> 64); return high ? 127 - __builtin_clzll(high) : 63 - __builtin_ctzll(low); } namespace internal { template using is_signed_int128 = typename std::conditional::value || std::is_same::value, std::true_type, std::false_type>::type; template using is_unsigned_int128 = typename std::conditional::value || std::is_same::value, std::true_type, std::false_type>::type; template using make_unsigned_int128 = typename std::conditional::value, __uint128_t, unsigned __int128>; template using is_integral = typename std::conditional::value || internal::is_signed_int128::value || internal::is_unsigned_int128::value, std::true_type, std::false_type>::type; template using is_signed_int = typename std::conditional<(is_integral::value && std::is_signed::value) || is_signed_int128::value, std::true_type, std::false_type>::type; template using is_unsigned_int = typename std::conditional<(is_integral::value && std::is_unsigned::value) || is_unsigned_int128::value, std::true_type, std::false_type>::type; template using to_unsigned = typename std::conditional< is_signed_int128::value, make_unsigned_int128, typename std::conditional::value, std::make_unsigned, std::common_type>::type>::type; template using is_integral_t = std::enable_if_t::value>; template using is_signed_int_t = std::enable_if_t::value>; template using is_unsigned_int_t = std::enable_if_t::value>; template using to_unsigned_t = typename to_unsigned::type; } // namespace internal struct Scanner { public: Scanner(const Scanner&) = delete; Scanner& operator=(const Scanner&) = delete; Scanner(FILE* fp) : fd(fileno(fp)) {} void read() {} template void read(H& h, T&... t) { bool f = read_single(h); assert(f); read(t...); } int read_unsafe() { return 0; } template int read_unsafe(H& h, T&... t) { bool f = read_single(h); if (!f) return 0; return 1 + read_unsafe(t...); } int close() { return ::close(fd); } private: static constexpr int SIZE = 1 << 15; int fd = -1; std::array line; int st = 0, ed = 0; bool eof = false; bool read_single(std::string& ref) { if (!skip_space()) return false; ref = ""; while (true) { char c = top(); if (c <= ' ') break; ref += c; st++; } return true; } bool read_single(double& ref) { std::string s; if (!read_single(s)) return false; ref = std::stod(s); return true; } template ::value>* = nullptr> bool read_single(T& ref) { if (!skip_space<50>()) return false; ref = top(); st++; return true; } template * = nullptr, std::enable_if_t::value>* = nullptr> bool read_single(T& sref) { using U = internal::to_unsigned_t; if (!skip_space<50>()) return false; bool neg = false; if (line[st] == '-') { neg = true; st++; } U ref = 0; do { ref = 10 * ref + (line[st++] & 0x0f); } while (line[st] >= '0'); sref = neg ? -ref : ref; return true; } template * = nullptr, std::enable_if_t::value>* = nullptr> bool read_single(U& ref) { if (!skip_space<50>()) return false; ref = 0; do { ref = 10 * ref + (line[st++] & 0x0f); } while (line[st] >= '0'); return true; } bool reread() { if (ed - st >= 50) return true; if (st > SIZE / 2) { std::memmove(line.data(), line.data() + st, ed - st); ed -= st; st = 0; } if (eof) return false; auto u = ::read(fd, line.data() + ed, SIZE - ed); if (u == 0) { eof = true; line[ed] = '\0'; u = 1; } ed += int(u); line[ed] = char(127); return true; } char top() { if (st == ed) { bool f = reread(); assert(f); } return line[st]; } template bool skip_space() { while (true) { while (line[st] <= ' ') st++; if (ed - st > TOKEN_LEN) return true; if (st > ed) st = ed; for (auto i = st; i < ed; i++) { if (line[i] <= ' ') return true; } if (!reread()) return false; } } }; //fast Output by ei1333 /** * @brief Printer(高速出力) */ struct Printer { public: explicit Printer(FILE *fp) : fp(fp) {} ~Printer() { flush(); } template< bool f = false, typename T, typename... E > void write(const T &t, const E &... e) { if(f) write_single(' '); write_single(t); write< true >(e...); } template< typename... T > void writeln(const T &...t) { write(t...); write_single('\n'); } void flush() { fwrite(line, 1, st - line, fp); st = line; } private: FILE *fp = nullptr; static constexpr size_t line_size = 1 << 16; static constexpr size_t int_digits = 20; char line[line_size + 1] = {}; char small[32] = {}; char *st = line; template< bool f = false > void write() {} void write_single(const char &t) { if(st + 1 >= line + line_size) flush(); *st++ = t; } template< typename T, enable_if_t< is_integral< T >::value, int > = 0 > void write_single(T s) { if(st + int_digits >= line + line_size) flush(); if(s == 0) { write_single('0'); return; } if(s < 0) { write_single('-'); s = -s; } char *mp = small + sizeof(small); typename make_unsigned< T >::type y = s; size_t len = 0; while(y > 0) { *--mp = y % 10 + '0'; y /= 10; ++len; } memmove(st, mp, len); st += len; } void write_single(const string &s) { for(auto &c : s) write_single(c); } void write_single(const char *s) { while(*s != 0) write_single(*s++); } template< typename T > void write_single(const vector< T > &s) { for(size_t i = 0; i < s.size(); i++) { if(i) write_single(' '); write_single(s[i]); } } }; }; //namespace fastio using mint=MontgomeryModInt<998244353>; int main(){ fastio::Scanner sc(stdin); fastio::Printer pr(stdout); #define in(...) sc.read(__VA_ARGS__) #define LL(...) ll __VA_ARGS__;in(__VA_ARGS__) #define INT(...) int __VA_ARGS__;in(__VA_ARGS__) #define STR(...) string __VA_ARGS__;in(__VA_ARGS__) #define out(...) pr.write(__VA_ARGS__) #define outln(...) pr.writeln(__VA_ARGS__) #define outspace(...) pr.write(__VA_ARGS__);pr.write(' ') #define rall(v) (v).rbegin(), (v).rend() #define fi first #define se second /* 草 好き */ INT(t); while(t--){ LL(n,x); if(x==1){ mint ans=n/2+(n%2?1:0); outln(ans.get()); }else{ mint ans=n/(x+1)*((x+1)/2+((x+1)%2?1:0)); outln(ans.get()); } } }