#pragma GCC target("avx2") #pragma GCC optimize("O3") #pragma GCC optimize("unroll-loops") #pragma GCC optimize("Ofast") #include #include using namespace std; #if __has_include() #include using namespace atcoder; #endif using ll = long long; using ld = long double; using ull = long long; #define endl "\n" typedef pair Pii; #define REP3(i, m, n) for (int i = (m); (i) < int(n); ++ (i)) #define rep(i,a,b) for(int i=(int)(a);i<(int)(b);i++) #define ALL(x) begin(x), end(x) #define all(s) (s).begin(),(s).end() //#define rep2(i, m, n) for (int i = (m); i < (n); ++i) //#define rep(i, n) rep2(i, 0, n) #define PB push_back #define drep2(i, m, n) for (int i = (m)-1; i >= (n); --i) #define drep(i, n) drep2(i, n, 0) #define rever(vec) reverse(vec.begin(), vec.end()) #define sor(vec) sort(vec.begin(), vec.end()) //#define FOR(i,a,b) for(ll i=a;i<=(ll)(b);i++) #define fi first #define se second #define pb push_back #define P pair #define NP next_permutation //const ll mod = 998244353; const ll mod = 1000000007; const ll inf = 4500000000000000000ll; const ld eps = ld(0.000000001); static const long double pi = 3.141592653589793; templatevoid vcin(vector &n){for(int i=0;i>n[i];} templatevoid vcin(vector &n,vector &m){for(int i=0;i>n[i]>>m[i];} templatevoid vcout(vector &n){for(int i=0;ivoid vcin(vector> &n){for(int i=0;i>n[i][j];}}} templatevoid vcout(vector> &n){for(int i=0;iauto min(const T& a){ return *min_element(all(a)); } templateauto max(const T& a){ return *max_element(all(a)); } templatevoid print(pair a){cout<bool chmax(T &a, const T &b) { if (abool chmin(T &a, const T &b) { if (b void ifmin(T t,T u){if(t>u){cout<<-1< void ifmax(T t,T u){if(t>u){cout<<-1<auto make_vector(T x,int arg,Args ...args){if constexpr(sizeof...(args)==0)return vector(arg,x);else return vector(arg,make_vector(x,args...));} ll modPow(ll a, ll n, ll mod) { if(mod==1) return 0;ll ret = 1; ll p = a % mod; while (n) { if (n & 1) ret = ret * p % mod; p = p * p % mod; n >>= 1; } return ret; } vector divisor(ll x){ vector ans; for(ll i = 1; i * i <= x; i++){ if(x % i == 0) {ans.push_back(i); if(i*i!=x){ ans.push_back(x / ans[i]);}}}sor(ans); return ans; } ll pop(ll a){ll res=0;while(a){res+=a%2;a/=2;}return res;} template struct Sum{ vector data; Sum(const vector& v):data(v.size()+1){ for(ll i=0;i struct Sum2{ vector> data; Sum2(const vector> &v):data(v.size()+1,vector(v[0].size()+1)){ for(int i=0;i vector NTT(vector a,vector b){ ll nmod=T::mod(); int n=a.size(); int m=b.size(); vector x1(n); vector y1(m); for(int i=0;i(x1,y1); auto z2=convolution<469762049>(x1,y1); auto z3=convolution<1224736769>(x1,y1); vector res(n+m-1); ll m1=167772161; ll m2=469762049; ll m3=1224736769; ll m1m2=104391568; ll m1m2m3=721017874; ll mm12=m1*m2%nmod; for(int i=0;i struct FormalPowerSeries : std::vector { using std::vector::vector; using std::vector::size; using std::vector::resize; using F = FormalPowerSeries; F &operator+=(const F &g){ for(int i=0;i>=(const int d) { int n=(*this).size(); (*this).erase((*this).begin(),(*this).begin()+min(n, d)); (*this).resize(n); return *this; } F &operator=(const std::vector &v) { int n = (*this).size(); for(int i = 0; i < n; ++i) (*this)[i] = v[i]; return *this; } F operator-() const { F ret = *this; return ret * -1; } F &operator*=(const F &g) { if(mode==FAST) { int n=(*this).size(); auto tmp=atcoder::convolution(*this,g); int k=tmp.size(); (*this).resize(k); for(int i=0;i>(const int d) const { return F(*this)>>=d;} void onemul(const int d,const T c){ int n=(*this).size(); for(int i=n-d-1;i>=0;i--){ (*this)[i+d]+=(*this)[i]*c; } } void onediv(const int d,const T c){ int n=(*this).size(); for(int i=0;i 0); F res{(*this)[0].inv()}; while(int(res.size()) < deg) { int m = res.size(); F f((*this).begin(), (*this).begin() + std::min(n, m * 2)), r(res); f.resize(m * 2), atcoder::internal::butterfly(f); r.resize(m * 2), atcoder::internal::butterfly(r); for(int i = 0; i < m * 2; ++i) f[i] *= r[i]; atcoder::internal::butterfly_inv(f); f.erase(f.begin(), f.begin() + m); f.resize(m * 2), atcoder::internal::butterfly(f); for(int i = 0; i < m * 2; ++i) f[i] *= r[i]; atcoder::internal::butterfly_inv(f); T iz = T(m * 2).inv(); iz *= -iz; for(int i = 0; i < m; ++i) f[i] *= iz; res.insert(res.end(), f.begin(), f.begin() + m); } res.resize(deg); return res; } else{ assert(n!=0&&(*this)[0]!=0); if(deg==-1) deg=n; assert(deg>0); F res{(*this)[0].inv()}; while(res.size() s=NTT(f,r); s.resize(2*m); for(int i=0;i<2*m;i++){ s[i]=-s[i]; } s[0]+=2; vector g=NTT(s,r); g.resize(2*m); swap(res,g); } res.resize(n); return res; } } F &diff_inplace() { int n = (*this).size(); for(int i = 1; i < n; ++i) (*this)[i - 1] = (*this)[i] * i; (*this)[n - 1] = 0; return *this; } F diff() const { F(*this).diff_inplace();} F &integral_inplace() { int n = (*this).size(), mod = T::mod(); std::vector inv(n); { inv[1] = 1; for(int i = 2; i < n; ++i) inv[i] = T(mod) - inv[mod % i] * (mod / i); } for(int i = n - 2; i >= 0; --i) (*this)[i + 1] = (*this)[i] * inv[i + 1]; (*this)[0] = 0; return *this; } F integral() const { return F(*this).integral_inplace(); } F &log_inplace() { int n = (*this).size(); assert(n and (*this)[0] == 1); F f_inv = (*this).inv(); (*this).diff_inplace(); (*this) *= f_inv; (*this).integral_inplace(); return *this; } F log() const { return F(*this).log_inplace(); } F &deriv_inplace() { int n = (*this).size(); assert(n); 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(); } F &exp_inplace() { int n = (*this).size(); assert(n and (*this)[0] == 0); F g{1}; (*this)[0] = 1; F h_drv((*this).deriv()); for(int m = 1; m < n; m *= 2) { F f((*this).begin(), (*this).begin() + m); f.resize(2 * m), atcoder::internal::butterfly(f); auto mult_f = [&](F &p) { p.resize(2 * m); atcoder::internal::butterfly(p); for(int i = 0; i < 2 * m; ++i) p[i] *= f[i]; atcoder::internal::butterfly_inv(p); p /= 2 * m; }; if(m > 1) { F g_(g); g_.resize(2 * m), atcoder::internal::butterfly(g_); for(int i = 0; i < 2 * m; ++i) g_[i] *= g_[i] * f[i]; atcoder::internal::butterfly_inv(g_); T iz = T(-2 * m).inv(); g_ *= iz; g.insert(g.end(), g_.begin() + m / 2, g_.begin() + m); } F t((*this).begin(), (*this).begin() + m); t.deriv_inplace(); { F r{h_drv.begin(), h_drv.begin() + m - 1}; mult_f(r); for(int i = 0; i < m; ++i) t[i] -= r[i] + r[m + i]; } t.insert(t.begin(), t.back()); t.pop_back(); t *= g; F v((*this).begin() + m, (*this).begin() + std::min(n, 2 * m)); v.resize(m); t.insert(t.begin(), m - 1, 0); t.push_back(0); t.integral_inplace(); for(int i = 0; i < m; ++i) v[i] -= t[m + i]; mult_f(v); for(int i = 0; i < std::min(n - m, m); ++i) (*this)[m + i] = v[i]; } return *this; } F exp() const { return F(*this).exp_inplace(); } F &pow_inplace(long long k) { int n = (*this).size(), l = 0; assert(k >= 0); if(!k){ for(int i = 0; i < n; ++i) (*this)[i] = !i; return *this; } while(l < n and (*this)[l] == 0) ++l; if(l > (n - 1) / k or l == n) return *this = F(n); T c = (*this)[l]; (*this).erase((*this).begin(), (*this).begin() + l); (*this) /= c; (*this).log_inplace(); (*this).resize(n - 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 { return F(*this).pow_inplace(); } void manymul(vector> g) { int n = (*this).size(); auto [d, c] = g.front(); if (d == 0) g.erase(g.begin()); else c = 0; drep(i, n) { (*this)[i] *= c; for (auto &[j, b] : g) { if (j > i) break; (*this)[i] += (*this)[i-j] * b; } } } void manydiv(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()); rep(i, 0,n) { for (auto &[j, b] : g) { if (j > i) break; (*this)[i] -= (*this)[i-j] * b; } (*this)[i] *= ic; } } }; template void GaussJordan(vector> &A,bool is_extended = false){ ll m=A.size(),n=A[0].size(); ll rank=0; for(int i=0;i void linear_equation(vector> a, vector b, vector &res) { ll m=a.size(),n=a[0].size(); vector> M(m,vector(n+1)); for(int i=0;i pair Characteristic_equation(const F &a) { using T=typename F::value_type; ll n=a.size(); ll p=n/2; ll u=p+(p+1); vector> f(u,vector(u)); f[0][0]=1; for(int i=1;i<=p;i++){ f[i][i-1]=-1; } for(int i=p;i b(u); b[0]=1; vector res(u); linear_equation(f,b,res); F X(p),Y(p+1); for(int i=0;i T getK(FormalPowerSeries p, FormalPowerSeries q,ll k){ if(k<0) return T(0); ll d=q.size(); while(k){ auto qn=q; for(int i=1;i;*/ using mint = modint998244353; using fps = FormalPowerSeries; int main() { cincout(); ll n; cin>>n; ll k; cin>>k; k=n*(n-1)/2-k; vector a(n); vcin(a); map m; for(auto e:a) m[e]++; vector u; for(auto e:m) u.pb(e.se); fps f(n*n); ll sum=0; f[0]=1; for(int i=1;i