結果
問題 | No.2215 Slide Subset Sum |
ユーザー | kaichou243 |
提出日時 | 2023-01-28 13:22:14 |
言語 | C++17 (gcc 12.3.0 + boost 1.83.0) |
結果 |
WA
|
実行時間 | - |
コード長 | 59,605 bytes |
コンパイル時間 | 3,560 ms |
コンパイル使用メモリ | 321,864 KB |
実行使用メモリ | 91,352 KB |
最終ジャッジ日時 | 2024-07-04 06:48:53 |
合計ジャッジ時間 | 16,938 ms |
ジャッジサーバーID (参考情報) |
judge5 / judge4 |
(要ログイン)
テストケース
テストケース表示入力 | 結果 | 実行時間 実行使用メモリ |
---|---|---|
testcase_00 | WA | - |
testcase_01 | WA | - |
testcase_02 | WA | - |
testcase_03 | WA | - |
testcase_04 | WA | - |
testcase_05 | AC | 2 ms
6,944 KB |
testcase_06 | AC | 3 ms
6,940 KB |
testcase_07 | AC | 2 ms
6,940 KB |
testcase_08 | AC | 2 ms
6,944 KB |
testcase_09 | AC | 2 ms
6,944 KB |
testcase_10 | AC | 2 ms
6,944 KB |
testcase_11 | AC | 2 ms
6,940 KB |
testcase_12 | AC | 3 ms
6,940 KB |
testcase_13 | AC | 2 ms
6,940 KB |
testcase_14 | AC | 3 ms
6,940 KB |
testcase_15 | AC | 2 ms
6,944 KB |
testcase_16 | AC | 3 ms
6,940 KB |
testcase_17 | WA | - |
testcase_18 | WA | - |
testcase_19 | WA | - |
testcase_20 | AC | 489 ms
51,048 KB |
testcase_21 | WA | - |
testcase_22 | AC | 509 ms
54,484 KB |
testcase_23 | AC | 612 ms
89,328 KB |
testcase_24 | AC | 613 ms
90,556 KB |
testcase_25 | AC | 584 ms
81,488 KB |
testcase_26 | AC | 506 ms
54,792 KB |
testcase_27 | WA | - |
testcase_28 | AC | 28 ms
7,052 KB |
testcase_29 | WA | - |
testcase_30 | WA | - |
testcase_31 | AC | 264 ms
35,784 KB |
testcase_32 | WA | - |
testcase_33 | WA | - |
testcase_34 | WA | - |
testcase_35 | WA | - |
testcase_36 | AC | 135 ms
17,080 KB |
testcase_37 | AC | 355 ms
90,876 KB |
testcase_38 | AC | 59 ms
15,992 KB |
testcase_39 | WA | - |
testcase_40 | WA | - |
testcase_41 | AC | 2 ms
6,944 KB |
testcase_42 | AC | 487 ms
50,924 KB |
testcase_43 | WA | - |
testcase_44 | WA | - |
testcase_45 | AC | 566 ms
76,840 KB |
testcase_46 | AC | 561 ms
74,240 KB |
ソースコード
#include<bits/stdc++.h> #include <immintrin.h> #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<ll,ll>; const long double PI=acos(-1); const ll INF=1e18; const int inf=1e9; struct Edge { ll to; ll cost; }; using Graph=vector<vector<Edge>>; template <typename T> bool chmax(T &a,const T& b){ if (a<b){ a=b; return true; } return false; } template <typename T> bool chmin(T &a,const T& b){ if (a>b){ a=b; return true; } return false; } template<int MOD> 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<MOD>& 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<MOD>& x) noexcept { return os << x.val; } friend constexpr Fp<MOD> modpow(const Fp<MOD>& a, long long n) noexcept { Fp<MOD> res=1,r=a; while(n){ if(n&1) res*=r; r*=r; n>>=1; } return res; } friend constexpr Fp<MOD> modinv(const Fp<MOD>& 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<MOD>(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<mod> modpow(const MontgomeryModInt<mod> &x,u64 n) noexcept { MontgomeryModInt<mod> ret(1), mul(x); while(n > 0) { if(n & 1) ret *= mul; mul *= mul; n >>= 1; } return ret; } friend constexpr MontgomeryModInt<mod> modinv(const MontgomeryModInt<mod> &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<mod>(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<class T,T (*op)(T,T),T (*e)()> struct SegmentTree{ int n; vector<T> dat; SegmentTree(int N){ n=1; while(n<N)n*=2; dat.assign(2*n,e()); } void add(int k,T x){ k+=n; dat[k]+=x; while(k){ k>>=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<r){ if(l&1) prodl=op(prodl,dat[l++]); if(r&1) prodr=op(dat[--r],prodr); l>>=1; r>>=1; } return op(prodl,prodr); } }; struct FenwickTree{ int n; vector<ll> dat; FenwickTree(int N){ n=1; while(n<N) n*=2; dat.assign(n,0); } void add(int k,ll x){ k+=1; while(k<=n){ dat[k]+=x; k+=k&-k; } } //sum[0,k) ll sum(int k){ ll ans=0; while(k){ ans+=dat[k]; k-=k&-k; } return ans; } //sum[l,r) ll sum(int l,int r){ return sum(r)-sum(l); } }; template<class S,S (*op)(S,S),S (*e)(),class F,S (*mapping)(F,S),F (*composition)(F,F),F (*id)()> struct LazySegTree{ private: int _n,size=1,idx=0; vector<S>seq; vector<F>lazy; 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<size)lazy[k]=composition(f,lazy[k]); } void eval(int k){ all_apply(2*k,lazy[k]); all_apply(2*k+1,lazy[k]); lazy[k]=id(); } public: LazySegTree(int n):LazySegTree(vector<S>(n,e())){} LazySegTree(const vector<S>&v):_n(int(v.size())){ while(size<_n)size<<=1,idx++; seq=vector<S>(2*size,e()); lazy=vector<F>(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)!=l)eval(l>>i); if(((r>>i)<<i)!=r)eval(r>>i); } while(l<r){ if(l&1)sml=op(sml,seq[l++]); if(r&1)smr=op(seq[--r],smr); 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)!=l)eval(l>>i); if(((r>>i)<<i)!=r)eval((r-1)>>i); } int l2=l,r2=r; while(l<r){ if(l&1)all_apply(l++,f); if(r&1)all_apply(--r,f); l>>=1; r>>=1; } l=l2,r=r2; for(int i=1;i<=idx;i++){ if(((l>>i)<<i)!=l)update(l>>i); if(((r>>i)<<i)!=r)update((r-1)>>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<P> prime_factorize(ll N) { vector<P> 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<ll> &b, const vector<ll> &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 <class mint> struct fft_info{ static constexpr int rank2 = bsf_constexpr(mint::getmod() - 1); std::array<mint, rank2 + 1> root; // root[i]^(2^i) == 1 std::array<mint, rank2 + 1> iroot; // root[i] * iroot[i] == 1 std::array<mint, std::max(0, rank2 - 2 + 1)> rate2; std::array<mint, std::max(0, rank2 - 2 + 1)> irate2; std::array<mint, std::max(0, rank2 - 3 + 1)> rate3; std::array<mint, std::max(0, rank2 - 3 + 1)> 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 <class mint> void trans(std::vector<mint>& a) { int n = int(a.size()); int h = ceil_pow2(n); int MOD=a[0].getmod(); static const fft_info<mint> 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 <class mint> void trans_inv(std::vector<mint>& a) { int n = int(a.size()); int h = ceil_pow2(n); static const fft_info<mint> 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 <typename mint> 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<mint> &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<mint> &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<mint> multiply(const vector<mint> &a, const vector<mint> &b) { if (a.size() == 0 && b.size() == 0) return vector<mint>{}; int l = a.size() + b.size() - 1; if (min<int>(a.size(), b.size()) <= 40) { vector<mint> 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<mint> 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<mint> &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<m0>; using mint1 = MontgomeryModInt<m1>; using mint2 = MontgomeryModInt<m2>; 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<mint998> ntt998; NumberTheoreticTransform<mint0> ntt0; NumberTheoreticTransform<mint1> ntt1; NumberTheoreticTransform<mint2> ntt2; // small case (T = mint, long long) template<class T> vector<T> naive_mul (const vector<T> &A, const vector<T> &B) { if (A.empty() || B.empty()) return {}; int N = (int)A.size(), M = (int)B.size(); vector<T> 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<class mint> vector<mint> mul(vector<mint> A,vector<mint> 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<mint998> a(n),b(m); for(int i=0;i<n;i++) a[i]=mint998(A[i].get()); for(int i=0;i<m;i++) b[i]=mint998(B[i].get()); vector<mint998> c=ntt998.multiply(a,b); vector<mint> res(n+m-1); for(int i=0;i<n+m-1;i++) res[i]=c[i].get(); return res; } vector<mint0> a0(n), b0(m); vector<mint1> a1(n), b1(m); vector<mint2> 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<mint0> c0=ntt0.multiply(a0,b0); vector<mint1> c1=ntt1.multiply(a1,b1); vector<mint2> c2=ntt2.multiply(a2,b2); vector<mint> 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 <typename mint> struct FPS : vector<mint> { using vector<mint>::vector; /* template<class...Args> FPS(Args...args) : vector<mint>(args...){} */ // constructor FPS(const vector<mint>& r) : vector<mint>(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(n<m){ (*this).clear(); return *this; } assert(r.back() != 0); this->normalize(); 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(n<m) return (*this); assert(r.back() != 0); this->normalize(); 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<mint> 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 <typename mint> FPS<mint> product(vector<FPS<mint>> a){ int siz=1; while(siz<int(a.size())){ siz<<=1; } vector<FPS<mint>> res(siz*2-1,{1}); for(int i=0;i<int(a.size());++i){ res[i+siz-1]=a[i]; } for(int i=siz-2;i>=0;--i){ res[i]=res[2*i+1]*res[2*i+2]; } return res[0]; } template<typename mint> FPS<mint> inv_sum(vector<FPS<mint>> f){ int siz=1; while(siz<int(f.size())){ siz<<=1; } vector<FPS<mint>> mol(siz*2-1),dem(siz*2-1,{1}); for(size_t i=0;i<f.size();++i){ mol[i+siz-1]={1}; dem[i+siz-1]=f[i]; } for(int i=siz-2;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 <typename mint> FPS<mint> rev(FPS<mint> p) { reverse(p.begin(),p.end()); return p; } template <typename mint> FPS<mint> RSZ(FPS<mint> p, int x) { p.resize(x); return p; } template<typename mint> struct subproduct_tree{ using poly=FPS<mint>; vector<poly> tree; int n,siz; subproduct_tree(const vector<mint> &x){ n=1; siz=sz(x); while(n<siz) n*=2;; tree.resize(2*n,{mint(1)}); for(int i=0;i<siz;i++) tree[i+n]={-x[i],mint(1)}; for(int i=n-1;i>0;i--) tree[i]=tree[2*i]*tree[2*i+1]; } vector<mint> multieval(const poly &f){ vector<poly> remainder(2*n); remainder[1]=f%tree[1]; for(int i=1;i<n;i++){ remainder[2*i]=remainder[i]%tree[2*i]; remainder[2*i+1]=remainder[i]%tree[2*i+1]; } vector<mint> ret(siz); for(int i=0;i<siz;i++){ if(remainder[i+n].empty()) ret[i]=0; else ret[i]=remainder[i+n][0]; } return ret; } poly interpolate(const vector<mint> &y){ poly g=diff(tree[1]); vector<mint> evaled=multieval(g); vector<poly> mol(2*n),dem(2*n,{1}); for(int i=0;i<siz;++i){ mol[i+n]={y[i]}; dem[i+n]=tree[i+n]*evaled[i]; } for(int i=n-1;i>0;--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 <typename mint> vector<mint> multieval(const FPS<mint> &f,const vector<mint> &x){ subproduct_tree<mint> tree(x); return tree.multieval(f); } template <typename mint> FPS<mint> interpolate(const vector<mint> &x,const vector<mint> &y){ assert(sz(x)==sz(y)); if(sz(x)==1) return {y[0]}; subproduct_tree<mint> tree(x); return tree.interpolate(y); } template<typename mint> mint Bostan_Mori(ll n,FPS<mint> P,FPS<mint> Q){ while(n){ auto C=Q; for(int i=1;i<C.size();i+=2) C[i]*=-1; P*=C; Q*=C; FPS<mint> H; for(int i=(n&1ll);i<P.size();i+=2) H.push_back(P[i]); P=H; FPS<mint> L; for(int i=0;i<Q.size();i+=2) L.push_back(Q[i]); Q=L; n>>=1; } return P[0]/Q[0]; } using mint=MontgomeryModInt<998244353>; int main(){ #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,k; cin >> n >> m >> k; assert(1<=n&&n<=200000&&1<=m&&m<=n&&1<=k&&k<=100); vector<int> a(n); for(int i=0;i<n;i++){ cin >> a[i]; assert(0<=a[i]&&a[i]<k); } vector<vector<mint>> dp1,dp2; vector<int> s1,s2; vector<mint> v(k); v[0]=1; dp1.push_back(v); dp2.push_back(v); for(int i=0;i<m;i++){ vector<mint> ndp(k); for(int j=0;j<k;j++){ ndp[(j+a[i])%k]=dp1[sz(dp1)-1][(j+a[i])%k]+dp1[sz(dp1)-1][j]; } dp1.push_back(ndp); s1.push_back(a[i]); } cout << (dp1[sz(dp1)-1][0]-1).get() << '\n'; for(int i=0;i<n-m;i++){ if(s2.size()==0){ while(s1.size()){ int x=s1.back(); s1.pop_back(); if(s1.size()>0) s2.push_back(x); dp1.pop_back(); } for(int j=0;j<s2.size();j++){ vector<mint> ndp(k); for(int t=0;t<k;t++){ ndp[(t+s2[j])%k]=dp2[sz(dp2)-1][(t+s2[j])%k]+dp2[sz(dp2)-1][t]; } dp2.push_back(ndp); } }else{ s2.pop_back(); dp2.pop_back(); } vector<mint> ndp(k); for(int j=0;j<k;j++){ ndp[(j+a[i+m])%k]=dp1[sz(dp1)-1][(j+a[i+m])%k]+dp1[sz(dp1)-1][j]; } dp1.push_back(ndp); mint ans=dp1[sz(dp1)-1][0]*dp2[sz(dp2)-1][0]-1; for(int j=1;j<k;j++){ ans+=dp1[sz(dp1)-1][j]*dp2[sz(dp2)-1][k-j]; } cout << ans.get() << '\n'; } }