#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; } template struct subproduct_tree{ using poly=FPS; vector tree; int n,siz; subproduct_tree(const vector &x){ n=1; siz=sz(x); while(n0;i--) tree[i]=tree[2*i]*tree[2*i+1]; } vector multieval(const poly &f){ vector remainder(2*n); remainder[1]=f%tree[1]; for(int i=1;i ret(siz); for(int i=0;i &y){ poly g=diff(tree[1]); vector evaled=multieval(g); vector mol(2*n),dem(2*n,{1}); for(int i=0;i0;--i){ dem[i]=dem[2*i]*dem[2*i+1]; mol[i]=mol[2*i]*dem[2*i+1]+mol[2*i+1]*dem[2*i]; } mol[1]*=inv(dem[1]); return RSZ(tree[1]*mol[1],siz); } }; template vector multieval(const FPS &f,const vector &x){ subproduct_tree tree(x); return tree.multieval(f); } template FPS interpolate(const vector &x,const vector &y){ assert(sz(x)==sz(y)); if(sz(x)==1) return {y[0]}; subproduct_tree tree(x); return tree.interpolate(y); } template mint Bostan_Mori(ll n,FPS P,FPS Q){ while(n){ auto C=Q; for(int i=1;i H; for(int i=(n&1ll);i L; for(int i=0;i>=1; } return P[0]/Q[0]; } template< class T > struct Matrix { vector< vector< T > > A; Matrix() {} Matrix(size_t n, size_t m) : A(n, vector< T >(m, 0)) {} Matrix(size_t n) : A(n, vector< T >(n, 0)) {}; size_t size() const { if(A.empty()) return 0; assert(A.size() == A[0].size()); return A.size(); } size_t height() const { return (A.size()); } size_t width() const { return (A[0].size()); } inline const vector< T > &operator[](int k) const { return (A.at(k)); } inline vector< T > &operator[](int k) { return (A.at(k)); } static Matrix I(size_t n) { Matrix mat(n); for(int i = 0; i < n; i++) mat[i][i] = 1; return (mat); } Matrix &operator+=(const Matrix &B) { size_t n = height(), m = width(); assert(n == B.height() && m == B.width()); for(int i = 0; i < n; i++) for(int j = 0; j < m; j++) (*this)[i][j] += B[i][j]; return (*this); } Matrix &operator-=(const Matrix &B) { size_t n = height(), m = width(); assert(n == B.height() && m == B.width()); for(int i = 0; i < n; i++) for(int j = 0; j < m; j++) (*this)[i][j] -= B[i][j]; return (*this); } Matrix &operator*=(const Matrix &B) { size_t n = height(), m = B.width(), p = width(); assert(p == B.height()); vector< vector< T > > C(n, vector< T >(m, 0)); for(int i = 0; i < n; i++) for(int j = 0; j < m; j++) for(int k = 0; k < p; k++) C[i][j] = (C[i][j] + (*this)[i][k] * B[k][j]); A.swap(C); return (*this); } Matrix &operator^=(long long k) { Matrix B = Matrix::I(height()); while(k > 0) { if(k & 1) B *= *this; *this *= *this; k >>= 1LL; } A.swap(B.A); return (*this); } Matrix operator+(const Matrix &B) const { return (Matrix(*this) += B); } Matrix operator-(const Matrix &B) const { return (Matrix(*this) -= B); } Matrix operator*(const Matrix &B) const { return (Matrix(*this) *= B); } friend Matrix pow(const Matrix& B,const long long k) { return (Matrix(B) ^= k); } friend ostream &operator<<(ostream &os, Matrix &p) { size_t n = p.height(), m = p.width(); for(int i = 0; i < n; i++) { os << "["; for(int j = 0; j < m; j++) { os << p[i][j] << (j + 1 == m ? "]\n" : ","); } } return (os); } T determinant() { Matrix B(*this); assert(width() == height()); T ret = 1; for(int i = 0; i < width(); i++) { int idx = -1; for(int j = i; j < width(); j++) { if(B[j][i] != 0) idx = j; } if(idx == -1) return (0); if(i != idx) { ret *= -1; swap(B[i], B[idx]); } ret *= B[i][i]; T vv = B[i][i]; for(int j = 0; j < width(); j++) { B[i][j] /= vv; } for(int j = i + 1; j < width(); j++) { T a = B[j][i]; for(int k = 0; k < width(); k++) { B[j][k] -= B[i][k] * a; } } } return (ret); } }; 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(n,m); ll n; cin >> n; Matrix o(1,2),m(2,2); m[0][0]=m[0][1]=m[1][0]=1; m=pow(m,n-1); o[0][0]=1; o*=m; cout << (o[0][0]+o[0][1]-1).get() << endl; }