#ifdef NACHIA #define _GLIBCXX_DEBUG #else #define NDEBUG #endif #include #include #include #include using i64 = long long; using u64 = unsigned long long; #define rep(i,n) for(int i=0; i void chmin(A& l, const A& r){ if(r < l) l = r; } template void chmax(A& l, const A& r){ if(l < r) l = r; } #include using Modint = atcoder::static_modint<998244353>; #include namespace nachia{ template struct PrimitiveRoot{ using u64 = unsigned long long; static constexpr u64 powm(u64 a, u64 i) { u64 res = 1, aa = a; for( ; i; i /= 2){ if(i & 1) res = res * aa % MOD; aa = aa * aa % MOD; } return res; } static constexpr bool ExamineVal(unsigned int g){ u64 t = MOD - 1; for(u64 d=2; d*d<=t; d+=1+(d&1)) if(t % d == 0){ if(powm(g, (MOD - 1) / d) == 1) return false; while(t % d == 0) t /= d; } if(t != 1) if(powm(g, (MOD - 1) / t) == 1) return false; return true; } static constexpr unsigned int GetVal(){ for(u64 x=2; x class Comb{ private: std::vector F; std::vector iF; public: void extend(int newN){ int prevN = (int)F.size() - 1; if(prevN >= newN) return; F.resize(newN+1); iF.resize(newN+1); for(int i=prevN+1; i<=newN; i++) F[i] = F[i-1] * Modint::raw(i); iF[newN] = F[newN].inv(); for(int i=newN; i>prevN; i--) iF[i-1] = iF[i] * Modint::raw(i); } Comb(int n = 1){ F.assign(2, Modint(1)); iF.assign(2, Modint(1)); extend(n); } Modint factorial(int n) const { return F[n]; } Modint invFactorial(int n) const { return iF[n]; } Modint invOf(int n) const { return iF[n] * F[n-1]; } Modint comb(int n, int r) const { if(n < 0 || n < r || r < 0) return Modint(0); return F[n] * iF[r] * iF[n-r]; } Modint invComb(int n, int r) const { if(n < 0 || n < r || r < 0) return Modint(0); return iF[n] * F[r] * F[n-r]; } Modint perm(int n, int r) const { if(n < 0 || n < r || r < 0) return Modint(0); return F[n] * iF[n-r]; } Modint invPerm(int n, int r) const { if(n < 0 || n < r || r < 0) return Modint(0); return iF[n] * F[n-r]; } Modint operator()(int n, int r) const { return comb(n,r); } }; } // namespace nachia namespace nachia{ int Popcount(unsigned long long c) noexcept { #ifdef __GNUC__ return __builtin_popcountll(c); #else c = (c & (~0ull/3)) + ((c >> 1) & (~0ull/3)); c = (c & (~0ull/5)) + ((c >> 2) & (~0ull/5)); c = (c & (~0ull/17)) + ((c >> 4) & (~0ull/17)); c = (c * (~0ull/257)) >> 56; return c; #endif } // please ensure x != 0 int MsbIndex(unsigned long long x) noexcept { #ifdef __GNUC__ return 63 - __builtin_clzll(x); #else using u64 = unsigned long long; int q = (x >> 32) ? 32 : 0; auto m = x >> q; constexpr u64 hi = 0x88888888; constexpr u64 mi = 0x11111111; m = (((m | ~(hi - (m & ~hi))) & hi) * mi) >> 35; m = (((m | ~(hi - (x & ~hi))) & hi) * mi) >> 31; q += (m & 0xf) << 2; q += 0x3333333322221100 >> (((x >> q) & 0xf) << 2) & 0xf; return q; #endif } // please ensure x != 0 int LsbIndex(unsigned long long x) noexcept { #ifdef __GNUC__ return __builtin_ctzll(x); #else return MsbIndex(x & -x); #endif } } namespace nachia { template struct NttInterface{ template void Butterfly(Iter, int) const {} template void IButterfly(Iter, int) const {} template void BitReversal(Iter a, int N) const { for(int i=0, j=0; j>1; k > (i^=k); k>>=1); } } }; } // namespace nachia #include #include namespace nachia{ template struct Ntt : NttInterface { using u32 = unsigned int; using u64 = unsigned long long; static int ceil_pow2(int n) { int x = 0; while ((1U << x) < (u32)(n)) x++; return x; } static constexpr int bsf_constexpr(unsigned int n) { int x = 0; while (!(n & (1 << x))) x++; return x; } struct fft_info { static constexpr u32 g = nachia::PrimitiveRoot::val; static constexpr int rank2 = bsf_constexpr(mint::mod()-1); using RootTable = std::array; RootTable root, iroot, rate3, irate3; fft_info(){ root[rank2] = mint(g).pow((mint::mod() - 1) >> rank2); iroot[rank2] = root[rank2].inv(); 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-3; i++){ rate3[i] = root[i+3] * prod; irate3[i] = iroot[i+3] * iprod; prod *= iroot[i+3]; iprod *= root[i+3]; } } }; template void ButterflyLayered(RandomAccessIterator a, int n, int stride, int repeat) const { static const fft_info info; int h = n * stride; while(repeat--){ int len = 1; int p = h; if(ceil_pow2(n)%2 == 1){ p >>= 1; for(int i=0; i stride; ){ p >>= 2; mint rot = 1, imag = info.root[2]; u64 mod2 = u64(mint::mod()) * mint::mod(); int offset = p; for(int s=0; s void Butterfly(RandomAccessIterator a, int n) const { ButterflyLayered(a, n, 1, 1); } template void IButterflyLayered(RandomAccessIterator a, int n, int stride, int repeat) const { static const fft_info info; constexpr int MOD = mint::mod(); while(repeat--){ int len = n; int p = stride; for( ; 2 < len; ){ len >>= 2; mint irot = 1, iimag = info.iroot[2]; int offset = p; for(int s=0; s void IButterfly(RandomAccessIterator a, int n) const { IButterflyLayered(a, n, 1, 1); } }; } // namespace nachia namespace nachia { template> struct FpsNtt { public: using Fps = FpsNtt; using ElemTy = Elem; static constexpr unsigned int MOD = Elem::mod(); static constexpr int CONV_THRES = 30; static const NttInst nttInst; static const unsigned int zeta = nachia::PrimitiveRoot::GetVal(); private: using u32 = unsigned int; static Elem ZeroElem() noexcept { return Elem(0); } static Elem OneElem() noexcept { return Elem(1); } static Comb comb; std::vector a; int RSZ(int& sz) const { return sz = (sz < 0 ? size() : sz); } public: int size() const noexcept { return a.size(); } Elem& operator[](int x) noexcept { return a[x]; } const Elem& operator[](int x) const noexcept { return a[x]; } Elem getCoeff(int x) const noexcept { return (0 <= x && x < size()) ? a[x] : ZeroElem(); } static Comb& GetComb() { return comb; } static int BestNttSize(int x) noexcept { assert(x); return 1 << MsbIndex(x*2-1); } Fps move(){ return std::move(*this); } Fps& set(int i, Elem c){ a[i] = c; return *this; } Fps& removeLeadingZeros(){ int newsz = size(); while(newsz && a[newsz-1].val() == 0) newsz--; a.resize(newsz); if((int)a.capacity() / 4 > newsz) a.shrink_to_fit(); return *this; } FpsNtt(){} FpsNtt(int sz) : a(sz, ZeroElem()) {} FpsNtt(int sz, Elem e) : a(sz, e) {} FpsNtt(std::vector&& src) : a(std::move(src)) {} FpsNtt(const std::vector& src) : a(src) {} Fps& ntt() { capSize(BestNttSize(size())); nttInst.Butterfly(a.begin(), size()); return *this; } Fps& intt() { nttInst.IButterfly(a.begin(), a.size()); return times(Elem::raw(size()).inv()); } Fps nttDouble(Fps vanilla) const { int n = size(); assert(n != 0 && n == (n&-n)); // n is a power of 2 Elem q = Elem::raw(zeta).pow((Elem::mod() - 1) / (n*2)); Elem qq = OneElem(); for(int i=0; i srcR = max(srcL, size()); // if resSz is unspecified -> resSz = destL + srcR - srcL Fps clip(int srcL, int srcR = -1, int destL = 0, int resSz = -1) const { srcR = RSZ(srcR); if(resSz < 0) resSz = destL + srcR - srcL; int rj = std::min(std::min(srcR, size()) - srcL, resSz - destL); Fps res(resSz); for(int j=std::max(0, -srcL); j b.size()) return convolution(b, a, sz); if(sz < 0) sz = std::max(0, a.size() + b.size() - 1); std::vector res(sz); for(int i=0; i=1; i--) a[i] = a[i-1] * comb.invOf(i); return set(0, ZeroElem()); } Fps log(int sz = -1){ RSZ(sz); assert(sz != 0); assert(a[0].val() == 1); return convolution(inv(sz), clip().difference(), sz-1).integral(); } Fps exp(int sz = -1){ RSZ(sz); Fps res = Fps(1).set(0, OneElem()); while(res.size() < sz){ auto z = res.size(); auto tmp = res.capSize(z*2).log().set(0, -OneElem()).move(); for(int i=0; i= (n-1) / k + 1) return Fps(n); Elem a0 = a[ctz]; return clip(ctz, ctz+n-ctz*k).times(a0.inv()).log().times(Elem(k)).exp().times(a0.pow(k)).clip(0, -1, ctz*k); } auto begin(){ return a.begin(); } auto end(){ return a.end(); } auto begin() const { return a.begin(); } auto end() const { return a.end(); } std::string toString(std::string beg = "[ ", std::string delim = " ", std::string en = " ]") const { std::string res = beg; bool f = false; for(auto x : a){ if(f){ res += delim; } f = true; res += std::to_string(x.val()); } res += en; return res; } std::vector getVectorMoved(){ return std::move(a); } Fps& operator+=(const Fps& r){ capSize(std::max(size(), r.size())); for(int i=0; i=0; i--) res = res * x + a[i]; return res; } }; template Comb FpsNtt::comb; template const NttInst FpsNtt::nttInst; } // namespace nachia using Fps = nachia::FpsNtt; namespace nachia{ namespace prime_sieve_explicit_internal{ std::vector isprime = { false }; // a[x] := isprime(2x+1) void CalcIsPrime(int z){ if((int)isprime.size() *2+1 < z+1){ int new_z = isprime.size(); while(new_z*2+1 < z+1) new_z *= 2; z = new_z-1; isprime.resize(z+1, true); for(int i=1; i*(i+1)*2<=z; i++) if(isprime[i]){ for(int j=i*(i+1)*2; j<=z; j+=i*2+1) isprime[j] = false; } } } std::vector prime_list = {2}; int prime_list_max = 0; void CalcPrimeList(int z){ while((int)prime_list.size() < z){ if((int)isprime.size() <= prime_list_max + 1) CalcIsPrime(prime_list_max * 2 + 10); for(int p=prime_list_max+1; p<(int)isprime.size(); p++){ if(isprime[p]) prime_list.push_back(p*2+1); } prime_list_max = isprime.size() - 1; } } void CalcPrimeListUntil(int z){ if(prime_list_max < z){ CalcIsPrime(z); for(int p=prime_list_max+1; p<(int)isprime.size(); p++){ if(isprime[p]) prime_list.push_back(p*2+1); } prime_list_max = isprime.size() - 1; } } } bool IsprimeExplicit(int n){ using namespace prime_sieve_explicit_internal; if(n == 2) return true; if(n % 2 == 0) return false; CalcIsPrime(n); return isprime[(n-1)/2]; } int NthPrimeExplicit(int n){ using namespace prime_sieve_explicit_internal; CalcPrimeList(n); return prime_list[n]; } int PrimeCountingExplicit(int n){ using namespace prime_sieve_explicit_internal; if(n < 2) return 0; CalcPrimeListUntil(n); auto res = std::upper_bound(prime_list.begin(), prime_list.end(), n) - prime_list.begin(); return (int)res; } // [l, r) std::vector SegmentedSieveExplicit(long long l, long long r){ assert(0 <= l); assert(l <= r); long long d = r - l; if(d == 0) return {}; std::vector res(d, true); for(long long p=2; p*p<=r; p++) if(IsprimeExplicit(p)){ long long il = (l+p-1)/p, ir = (r+p-1)/p; if(il <= p) il = p; for(long long i=il; i void DivisorZeta(std::vector& a){ int n = a.size() - 1; for(int d=2; d<=n; d++) if(IsprimeExplicit(d)) for(int i=1; i*d<=n; i++) a[i*d] += a[i]; } template void DivisorReversedZeta(std::vector& a){ int n = a.size() - 1; for(int d=2; d<=n; d++) if(IsprimeExplicit(d)) for(int i=n/d; i>=1; i--) a[i] += a[i*d]; } template void DivisorMobius(std::vector& a){ int n = a.size() - 1; for(int d=2; d<=n; d++) if(IsprimeExplicit(d)) for(int i=n/d; i>=1; i--) a[i*d] -= a[i]; } template void DivisorReversedMobius(std::vector& a){ int n = a.size() - 1; for(int d=2; d<=n; d++) if(IsprimeExplicit(d)) for(int i=1; i*d<=n; i++) a[i] -= a[i*d]; } template std::vector GcdConvolution(std::vector a, std::vector b){ assert(a.size() == b.size()); assert(1 <= a.size()); DivisorReversedZeta(a); DivisorReversedZeta(b); for(int i=1; i<(int)a.size(); i++) a[i] *= b[i]; DivisorReversedMobius(a); return a; } template std::vector LcmConvolution(std::vector a, std::vector b){ assert(a.size() == b.size()); assert(1 <= a.size()); DivisorZeta(a); DivisorZeta(b); for(int i=1; i<(int)a.size(); i++) a[i] *= b[i]; DivisorMobius(a); return a; } template void SumForCoprimeIndex(std::vector& f){ if((int)f.size() <= 1) return; Elem q = f[1]; for(int i=2; i<(int)f.size(); i++) q += f[i]; std::vector F(f.size()); F[1] = -1; DivisorMobius(F); DivisorReversedZeta(f); f[1] -= f[1]; Elem t = f[1]; for(int i=2; i<(int)f.size(); i++){ if(F[i] == 0) f[i] = f[1]; if(F[i] == -1){ t = f[1]; t -= f[i]; f[i] = t; } } DivisorZeta(f); for(int i=1; i<(int)f.size(); i++){ t = q; t -= f[i]; f[i] = t; } } } // namespace nachia using namespace std; Modint solve(int c, vector F, const vector& logPrefExp){ auto f = Fps(c+1); for(int i=1; i<=c; i++){ Modint t = F[i]; for(int j=0; j<=c; j++) f[j] += logPrefExp[i][j] * t; } return f.exp().timesFactorial().getCoeff(c); } void testcase(){ int N, K; cin >> N >> K; vector A(K+1); rep(i,N){ int a; cin >> a; A[min(a,K)] += 1; } vector I(K+1); for(int k=1; k<=K; k++) if(K%k == 0) I[k] = k; nachia::DivisorMobius(I); Modint ans = 0; vector logPrefExp(K+1, Fps(K+1)); rep(i,K+1){ rep(j,i+1) logPrefExp[i][j] = 1; logPrefExp[i] = logPrefExp[i].timesInvFactorial().log(); } for(int k=1; k<=K; k++) if(K%k == 0){ int c = K / k; vector X(c+1); for(int i=k; i<=K; i++) X[i/k] += A[i]; auto ansk = solve(c, move(X), logPrefExp); ans += ansk * I[k]; } ans /= K; cout << ans.val() << endl; } int main(){ ios::sync_with_stdio(false); cin.tie(nullptr); testcase(); return 0; }