#include #include using namespace std; using namespace atcoder; //using mint = modint1000000007; //const int mod = 1000000007; using mint = modint998244353; const int mod = 998244353; //const int INF = 1e9; //const long long LINF = 1e18; #define rep(i, n) for (int i = 0; i < (n); ++i) #define rep2(i,l,r)for(int i=(l);i<(r);++i) #define rrep(i, n) for (int i = (n) - 1; i >= 0; --i) #define rrep2(i,l,r)for(int i=(r) - 1;i>=(l);--i) #define all(x) (x).begin(),(x).end() #define allR(x) (x).rbegin(),(x).rend() #define P pair template inline bool chmax(A & a, const B & b) { if (a < b) { a = b; return true; } return false; } template inline bool chmin(A & a, const B & b) { if (a > b) { a = b; return true; } return false; } // combination mod prime // https://www.youtube.com/watch?v=8uowVvQ_-Mo&feature=youtu.be&t=1619 struct combination { std::vector fact, ifact; combination(int n) :fact(n + 1), ifact(n + 1) { assert(n < mod); fact[0] = 1; for (int i = 1; i <= n; ++i) fact[i] = fact[i - 1] * i; ifact[n] = fact[n].inv(); for (int i = n; i >= 1; --i) ifact[i - 1] = ifact[i] * i; } mint operator()(int n, int k) { return com(n, k); } mint com(int n, int k) { //負の二項係数を考慮する場合にコメントアウトを外す //if (n < 0) return com(-n, k) * (k % 2 ? -1 : 1); if (k < 0 || k > n) return 0; return fact[n] * ifact[k] * ifact[n - k]; } mint comsub(long long n, long long k) { if (n - k < k) k = n - k; assert(k < (int)fact.size()); mint val = ifact[k]; for (int i = 0; i < k; ++i) val *= n - i; return val; } template ...>, std::nullptr_t> = nullptr> mint polynom(const int n, const Ms & ...ms) { mint res = fact[n]; int sum = 0; for (int m : { ms... }) { if (m < 0 or m > n) return 0; res *= ifact[m]; sum += m; } if (sum > n)return 0; res *= ifact[n - sum]; return res; } mint div(int x) { if (x >= (int)fact.size())return mint(x).inv(); return fact[x - 1] * ifact[x]; } mint inv(int n, int k) { //if (n < 0) return inv(-n, k) * (k % 2 ? -1 : 1); if (k < 0 || k > n) return 0; return ifact[n] * fact[k] * fact[n - k]; } mint p(int n, int k) { return fact[n] * ifact[n - k]; } }com(2000006); #include #include #include // https://opt-cp.com/fps-implementation/ // verified by: // https://judge.yosupo.jp/problem/convolution_mod // https://judge.yosupo.jp/problem/inv_of_formal_power_series // https://judge.yosupo.jp/problem/log_of_formal_power_series // https://judge.yosupo.jp/problem/exp_of_formal_power_series // https://judge.yosupo.jp/problem/pow_of_formal_power_series // https://judge.yosupo.jp/problem/polynomial_taylor_shift // https://judge.yosupo.jp/problem/bernoulong longi_number // https://judge.yosupo.jp/problem/sharp_p_subset_sum //using mint = modint998244353; template struct Factorial { int MAX; std::vector fac, finv; Factorial(int m = 0) : MAX(m), fac(m + 1, 1), finv(m + 1, 1) { for (int i = 2; i < MAX + 1; ++i) fac[i] = fac[i - 1] * i; finv[MAX] /= fac[MAX]; for (int i = MAX; i >= 3; --i)finv[i - 1] = finv[i] * i; } T binom(int n, int k) { if (k < 0 || n < k) return 0; return fac[n] * finv[k] * finv[n - k]; } T perm(int n, int k) { if (k < 0 || n < k) return 0; return fac[n] * finv[n - k]; } }; Factorial fc; std::istream &operator>>(std::istream &is, modint998244353 &a) { long long v; is >> v; a = v; return is; } std::ostream &operator<<(std::ostream &os, const modint998244353 &a) { return os << a.val(); } std::istream &operator>>(std::istream &is, modint1000000007 &a) { long long v; is >> v; a = v; return is; } std::ostream &operator<<(std::ostream &os, const modint1000000007 &a) { return os << a.val(); } template std::istream &operator>>(std::istream &is, static_modint &a) { long long v; is >> v; a = v; return is; } template std::istream &operator>>(std::istream &is, dynamic_modint &a) { long long v; is >> v; a = v; return is; } template std::ostream &operator<<(std::ostream &os, const static_modint &a) { return os << a.val(); } template std::ostream &operator<<(std::ostream &os, const dynamic_modint &a) { return os << a.val(); } template std::istream &operator>>(std::istream &is, std::vector &v) { for (auto &e : v) is >> e; return is; } template std::ostream &operator<<(std::ostream &os, const std::vector &v) { for (auto &e : v) os << e << ' '; return os; } template struct FormalPowerSeries : std::vector { using std::vector::vector; using std::vector::operator=; using F = FormalPowerSeries; F operator-() const { F res(*this); for (auto &e : res) e = -e; return res; } F &operator*=(const T &g) { for (auto &e : *this) e *= g; return *this; } F &operator/=(const T &g) { assert(g != T(0)); *this *= g.inv(); return *this; } F &operator+=(const F &g) { int n = this->size(), m = g.size(); for (int i = 0; i < std::min(n, m); ++i) (*this)[i] += g[i]; return *this; } F &operator-=(const F &g) { int n = this->size(), m = g.size(); for (int i = 0; i < std::min(n, m); ++i)(*this)[i] -= g[i]; return *this; } F &operator<<=(const int d) { int n = this->size(); if (d >= n) *this = F(n); this->insert(this->begin(), d, 0); this->resize(n); return *this; } F &operator>>=(const int d) { int n = this->size(); this->erase(this->begin(), this->begin() + min(n, d)); this->resize(n); return *this; } // O(n log n) F inv(int d = -1) const { int n = this->size(); assert(n != 0 && (*this)[0] != 0); if (d == -1) d = n; assert(d >= 0); F res{ (*this)[0].inv() }; for (int m = 1; m < d; m *= 2) { F f(this->begin(), this->begin() + std::min(n, 2 * m)); F g(res); f.resize(2 * m), internal::butterfly(f); g.resize(2 * m), internal::butterfly(g); for (int i = 0; i < 2 * m; ++i) f[i] *= g[i]; internal::butterfly_inv(f); f.erase(f.begin(), f.begin() + m); f.resize(2 * m), internal::butterfly(f); for (int i = 0; i < 2 * m; ++i) f[i] *= g[i]; internal::butterfly_inv(f); T iz = T(2 * m).inv(); iz *= -iz; for (int i = 0; i < m; ++i) f[i] *= iz; res.insert(res.end(), f.begin(), f.begin() + m); } res.resize(d); return res; } // fast: FMT-friendly modulus only // O(n log n) F &multiply_inplace(const F &g, int d = -1) { int n = this->size(); if (d == -1) d = n; assert(d >= 0); *this = convolution(move(*this), g); this->resize(d); return *this; } F multiply(const F &g, const int d = -1) const { return F(*this).multiply_inplace(g, d); } // O(n log n) F ÷_inplace(const F &g, int d = -1) { int n = this->size(); if (d == -1) d = n; assert(d >= 0); *this = convolution(move(*this), g.inv(d)); this->resize(d); return *this; } F divide(const F &g, const int d = -1) const { return F(*this).divide_inplace(g, d); } // // naive // // O(n^2) // F &multiply_inplace(const F &g) { // int n = this->size(), m = g.size(); // rrep(i, n) { // (*this)[i] *= g[0]; // rep2(j, 1, min(i+1, m)) (*this)[i] += (*this)[i-j] * g[j]; // } // return *this; // } // F multiply(const F &g) const { return F(*this).multiply_inplace(g); } // // O(n^2) // F ÷_inplace(const F &g) { // assert(g[0] != T(0)); // T ig0 = g[0].inv(); // int n = this->size(), m = g.size(); // rep(i, n) { // rep2(j, 1, min(i+1, m)) (*this)[i] -= (*this)[i-j] * g[j]; // (*this)[i] *= ig0; // } // return *this; // } // F divide(const F &g) const { return F(*this).divide_inplace(g); } // sparse // O(nk) F &multiply_inplace(std::vector> g) { int n = this->size(); auto[d, c] = g.front(); if (d == 0) g.erase(g.begin()); else c = 0; for (int i = (n - 1); i >= 0; --i) { (*this)[i] *= c; for (auto &[j, b] : g) { if (j > i) break; (*this)[i] += (*this)[i - j] * b; } } return *this; } F multiply(const std::vector> &g) const { return F(*this).multiply_inplace(g); } // O(nk) F ÷_inplace(std::vector> g) { int n = this->size(); auto[d, c] = g.front(); assert(d == 0 && c != T(0)); T ic = c.inv(); g.erase(g.begin()); for (int i = 0; i < n; ++i) { for (auto &[j, b] : g) { if (j > i) break; (*this)[i] -= (*this)[i - j] * b; } (*this)[i] *= ic; } return *this; } F divide(const std::vector> &g) const { return F(*this).divide_inplace(g); } // multiply and divide (1 + cz^d) // O(n) void multiply_inplace(const int d, const T c) { int n = this->size(); if (c == T(1)) for (int i = (n - d - 1); i >= 0; --i) (*this)[i + d] += (*this)[i]; else if (c == T(-1)) for (int i = (n - d - 1); i >= 0; --i) (*this)[i + d] -= (*this)[i]; else for (int i = (n - d - 1); i >= 0; --i) (*this)[i + d] += (*this)[i] * c; } // O(n) void divide_inplace(const int d, const T c) { int n = this->size(); if (c == T(1)) for (int i = 0; i < n - d; ++i) (*this)[i + d] -= (*this)[i]; else if (c == T(-1)) for (int i = 0; i < n - d; ++i) (*this)[i + d] += (*this)[i]; else for (int i = 0; i < n - d; ++i) (*this)[i + d] -= (*this)[i] * c; } // O(n) T eval(const T &a) const { T x(1), res(0); for (auto e : *this) res += e * x, x *= a; return res; } // O(n) F &integ_inplace() { int n = this->size(); assert(n > 0); if (n == 1) return *this = F{ 0 }; this->insert(this->begin(), 0); this->pop_back(); std::vector inv(n); inv[1] = 1; int p = T::mod(); for (int i = 2; i < n; ++i) inv[i] = -inv[p%i] * (p / i); for (int i = 2; i < n; ++i) (*this)[i] *= inv[i]; return *this; } F integ() const { return F(*this).integ_inplace(); } // O(n) F &deriv_inplace() { int n = this->size(); assert(n > 0); for (int i = 2; i < n; ++i) (*this)[i] *= i; this->erase(this->begin()); this->push_back(0); return *this; } F deriv() const { return F(*this).deriv_inplace(); } // O(n log n) F &log_inplace(int d = -1) { int n = this->size(); assert(n > 0 && (*this)[0] == 1); if (d == -1) d = n; assert(d >= 0); if (d < n) this->resize(d); F f_inv = this->inv(); this->deriv_inplace(); this->multiply_inplace(f_inv); this->integ_inplace(); return *this; } F log(const int d = -1) const { return F(*this).log_inplace(d); } // O(n log n) // https://arxiv.org/abs/1301.5804 (Figure 1, right) F &exp_inplace(int d = -1) { int n = this->size(); assert(n > 0 && (*this)[0] == 0); if (d == -1) d = n; assert(d >= 0); F g{ 1 }, g_fft{ 1, 1 }; (*this)[0] = 1; this->resize(d); F h_drv(this->deriv()); for (int m = 2; m < d; m *= 2) { // prepare F f_fft(this->begin(), this->begin() + m); f_fft.resize(2 * m), internal::butterfly(f_fft); // Step 2.a' // { F _g(m); for (int i = 0; i < m; ++i) _g[i] = f_fft[i] * g_fft[i]; internal::butterfly_inv(_g); _g.erase(_g.begin(), _g.begin() + m / 2); _g.resize(m), internal::butterfly(_g); for (int i = 0; i < m; ++i) _g[i] *= g_fft[i]; internal::butterfly_inv(_g); _g.resize(m / 2); _g /= T(-m) * m; g.insert(g.end(), _g.begin(), _g.begin() + m / 2); // } // Step 2.b'--d' F t(this->begin(), this->begin() + m); t.deriv_inplace(); // { // Step 2.b' F r{ h_drv.begin(), h_drv.begin() + m - 1 }; // Step 2.c' r.resize(m); internal::butterfly(r); for (int i = 0; i < m; ++i) r[i] *= f_fft[i]; internal::butterfly_inv(r); r /= -m; // Step 2.d' t += r; t.insert(t.begin(), t.back()); t.pop_back(); // } // Step 2.e' if (2 * m < d) { t.resize(2 * m); internal::butterfly(t); g_fft = g; g_fft.resize(2 * m); internal::butterfly(g_fft); for (int i = 0; i < 2 * m; ++i) t[i] *= g_fft[i]; internal::butterfly_inv(t); t.resize(m); t /= 2 * m; } else { // この場合分けをしても数パーセントしか速くならない F g1(g.begin() + m / 2, g.end()); F s1(t.begin() + m / 2, t.end()); t.resize(m / 2); g1.resize(m), internal::butterfly(g1); t.resize(m), internal::butterfly(t); s1.resize(m), internal::butterfly(s1); for (int i = 0; i < m; ++i) s1[i] = g_fft[i] * s1[i] + g1[i] * t[i]; for (int i = 0; i < m; ++i) t[i] *= g_fft[i]; internal::butterfly_inv(t); internal::butterfly_inv(s1); for (int i = 0; i < m / 2; ++i) t[i + m / 2] += s1[i]; t /= m; } // Step 2.f' F v(this->begin() + m, this->begin() + std::min(d, 2 * m)); v.resize(m); t.insert(t.begin(), m - 1, 0); t.push_back(0); t.integ_inplace(); for (int i = 0; i < m; ++i) v[i] -= t[m + i]; // Step 2.g' v.resize(2 * m); internal::butterfly(v); for (int i = 0; i < 2 * m; ++i) v[i] *= f_fft[i]; internal::butterfly_inv(v); v.resize(m); v /= 2 * m; // Step 2.h' for (int i = 0; i < std::min(d - m, m); ++i)(*this)[m + i] = v[i]; } return *this; } F exp(const int d = -1) const { return F(*this).exp_inplace(d); } // O(n log n) F &pow_inplace(const long long k, int d = -1) { int n = this->size(); if (d == -1) d = n; assert(d >= 0 && k >= 0); if (k == 0) { *this = F(d); if (d > 0) (*this)[0] = 1; return *this; } int l = 0; while (l < n && (*this)[l] == 0) ++l; if (l > (d - 1) / k || l == n) return *this = F(d); T c = (*this)[l]; this->erase(this->begin(), this->begin() + l); *this /= c; this->log_inplace(d - l * k); *this *= k; this->exp_inplace(); *this *= c.pow(k); this->insert(this->begin(), l*k, 0); return *this; } F pow(const long long k, const int d = -1) const { return F(*this).pow_inplace(k, d); } // O(n log n) F &shift_inplace(const T c) { int n = this->size(); fc = Factorial(n); for (int i = 0; i < n; ++i) (*this)[i] *= fc.fac[i]; reverse(this->begin(), this->end()); F g(n); T cp = 1; for (int i = 0; i < n; ++i) g[i] = cp * fc.finv[i], cp *= c; this->multiply_inplace(g, n); reverse(this->begin(), this->end()); for (int i = 0; i < n; ++i) (*this)[i] *= fc.finv[i]; return *this; } F shift(const T c) const { return F(*this).shift_inplace(c); } F operator*(const T &g) const { return F(*this) *= g; } F operator/(const T &g) const { return F(*this) /= g; } F operator+(const F &g) const { return F(*this) += g; } F operator-(const F &g) const { return F(*this) -= g; } F operator<<(const int d) const { return F(*this) <<= d; } F operator>>(const int d) const { return F(*this) >>= d; } F operator*(std::vector> g) const { return F(*this) *= g; } F operator/(std::vector> g) const { return F(*this) /= g; } }; using fps = FormalPowerSeries; using sfps = std::vector>; // https://judge.yosupo.jp/problem/inv_of_formal_power_series // f(x)に対して 1/f(x)の先頭n項を求める void yosupo_inv() { int n; std::cin >> n; fps a(n); std::cin >> a; std::cout << a.inv() << '\n'; } // https://judge.yosupo.jp/problem/log_of_formal_power_series // f(x)に対して log(f(x))の先頭n項を求める void yosupo_log() { int n; std::cin >> n; fps a(n); std::cin >> a; std::cout << a.log_inplace() << '\n'; } // https://judge.yosupo.jp/problem/exp_of_formal_power_series // f(x)に対して exp(f(x))の先頭n項を求める void yosupo_exp() { int n; std::cin >> n; fps a(n); std::cin >> a; std::cout << a.exp_inplace() << '\n'; } // https://judge.yosupo.jp/problem/pow_of_formal_power_series // f(x)に対してf(x)^mの先頭n項を求める void yosupo_pow() { int n, m; std::cin >> n >> m; fps a(n); std::cin >> a; std::cout << a.pow_inplace(m) << '\n'; } // https://judge.yosupo.jp/problem/polynomial_taylor_shift // f(x)に対してf(x+c)の先頭n項を求める void yosupo_shift() { int n, c; std::cin >> n >> c; fps a(n); std::cin >> a; std::cout << a.shift_inplace(c) << '\n'; } // https://judge.yosupo.jp/problem/bernoulli_number // ベルヌーイ数の先頭n項を求める void yosupo_bernoulli() { int n; std::cin >> n; ++n; fc = Factorial(n); fps f((fc.finv).begin() + 1, (fc.finv).end()); f = f.inv(); for (int i = 0; i < n; ++i) f[i] *= fc.fac[i]; std::cout << f << '\n'; } // https://judge.yosupo.jp/problem/sharp_p_subset_sum // Π(1+x^a)の先頭t+1項を求める void yosupo_count_subset_sum() { int n, t; std::cin >> n >> t; ++t; std::vector c(t); for (int i = 0; i < n; ++i) { int s; std::cin >> s; ++c[s]; } std::vector inv(t); inv[1] = 1; int p = mint::mod(); for (int i = 2; i < t; ++i) inv[i] = -inv[p%i] * (p / i); fps f(t); for (int i = 0; i < t; ++i) { if (c[i] == 0) continue; for (int j = 1, d = i; d < t; ++j, d += i) { if (j & 1) f[d] += c[i] * inv[j]; else f[d] -= c[i] * inv[j]; } } f.exp_inplace(); f.erase(f.begin()); std::cout << f << '\n'; } int main() { ios::sync_with_stdio(false); cin.tie(nullptr); int n, k; cin >> n >> k; vectorcnt(k + 1); rep(i, n) { int a; cin >> a; chmin(a, k); cnt[a]++; } // gcd周期になる vectorg_cnt(k + 1); rep(i, k) { int g = gcd(i, k); g_cnt[g]++; } vectorf(k + 1, fps(k + 1)); rep(i, k + 1) { if (i != 0)f[i] = f[i - 1]; f[i][i] = com.ifact[i]; } vectorlog_f(k + 1, fps(k + 1)); rep(i, k + 1) { log_f[i] = f[i].log(); } mint ans = 0; for (int p = 0; p < k + 1; ++p) { int c = g_cnt[p]; if (c == 0)continue; int q = k / p; vectors_cnt(p + 1); rep(i, k + 1) s_cnt[i / q] += cnt[i]; mint tmp = 0; { // exp^log(f) = f fps log(p + 1); rep(i, p + 1) log += log_f[i] * s_cnt[i]; log.exp_inplace(); tmp += log[p]; } ans += tmp * com.fact[p] * c; } ans /= k; cout << ans.val() << endl; return 0; }