結果
| 問題 |
No.2243 Coaching Schedule
|
| コンテスト | |
| ユーザー |
 GOTKAKO
|
| 提出日時 | 2025-05-24 06:20:09 |
| 言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 2,956 ms / 4,000 ms |
| コード長 | 40,276 bytes |
| コンパイル時間 | 5,398 ms |
| コンパイル使用メモリ | 257,848 KB |
| 実行使用メモリ | 72,512 KB |
| 最終ジャッジ日時 | 2025-05-24 06:21:27 |
| 合計ジャッジ時間 | 69,942 ms |
|
ジャッジサーバーID （参考情報） |
judge2 / judge3 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 37 |
ソースコード
#include <bits/stdc++.h>
using namespace std;
long long mod = 998244353;
//入力が必ず-mod<a<modの時.
struct mint{
long long v; //vにはアクセス出来る.
mint():v(0){} mint(int a){v = a<0?a+mod:a;} mint(long long a){v = a<0?a+mod:a;}
mint(unsigned int a):v(a){} mint(unsigned long long a):v(a){}
long long val(){return v;} //空気.
mint &operator=(const mint &b) = default;
mint operator-()const{return mint(0)-(*this);}
mint operator+(const mint b)const{return mint(v)+=b;}
mint operator-(const mint b)const{return mint(v)-=b;}
mint operator*(const mint b)const{return mint(v)*=b;}
mint operator/(const mint b)const{return mint(v)/=b;}
mint operator+=(const mint b){
v += b.v; if(v >= mod) v -= mod;
return *this;
}
mint operator-=(const mint b){
v -= b.v; if(v < 0) v += mod;
return *this;
}
mint operator*=(const mint b){v = v*b.v%mod; return *this;}
mint operator/=(mint b){ //素数modオンリー.
if(b == 0) assert(false);
long long left = mod-2;
while(left){if(left&1) *this *= b; b *= b; left >>= 1;}
return *this;
}
mint operator++(){*this += 1; return *this;}
mint operator--(){*this -= 1; return *this;}
mint operator++(int){*this += 1; return *this;}
mint operator--(int){*this -= 1; return *this;}
bool operator==(const mint b)const{return v == b.v;}
bool operator!=(const mint b)const{return v != b.v;}
bool operator>(const mint b)const{return v > b.v;}
bool operator>=(const mint b)const{return v >= b.v;}
bool operator<(const mint b)const{return v < b.v;}
bool operator<=(const mint b)const{return v <= b.v;}
mint pow(long long n)const{ //nが超デカい時はmod P-1で取るべき場合がある.
mint ret = 1,p = v;
if(n < 0) p = p.inv(),n = -n; //負も対応.
while(n){
if(n&1) ret *= p;
p *= p; n >>= 1;
}
return ret;
}
mint inv()const{return mint(1)/v;}
};
namespace to_fold{ //纏めて折り畳むためのやつ 苦労したが結局ほぼACLの写経.
__int128_t safemod(__int128_t a,long long m){a %= m; if(a < 0) a += m; return a;}
pair<long long,long long> invgcd(long long a,long long b){
//return {gcd(a,b),x} (xa≡g(mod b))
a = safemod(a,b);
if(a == 0) return {b,0};
long long x = 0,y = 1,memob = b;
while(a){
long long q = b/a;
b -= a*q;
swap(x,y); y -= q*x;
swap(a,b);
}
if(x < 0) x += memob/b;
return {b,x};
}
long long Garner(vector<long long> &A,vector<long long> &M){
__int128_t mulM = 1,x = A.at(0)%M.at(0); //Mの要素のペア互いに素必須.
for(int i=1; i<A.size(); i++){
//assert(gcd(mulM,M.at(i-1)) == 1);
mulM *= M.at(i-1); //2乗がオーバーフローする時__int128_t
long long t = safemod((A.at(i)-x)*invgcd(mulM,M.at(i)).second,M.at(i));
x += t*mulM;
}
return x%mod;
}
int countzero(unsigned long long x){
if(x == 0) return 64;
else return __popcount((x&-x)-1);
}
struct fftinfo{
private:
long long checkmod;
public:
long long getmod(){return checkmod;}
array<mint,30> root,iroot,rate2,irate2,rate3,irate3;
fftinfo(mint g){
checkmod = mod;
int rank2 = countzero(mod-1);
root[rank2] = g.pow((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 mul = 1,imul = 1;
for(int i=0; i<=rank2-2; i++){
rate2[i] = root[i+2]*mul;
irate2[i] = iroot[i+2]*imul;
mul *= iroot[i+2],imul *= root[i+2];
}
mul = 1,imul = 1;
for(int i=0; i<=rank2-3; i++){
rate3[i] = root[i+3]*mul;
irate3[i] = iroot[i+3]*imul;
mul *= iroot[i+3],imul *= root[i+3];
}
}
};
mint findroot(){
if(mod == 998244353) return mint(3);
else if(mod == 754974721) return mint(11);
else if(mod == 167772161) return mint(3);
else if(mod == 469762049) return mint(3);
assert(false); //バグあるので禁止.
long long p = mod-1;
vector<long long> P;
for(long long i=2; i*i<=p; i++){
if(p%i) continue;
while(p%i == 0) p /= i;
P.push_back(i);
}
if(p != 1) P.push_back(p);
p = mod-1;
random_device rnd; mt19937 mt(rnd());
while(true){
mint ret = mt()%mod;
if(ret == 1 || ret == 0) continue;
if(ret.pow(p) != 1) continue;
bool ok = true;
for(auto check : P) if(ret.pow(p/check) == 1){ok = false; break;}
if(ok) return ret;
} //O(√P)らしいです.
}
fftinfo info(3);
void NTT(vector<mint> &A){
if(info.getmod() != mod) info = fftinfo(findroot());
int N = A.size(),ln = countzero(N);
int dep = 0;
while(dep < ln){
if(ln-dep == 1 || true){
int p = 1<<(ln-dep-1);
mint rot = 1;
for(int o=0; o<(1<<dep); o++){
int offset = o<<(ln-dep);
for(int i=0; i<p; i++){
mint l = A.at(i+offset),r = A.at(i+offset+p)*rot;
A.at(i+offset) = l+r; A.at(i+offset+p) = l-r;
}
if(o+1 != (1<<dep)) rot *= info.rate2[countzero(~(unsigned int)o)];
}
dep++;
}
else{
int p = 1<<(ln-dep-2);
mint rot = 1,imag = info.root[2];
for(int o=0; o<(1<<dep); o++){
mint rot2 = rot*rot,rot3 = rot2*rot;
int offset = o<<(ln-dep);
for(int i=0; i<p; i++){
auto mod2 = 1ULL*mod*mod;
auto a0 = 1ULL*A.at(i+offset).v;
auto a1 = 1ULL*A.at(i+offset+p).v*rot.v;
auto a2 = 1ULL*A.at(i+offset+p+p).v*rot2.v;
auto a3 = 1ULL*A.at(i+offset+p+p+p).v*rot3.v;
auto a1m2a3im = 1ULL*(a1+mod2-a3+mod)%mod*imag.v;
auto m2a2 = mod2-a2;
A.at(i+offset) = (long long)((a0+a2+a1+a3)%mod);
A.at(i+offset+p) = (long long)((a0+a2+(2*mod2-(a1+a3)))%mod);
A.at(i+offset+p+p) = (long long)((a0+m2a2+a1m2a3im)%mod);
A.at(i+offset+p+p+p) = (long long)((a0+m2a2+(mod2-a1m2a3im))%mod);
}
if(o+1 != (1<<dep)) rot *= info.rate3[countzero(~(unsigned int)o)];
}
dep += 2;
}
}
}
void INTT(vector<mint> &A){
if(info.getmod() != mod) info = fftinfo(findroot());
int N = A.size(),ln = countzero(N),dep = ln;
while(dep){
if(dep == 1){
int p = 1<<(ln-dep);
mint irot = 1;
for(int o=0; o<(1<<(dep-1)); o++){
int offset = o<<(ln-dep+1);
for(int i=0; i<p; i++){
mint l = A.at(i+offset),r = A.at(i+offset+p);
A.at(i+offset) = l+r;
A.at(i+offset+p) = (long long)((unsigned long long)((mod+l.v-r.v)*irot.v)%mod);
}
if(o+1 != (1<<(dep-1))) irot *= info.irate2[countzero(~(unsigned int)o)];
}
dep--;
}
else{
int p = 1<<(ln-dep);
mint irot = 1,iimag = info.iroot[2];
for(int o=0; o<(1<<(dep-2)); o++){
mint irot2 = irot*irot,irot3 = irot2*irot;
int offset = o<<(ln-dep+2);
for(int i=0; i<p; i++){
auto a0 = 1ULL*A.at(i+offset).v,a1 = 1ULL*A.at(i+offset+p).v;
auto a2 = 1ULL*A.at(i+offset+p+p).v,a3 = 1ULL*A.at(i+offset+p+p+p).v;
auto gotyagotya = 1ULL*((mod+a2-a3)*iimag.v%mod);
A.at(i+offset) = (long long)((a0+a1+a2+a3)%mod);
A.at(i+offset+p) = (long long)((a0+(mod-a1)+gotyagotya)*irot.v%mod);
A.at(i+offset+p+p) = (long long)((a0+a1+(mod-a2)+(mod-a3))*irot2.v%mod);
A.at(i+offset+p+p+p) = (long long)((a0+(mod-a1)+(mod-gotyagotya))*irot3.v%mod);
}
if(o+1 != (1<<(dep-2))) irot *= info.irate3[countzero(~(unsigned int)o)];
}
dep -= 2;
}
}
}
vector<mint> convolution(const vector<mint> &A,const vector<mint> &B){
int siza = A.size(),sizb = B.size(),sizc = siza+sizb-1,N = 1;
if(siza == 0 || sizb == 0) return {};
if(min(siza,sizb) <= 60){ //naive.
vector<mint> ret(sizc);
if(siza >= sizb){for(int i=0; i<siza; i++) for(int k=0; k<sizb; k++) ret.at(i+k) += A.at(i)*B.at(k);}
else{for(int i=0; i<sizb; i++) for(int k=0; k<siza; k++) ret.at(i+k) += B.at(i)*A.at(k);}
return ret;
}
while(N < sizc) N <<= 1;
vector<mint> ca = A,cb = B;
ca.resize(N),cb.resize(N);
assert((mod-1)%N == 0);
mint root = findroot();
NTT(ca); NTT(cb);
for(int i=0; i<N; i++) ca.at(i) *= cb.at(i);
INTT(ca); ca.resize(sizc);
mint divN = mint(N).inv();
for(int i=0; i<sizc; i++) ca.at(i) *= divN;
return ca;
}
vector<long long> convolution_ll(const vector<long long> &A,const vector<long long> &B){ //long longに収まる範囲.
int siza = A.size(),sizb = B.size(),sizc = siza+sizb-1;
if(siza == 0 || sizb == 0) return {};
vector<long long> ret(sizc);
if(min(siza,sizb) <= 200){ //naive 200はやばい?.
vector<long long> ret(sizc);
if(siza >= sizb){for(int i=0; i<siza; i++) for(int k=0; k<sizb; k++) ret.at(i+k) += A.at(i)*B.at(k);}
else{for(int i=0; i<sizb; i++) for(int k=0; k<siza; k++) ret.at(i+k) += B.at(i)*A.at(k);}
return ret;
}
unsigned long long mod1 = 754974721,mod2 = 167772161,mod3 = 469762049;
unsigned long long m1m2 = mod1*mod2,m2m3 = mod2*mod3,m3m1 = mod3*mod1,m1m2m3 = mod1*mod2*mod3;
unsigned long long i1 = invgcd(m2m3,mod1).second,i2 = invgcd(m3m1,mod2).second,i3 = invgcd(m1m2,mod3).second;
assert(sizc <= (1<<24));
mod = mod1;
vector<mint> a(siza),b(sizb);
for(int i=0; i<siza; i++) a.at(i) = A.at(i)%mod;
for(int i=0; i<sizb; i++) b.at(i) = B.at(i)%mod;
vector<mint> C1 = convolution(a,b);
mod = mod2;
for(int i=0; i<siza; i++) a.at(i) = A.at(i)%mod;
for(int i=0; i<sizb; i++) b.at(i) = B.at(i)%mod;
vector<mint> C2 = convolution(a,b);
mod = mod3;
for(int i=0; i<siza; i++) a.at(i) = A.at(i)%mod;
for(int i=0; i<sizb; i++) b.at(i) = B.at(i)%mod;
vector<mint> C3 = convolution(a,b);
vector<unsigned long long> offset = {0,0,m1m2m3,2*m1m2m3,3*m1m2m3};
for(int i=0; i<sizc; i++){
unsigned long long x = 0;
x += (C1.at(i).v*i1)%mod1*m2m3;
x += (C2.at(i).v*i2)%mod2*m3m1;
x += (C3.at(i).v*i3)%mod3*m1m2;
long long diff = C1.at(i).v-((long long)x+(long long)mod1)%mod1;
if(diff < 0) diff += mod1;
x -= offset.at(diff%5);
ret.at(i) = x;
}
return ret;
}
vector<mint> convolution_llmod(vector<mint> &A,vector<mint> &B,long long memomod){
int siza = A.size(),sizb = B.size(),sizc = siza+sizb-1;
if(siza == 0 || sizb == 0) return {};
vector<mint> ret(sizc);
if(min(siza,sizb) <= 200){
for(int i=0; i<siza; i++) for(int k=0; k<sizb; k++) ret.at(i+k) += A.at(i)*B.at(k);
return ret;
}
long long mod1 = 754974721,mod2 = 167772161,mod3 = 469762049;
assert(sizc <= (1<<24));
mod = mod1;
vector<mint> a(siza),b(sizb);
for(int i=0; i<siza; i++) a.at(i) = A.at(i).v%mod;
for(int i=0; i<sizb; i++) b.at(i) = B.at(i).v%mod;
vector<mint> C1 = convolution(a,b);
mod = mod2;
for(int i=0; i<siza; i++) a.at(i) = A.at(i).v%mod;
for(int i=0; i<sizb; i++) b.at(i) = B.at(i).v%mod;
vector<mint> C2 = convolution(a,b);
mod = mod3;
for(int i=0; i<siza; i++) a.at(i) = A.at(i).v%mod;
for(int i=0; i<sizb; i++) b.at(i) = B.at(i).v%mod;
vector<mint> C3 = convolution(a,b);
mod = memomod;
for(int i=0; i<sizc; i++){
vector<long long> A = {C1.at(i).v,C2.at(i).v,C3.at(i).v};
vector<long long> M = {mod1,mod2,mod3};
ret.at(i) = Garner(A,M);
}
return ret;
}
vector<int> convolution_int(vector<int> &A,vector<int> &B){ //intに収まる範囲.
if(A.size() == 0 || B.size() == 0) return {};
vector<int> ret;
if(min(A.size(),B.size()) <= 60){
ret.resize(A.size()+B.size()-1);
for(int i=0; i<A.size(); i++) for(int k=0; k<B.size(); k++) ret.at(i+k) += A.at(i)*B.at(k);
}
else{
vector<mint> X,Y,Z;
for(auto &a : A) X.push_back(a);
for(auto &b : B) Y.push_back(b);
Z = convolution(X,Y);
for(auto &z : Z) ret.push_back(z.v);
}
return ret;
}
}
using namespace to_fold;
template<typename T> //実質mintだけ?.
struct FormalPowerSeries:vector<T>{ //NTT-friendly素数だけ じゃなくてもいいけど全部書き直せ!.
using vector<T>::vector;
using fps = FormalPowerSeries;
//重要なところは某のほぼパクリ.
fps operator++(){*this += 1; return *this;}
fps operator--(){*this -= 1; return *this;}
fps operator++(int){*this += 1; return *this;}
fps operator--(int){*this -= 1; return *this;}
fps operator+(const fps &b) const {return fps(*this)+=b;}
fps operator+(const T &b) const {return fps(*this)+=b;}
fps operator-(const fps &b) const {return fps(*this)-=b;}
fps operator-(const T &b) const {return fps(*this)-=b;}
fps operator*(const fps &b){return fps(*this)*=b;}
fps operator*(const T &b) const {return fps(*this)*=b;}
fps operator/(const fps &b) const {return fps(*this)/=b;}
fps operator%(const fps &b) const {return fps(*this)%=b;}
fps operator>>(const unsigned int b) const {return fps(*this)>>=b;}
fps operator<<(const unsigned int b) const {return fps(*this)<<=b;}
fps operator-() const {
fps ret = (*this);
for(auto &v : ret) v = -v;
return ret;
}
bool operator==(const fps &b)const{
if((*this).size() != b.size()) return false;
for(int i=0; i<(*this).size(); i++) if((*this).at(i) != b.at(i)) return false;
return true;
}
bool operator!=(const fps &b)const{return !((*this)==b);}
fps &operator+=(const fps &b){ //Cix^i = (Ai+Bi)x^i. O(n).
if((*this).size() < b.size()) (*this).resize(b.size(),0);
for(int i=0; i<b.size(); i++) (*this).at(i) += b.at(i);
return *this;
}
fps &operator+=(const T &b){ //x^0の係数に+b O(1).
if((*this).size() == 0) (*this).resize(1);
(*this).at(0) += b;
return *this;
}
fps &operator-=(const fps &b){ //Cix^i = (Ai-Bi)x^i. O(n).
int n = b.size();
if((*this).size() < n) (*this).resize(n,0);
for(int i=0; i<n; i++) (*this).at(i) -= b.at(i);
return *this;
}
fps &operator-=(const T &b){ //x^0の係数に-b O(1).
if((*this).size() == 0) (*this).resize(1);
(*this).at(0) -= b;
return *this;
}
fps &operator*=(const fps &b){ //C[i+k]x^(i+k) = Aix^i+Bkx^k. O(nlogn).
vector<T> f,g;
for(auto v : (*this)) f.emplace_back(v);
for(auto v : b) g.emplace_back(v);
vector<T> ret;
if(mod == 998244353) ret = convolution(f,g); //型を仕方なく合わせる C++わからん.
else ret = convolution_llmod(f,g,mod); //modがNTT-friendlyじゃない時用 直すのだるい 998以外はとりあえずこっち.
(*this).resize(ret.size());
for(int i=0; i<ret.size(); i++) (*this).at(i) = ret.at(i);
return (*this);
}
fps &operator*=(const T &b){ //x^iの係数に*b O(n).
for(auto &v : (*this)) v *= b;
return *this;
}
fps &operator/=(const T &b){ //x^iの係数に/b O(logmod+n).
T mul = b.inv();
for(auto &v : (*this)) v *= mul;
return *this;
}
fps &operator>>=(const unsigned int &b){//b<0は対象外. 先頭b項を削除. O(n)
if((*this).size() <= b) (*this).clear();
else (*this).erase((*this).begin(),(*this).begin()+b);
return *this;
}
fps &operator<<=(const unsigned int &b){//b<0は対象外. 先頭b項に0を挿入. O(n)
(*this).insert((*this).begin(),b,0);
return *this;
}
fps &operator%=(const fps &b){ //多項式の余り. O(nlogn)
(*this) -= (*this)/b*b;
del0();
return (*this);
}
fps &operator/=(const fps &b){ //多項式としての除算 O(nlogn).
assert(b.size() > 0); //分母の末尾0は駄目.
T check = b.back(); assert(check != 0);
del0(); //分子の末尾0は消して許容.
if((*this).size() < b.size()){
(*this).clear();
return *this;
}
int n = (*this).size()-b.size()+1;
if(b.size() <= 64){ //愚直.
fps G(b);
assert(G.size() > 0);
T div = G.back().inv();
for (auto &v : G) v *= div;
int deg = (*this).size()-G.size()+1;
fps Q(deg);
for(int i=deg-1; i>=0; i--) {
Q[i] = (*this).at(i+G.size()-1);
for(int k=0; k<G.size(); k++) (*this).at(i+k) -= Q.at(i)*G.at(k);
}
(*this) = Q*div; (*this).resize(n,0);
return *this;
}
//rev(f)/rev(g)≡rev(ret)(mod x^n)らしい わお.
//998244353以外は定数倍悪い方に突っ込む.
if(mod == 998244353) return (*this) = ((*this).rev().prefix(n) * b.rev().inv(n)).prefix(n).rev();
return (*this) = ((*this).rev().prefix(n) * b.rev().inv2(n)).prefix(n).rev();
}
fps multi_sparse(const fps &b)const{ //非0が少ない時の掛け算 愚直 O(n*|非0|+m).
int n = (*this).size(),m = b.size();
vector<pair<int,T>> P;
for(int i=0; i<m; i++) if(b.at(i) != 0) P.push_back({i,b.at(i)});
fps ret(n+m-1);
for(int i=0; i<n; i++) for(auto [k,v] : P) ret.at(i+k) += (*this).at(i)*v;
return ret;
};
fps inv_sparse(int deg=-1)const{ //非0が少ない時の1/fを返す O(deg*|非0|+m+logmod).
int n = (*this).size();
if(deg == -1) deg = n;
assert((*this).at(0) != 0);
T div = 1;
vector<pair<int,T>> P;
if((*this).at(0) != 1) div = T((*this).at(0)).inv();
for(int i=1; i<n; i++) if((*this).at(i) != 0) P.emplace_back(pair{i,(*this).at(i)*div});
fps ret(deg); ret.at(0) = 1;
for(int i=1; i<deg; i++) for(auto [k,v] : P) if(i >= k) ret.at(i) -= v*ret.at(i-k); //-xが+ret[i-1]に対応.
if(div != 1) for(auto &v : ret) v *= div;
return ret;
}
fps inv_sparse(const fps &b,int deg = -1)const{ //f/gを返す 1/fでは分母だがこれは分子に注意.
int n = (*this).size(),m = b.size();
if(deg == -1) deg = n;
assert(b.at(0) != 0);
T div = 1;
vector<pair<int,T>> P;
if(b.at(0) != 1) div = T(b.at(0)).inv();
for(int i=1; i<m; i++) if(b.at(i) != 0) P.emplace_back(pair{i,b.at(i)*div});
fps ret = (*this).prefix(deg);
for(int i=1; i<deg; i++) for(auto [k,v] : P) if(i >= k) ret.at(i) -= v*ret.at(i-k);
if(div != 1) for(auto &v : ret) v *= div;
return ret;
}
fps log_sparse(int deg = -1){ //log(f)を返す O(N*非0).
//logf = ∫(f'/f) inv,1/f*(f')がO(N*非0) 他はO(N).
assert((*this).size()&&(*this).at(0)==1);
if(deg == -1) deg = (*this).size();
fps ret = (*this).diff();
ret = ((*this).inv_sparse(deg)).multi_sparse(ret);
return ret.inte().prefix(deg);
}
fps exp_sparse(int deg = -1)const{ //exp(f)を返す O(N*非0).
//(expf)'=(f')*exp(f)より低次から決まる.
//[x^0]expf=1より左辺のx^0の係数が求まる->expfのx^1の係数が求まる->...
if(deg == -1) deg = (*this).size();
fps ret(deg);
if((*this).size() == 0){
if(deg > 0) ret.at(0) = 1;
return ret;
}
assert((*this).at(0) == 0);
if(deg == 1) return fps{1};
ret.at(0) = 1; ret.at(1) = 1;
for(int i=2; i<deg; i++) ret.at(i) = (-ret.at(mod%i)*(mod/i));
vector<pair<int,T>> P;
for(int i=1; i<(*this).size(); i++) if((*this).at(i) != 0) P.emplace_back(pair{i-1,(*this).at(i)*i});
for(int i=0; i<deg-1; i++){
T now = 0;
for(auto [k,v] : P){
if(i < k) break;
now += ret.at(i-k)*v;
}
ret.at(i+1) *= now;
}
return ret;
}
fps pow_sparse(long long K,int deg = -1)const{ //f^Kを返す O(N*非0).
//f^K = Fとする fF'= K*F*f'を使う.
//[x^0]f=1に因数分解 (ax^i)^Kは別で処理.
//n次はΣ[i=0~n]fi*F'[n-i]=K*Σ[i=1~n+1]F[n+1-i]*i*fi.
//左辺のi=n以外移項 f0=1より F'n=K*Σ[i=1~n+1]F[n+1-i]*i*fi-Σ[i=1~n]fi*(n+1-i)*F[n+1-i].
//Fのn次まで分かっていればF'nが求まりFn+1も求まる F0=1からスタート.
int n = (*this).size();
if(deg == -1) deg = (*this).size();
if(K == 0 || deg == 0){
fps ret(deg);
if(deg > 0) ret.at(0) = 1;
return ret;
}
for(int t=0; t<n; t++){
if((*this).at(t) != 0){
T div = T(1)/(*this).at(t);
vector<pair<int,T>> P;
for(int i=t+1; i<n; i++) if((*this).at(i) != 0) P.emplace_back(pair{i-t,(*this).at(i)*div});
fps ret(deg); ret.at(0) = 1;
if(deg > 1) ret.at(1) = 1;
for(int i=2; i<deg; i++) ret.at(i) = (-ret.at(mod%i)*(mod/i));
T mulK = K%mod;
for(n=0; n<deg-1; n++){ //ここでnの役割が変わる注意.
T now = 0;
for(auto [i,v] : P){
if(i > n+1) break;
if(i > 0 && i<=n+1) now += mulK*ret.at(n+1-i)*i*v;
if(i > 0 && i <= n) now -= v*(n+1-i)*ret.at(n+1-i);
}
ret.at(n+1) *= now;
}
ret *= T((*this).at(t)).pow(K);
return (ret<<(t*K)).prefix(deg);
}
if(K >= deg || (t+1)*K >= deg) break;
}
return fps(deg,0);
}
fps diff()const{ //微分 nx^(n-1) = (x^n)' O(n).
int n = (*this).size();
if(n <= 1) return {};
fps ret(n-1);
T multi = 1;
for(int i=1; i<n; i++) ret.at(i-1) = (*this).at(i)*multi,multi++;
return ret;
}
fps inte()const{ //積分. 1/(n+1)x^(n+1) = ∫x^n dx O(n).
int n = (*this).size();
fps ret(n+1);
ret.at(0) = 0;
if(n == 0) return ret;
ret.at(1) = 1;
for(int i=2; i<=n; i++) ret.at(i) = (-ret.at(mod%i)*(mod/i)); //ここでret[i]=1/iになる.
for(int i=0; i<n; i++) ret.at(i+1) *= (*this).at(i);
return ret;
}
fps inv(int deg=-1)const{ //f*g=1 mod x^degとなるgを求める O(nlogn).
//g.pre(2m)<- (2g.pre(m)-f*g.pre(m)*g.pre(m))mod x^2mとする.
//fgm≡1 mod x^mからfg2m≡1 mod x^2mが求まる ダブリング.
//g2m<- gm-((fgm-1)*gm)を求める.
if((*this).size() <= 100) return inv_sparse(deg);
assert((*this).at(0) != 0);
if(deg == -1) deg = (*this).size();
fps ret(deg);
ret.at(0) = T(1)/(*this).at(0); T div2 = T(1)/2,divN = 1;
for(int i=1; i<deg; i<<=1){ //NTT-friendlyじゃないと駄目らしい.
fps f(2*i),g(2*i); divN *= div2;
int n = min(2*i,(int)(*this).size());
for(int k=0; k<n; k++) f.at(k) = (*this).at(k);
for(int k=0; k<i; k++) g.at(k) = ret.at(k);
NTT(f),NTT(g);
for(int k=0; k<2*i; k++) f.at(k) *= g.at(k); //(fgm-1)を計算.
INTT(f);
for(int k=0; k<i; k++) f.at(k) = 0; //[m,2m)の項だけにする([0,m)は一致してない).
for(int k=i; k<2*i; k++) f.at(k) *= divN; //convの後のdivNを忘れない.
NTT(f);
for(int k=0; k<2*i; k++) f.at(k) *= g.at(k); //(fgm-1)*gmを計算.
INTT(f);
n = min(2*i,deg);
for(int k=i; k<n; k++) f.at(k) *= divN;
for(int k=i; k<n; k++) ret.at(k) = -f.at(k); //-(fgm^2-gm)が入る.
}
return ret; //ね、簡単でしょ？.
}
fps inv2(int deg=-1)const{ //NTT-friendlyじゃなくてもいい代わりに定数が重い convも変える O(nlogn).
//g.pre(2m)<- (2g.pre(m)-f*g.pre(m)*g.pre(m))mod x^2mとする.
//fgm≡1 mod x^mからfg2m≡1 mod x^2mが求まる ダブリング.
//こっちは愚直にやる NTT-friendlyじゃないならllmodのconvolutionに変更しないと使えない.
if((*this).size() <= 200) return inv_sparse(deg);
assert((*this).at(0) != 0);
if(deg == -1) deg = (*this).size();
fps ret(1); ret.at(0) = T(1)/(*this).at(0);
for(int i=1; i<deg; i<<=1) ret = (ret+ret-ret*ret*((*this).prefix(i<<1))).prefix(i<<1);
return ret.prefix(deg);
}
fps log(int deg=-1,bool NTTfriendly = false)const{ //logfを求める O(nlogn).
//log(1-f) = -Σ[n=1~∞](f^n/n)と定義する [x^0]f=1&n<modが条件.
//(logf)'=f'/f,log(fg)=logf+loggが成り立つ.
assert((*this).size()&&(*this).at(0)==1);
if(deg == -1) deg = (*this).size();
if(NTTfriendly) ((*this).diff()*(*this).inv(deg)).prefix(deg-1).inte(); //∫f'/fから求める.
return ((*this).diff()*(*this).inv2(deg)).prefix(deg-1).inte();
}
fps exp(int deg=-1)const{ //expfを求める O(nlogn).
//expf=Σ[n=0~∞]f^n/n!と定める [x^0]f=0,n<modが条件.
//(expf)'= f'expf,exp(f+g)=expf*expg,log(expf)=fを満たす (各々条件は忘れずに).
//ニュートン法で求める g2m<-gm*(1-log(gm)+f) mod x^2nに更新 納得してない.
assert((*this).size() == 0 || (*this)[0] == mint(0));
if (deg == -1) deg = (*this).size();
vector<T> divi; //invと衝突回避用 mintでiの逆元.
divi.resize(deg*2); divi.at(1) = 1;
for(int i=2; i<deg*2; i++) divi.at(i) = ((-divi.at(mod%i)))*(mod/i);
auto integral = [&](fps &f) -> void { //inplaceで積分.
int n = f.size();
f.insert(f.begin(),0);
for(int i=1; i<=n; i++) f.at(i) *= divi.at(i);
};
auto differential = [&](fps &f) -> void { //inplaceで微分.
if(f.size() == 0) return;
f.erase(f.begin());
T multi = 0;
for(int i=0; i<f.size(); i++) f.at(i) *= ++multi;
};
fps f({1}),g({1}),z1,z2({1,1}); //何やってるか本当に意味わからずコメント書けない 後で確認.
if((*this).size() > 1) f.push_back((*this).at(1));
else f.push_back(0);
for(int m=2; m<deg; m<<=1){
auto ff = f;
ff.resize(2*m);
NTT(ff); z1 = z2;
fps z(m);
for(int i=0; i<m; i++) z.at(i) = ff.at(i)*z1.at(i);
INTT(z);
for(int i=0; i<m/2; i++) z.at(i) = 0;
for(int i=m/2; i<m; i++) z.at(i) *= divi.at(m);
NTT(z);
for(int i=0; i<m; i++) z.at(i) *= -z1.at(i);
INTT(z);
for(int i=m/2; i<m; i++) g.emplace_back(z.at(i)*divi.at(m));
z2 = g;
z2.resize(2*m);
NTT(z2);
fps h = (*this).prefix(m);
differential(h);
h.emplace_back(T(0));
NTT(h);
for(int i=0; i<m; i++) h.at(i) *= ff.at(i);
INTT(h);
for(int i=0; i<m; i++) h.at(i) *= divi.at(m);
h -= f.diff();
h.resize(2*m);
for(int i=0; i<m-1; i++) h.at(m+i) = h.at(i),h.at(i) = 0;
NTT(h);
for(int i=0; i<2*m; i++) h.at(i) *= z2.at(i);
INTT(h);
for(int i=0; i<2*m; i++) h.at(i) *= divi.at(2*m);
h.pop_back();
integral(h);
int n = min((int)(*this).size(),2*m);
for(int i=m; i<n; i++) h.at(i) += (*this).at(i);
for(int i=0; i<m; i++) h.at(i) = 0;
NTT(h);
for(int i=0; i<2*m; i++) h.at(i) *= ff.at(i);
INTT(h);
for(int i=0; i<2*m; i++) h.at(i) *= divi.at(2*m);
for(int i=m; i<2*m; i++) f.emplace_back(h.at(i));
}
return f.prefix(deg);
}
fps exp2(int deg=-1)const{ //定数重いやつinvと同じ O(nlogn).
//NTT-friendlyじゃないとK=10^5で1sec以上かかる 畳み込みが遅い.
assert((*this).size() == 0 || (*this).at(0) == 0);
if(deg == -1) deg = (*this).size();
fps ret(1); ret.at(0) = 1;
for(int i=1; i<deg; i<<=1){
//ニュートン法 g2m<-gm*(1-log(gm)+f) mod x^2n 再掲.
ret = (ret*((*this).prefix(2*i)+T(1)-ret.log(2*i,false))).prefix(2*i); //logの引数普段と違う.
}
return ret.prefix(deg);
}
fps pow(long long K,int deg=-1,bool NTTfriendly = true)const{ //K乗を返すO(nlogn) mod10^9+7で3番目引数false.
//f^k = exp(klog(f))を使う.
//logfが計算できるように調整する.
int n = (*this).size();
if(deg == -1) deg = n;
if(K == 0){ //0乗はx^0だけ係数1.
fps ret(deg);
if(deg > 0) ret.at(0) = 1;
return ret;
}
for(int i=0; i<n; i++){
if((*this).at(i) != 0){ //0じゃなかったらOK.
if(NTTfriendly){
T div = T(1)/(*this).at(i); //*([x^i]f)^Kと *x^(i*k)の分は後回し.
fps ret = ((((*this)*div)>>i).log(deg)*(K%mod)).exp(deg);
ret *= T((*this).at(i)).pow(K); //[x^i]f^Kの分.
ret = (ret<<(i*K)).prefix(deg); //*x^(i*k)の分.
return ret;
}
else{
T div = T(1)/(*this).at(i); //*([x^i]f)^Kと *x^(i*k)の分は後回し.
fps ret = ((((*this)*div)>>i).log(deg,false)*(K%mod)).exp2(deg);
ret *= T((*this).at(i)).pow(K); //[x^i]f^Kの分.
ret = (ret<<(i*K)).prefix(deg); //*x^(i*k)の分.
return ret;
}
}
if(K >= deg || (i+1)*K >= deg) break; //((i+1)*K)乗未満は0確定 int128回避用にK>=deg(degがllはやばい).
}
return fps(deg,0); //fの係数全て0なら係数全て0.
}
fps prefix(int siz)const{ //先頭siz項を返す なかったら0埋め. O(siz).
fps ret((*this).begin(),(*this).begin()+min((int)(*this).size(),siz));
if(ret.size() < siz) ret.resize(siz,0);
return ret;
}
void del0(){ //末尾の0を消す O(n).
while((*this).size() && (*this).back() == 0) (*this).pop_back();
}
fps rev()const{ //ひっくり返す O(n).
fps ret(*this);
reverse(ret.begin(),ret.end());
return ret;
}
pair<fps,fps> getQR(const fps &b)const{ //多項式の商と余りを同時に得る O(nlogn).
fps Q = (*this)/b,R = (*this)-Q*b;
R.del0();
return {Q,R};
}
fps cumulativeNtimes(int N,T b,int deg=-1){ //1/(1-bx)^Nをdeg次まで返す 指定なしはN次まで.
//負の二項定理を使う. 1/(1-x)^N=Σ[i=0~∞]((n+i-1) choose i)(bx^i);
//fps{}.cumulativeNtime()で無理やり関数を呼び出す.
assert(N <= 0); //N=0も駄目? {1}を返すべき所{0}になる.
if(deg == -1) deg = N+1;
int Limit = N+deg; //Limit 必要なサイズ fac->x! facinv->1/x! inv->1/x.
vector<T> FAC(Limit+1,1); for(int i=1; i<=Limit; i++) FAC.at(i) = FAC.at(i-1)*i;
vector<T> FACinv(Limit+1); FACinv.at(min((int)mod-1,Limit)) = FAC.at(min((int)mod-1,Limit)).inv();
for(int i=min((int)mod-2,Limit-1); i>=0; i--) FACinv.at(i) = FACinv.at(i+1)*(i+1);
auto nCr = [&](int n, int r) -> T {
if(n < r || r < 0 || n < 0) return 0;
return FAC.at(n)*FACinv.at(r)*FACinv.at(n-r);
};
fps ret(deg); T value = 1;
for(int i=0; i<deg; i++) ret.at(i) = nCr(N+i-1,i)*value,value *= b;
return ret;
}
fps Taylorshift(long long C)const{ //Σ[i=0~n]Ai*(x+c)^iを返す O(nlogn).
//Σ[i=0~n]x^i/i! Σ[k=0~n-i](A[i+k]*(i+k)!)*(c^k/k!)になる.
//kの部分はXi=Ai*i!,Yi=c^i/i!として.Zi=ΣX[i+k]Yk.
//これは[x^(i+n)]X*rev(Y)で求められる.
mint c = (C%mod+mod)%mod,pc = 1;
if(c == 0) return (*this);
int deg = (*this).size(); //deg 必要なサイズ fac->x! facinv->1/x! inv->1/x.
vector<mint> fac(deg+1,1); for(int i=1; i<=deg; i++) fac.at(i) = fac.at(i-1)*i;
vector<mint> facinv(deg+1); facinv.at(min((int)mod-1,deg)) = fac.at(min((int)mod-1,deg)).inv();
for(int i=min((int)mod-2,deg-1); i>=0; i--) facinv.at(i) = facinv.at(i+1)*(i+1);
fps X(deg),Y(deg);
for(int i=0; i<deg; i++) X.at(i) = (*this).at(i)*fac.at(i),Y.at(deg-1-i) = pc*facinv.at(i),pc *= c;
X *= Y;
fps ret(deg);
for(int i=0; i<deg; i++) ret.at(i) = X.at(i+deg-1)*facinv.at(i);
return ret;
}
void NTTdoubling(){ //NTTの原理を忘れているため何やってるのか意味が分からない NTT-friendly専用.
//INTT->resize(2倍）->NTTの代わりにcopy->INTT->謎の操作->NTT->push sizeが小さい時は効率悪いらしいよ.
int n = (*this).size();
assert(n > 0 && (n&-n) == n); //2べきじゃないとダメ.
mint div = mint(1)/n;
fps b = (*this);
INTT(b);
for(auto &v : b) v *= div;
mint rot = 1,zeta = mint(findroot()).pow((mod-1)/(n<<1));
for(auto &v : b) v *= rot,rot *= zeta;
NTT(b);
(*this).reserve(n<<1);
for(auto &v : b) (*this).push_back(v);
}
};
using fps = FormalPowerSeries<mint>;
//mod998244353以外のNTT-friendlyの時は色々気を付ける 後で書き直すかも?.
ostream &operator<<(ostream &os,const fps &a){
for(auto v : a) cout << v.v << " ";
return os;
}
vector<mint> MultiPointEvalution(const fps &f,const vector<mint> &X){ //M点で次数Nの多項式F(xi)の値をO(Mlog^2M+NlogN)で求める.
//剰余の定理よりF(a)=F(x)mod(x-a)
//F(x)mod(x-xi)の問題に帰着->F(x)mod(x-xi)を葉 F(x)mod(Π(x-xi))を根として考える.
//modの多項式を求めに葉から分割統治 逆に根から葉で求めたいものを求める.
int N = f.size(),M = X.size(),m2 = 1;
assert(N != 0 && M != 0);
if(M == 1){
mint px = 1,sum = 0;
for(auto v : f) sum += v*px,px *= X.at(0);
return {sum};
}
while(m2 < M) m2 <<= 1;
if(mod == 998244353){
vector<fps> Mods(m2*2);
for(int i=0; i<m2; i++){
mint v = 0;
if(i < M) v = -X.at(i);
Mods.at(i+m2) = fps{v+1,v-1};
}
for(int i=m2-1; i>0; i--){ //葉から構築.
fps &L = Mods.at(i<<1),&R = Mods.at((i<<1)+1);
int n = L.size();
Mods.at(i).resize(n);
for(int k=0; k<n; k++) Mods.at(i).at(k) = L.at(k)*R.at(k);
if(i == 1){
for(int k=0; k<n; k++) Mods.at(i).at(k)--;
break;
}
Mods.at(i).NTTdoubling(); n += n;
for(int k=n/2; k<n; k++) Mods.at(i).at(k) -= 2; //???.
}
fps root = Mods.at(1); //こっから意味不明.
INTT(root);
mint divn = mint(1)/(int)root.size();
for(auto &v : root) v *= divn;
root.push_back(1);
reverse(root.begin(),root.end());
root = root.inv(N).rev()*f;
root.erase(root.begin(),root.begin()+N-1);
root.resize(m2);
vector<mint> ret(M);
auto dfs = [&](auto dfs,int L,int R,int pos,fps g,mint div) -> void {
if(pos >= m2){ret.at(pos-m2) = g.at(0); return;}
pos <<= 1;
int len = g.size(),mid = (L+R)>>1;
NTT(g);
fps &LM = Mods.at(pos+1);
for(int i=0; i<len; i++) LM.at(i) *= g.at(i);
INTT(LM);
for(int i=len>>1; i<len; i++) LM.at(i) *= div;
dfs(dfs,L,mid,pos,fps{LM.begin()+(len>>1),LM.end()},div+div);
if(mid >= M) return;
fps &RM = Mods.at(pos);
for(int i=0; i<len; i++) RM.at(i) *= g.at(i);
INTT(RM);
for(int i=len>>1; i<len; i++) RM.at(i) *= div;
dfs(dfs,mid,R,pos+1,fps{RM.begin()+(len>>1),RM.end()},div+div);
};
dfs(dfs,0,m2,1,root,divn);
return ret;
}
else{
vector<fps> Mods(m2*2,{1}),F(m2*2);
for(int i=m2; i<m2+M; i++) Mods.at(i) = {-X.at(i-m2),1};
for(int i=m2-1; i>0; i--) Mods.at(i) = Mods.at((i<<1))*Mods.at((i<<1)+1);
F.at(1) = f%Mods.at(1);
for(int i=1; i<m2; i++){
int pos = i<<1;
F.at(pos) = F.at(i)%Mods.at(pos);
F.at(pos+1) = F.at(i)%Mods.at(pos+1);
}
vector<mint> ret(M);
for(int i=0; i<M; i++) if(F.at(i+m2).size()) ret.at(i) = F.at(i+m2).at(0);
return ret;
}
}
fps Productfps(vector<fps> Fs){ //ΠFiを返す 返り値の次数をdをしてO(dlog^2d).
int n = Fs.size();
if(n == 0) return {1};
if(n == 1) return Fs.at(0);
auto comp = [&](auto a,auto b){return Fs.at(a).size()>Fs.at(b).size();};
priority_queue<int,vector<int>,decltype(comp)> PQ(comp);
for(int i=0; i<n; i++) PQ.push(i);
while(PQ.size() > 1){
int p = PQ.top(); PQ.pop();
int q = PQ.top(); PQ.pop();
Fs.at(p) *= Fs.at(q); PQ.push(p);
//Fs.at(q).clear();
}
return Fs.at(PQ.top());
}
vector<mint> Productmint(vector<vector<mint>> Fs){ //ΠFiを返す 返り値の次数をdをしてO(dlog^2d).
int n = Fs.size();
if(n == 0) return {1};
if(n == 1) return Fs.at(0);
auto comp = [&](auto a,auto b){return Fs.at(a).size()>Fs.at(b).size();};
priority_queue<int,vector<int>,decltype(comp)> PQ(comp);
for(int i=0; i<n; i++) PQ.push(i);
while(PQ.size() > 1){
int p = PQ.top(); PQ.pop();
int q = PQ.top(); PQ.pop();
Fs.at(p) = convolution(Fs.at(p),Fs.at(q)); PQ.push(p);
//Fs.at(q).clear();
}
return Fs.at(PQ.top());
}
int main(){
ios_base::sync_with_stdio(false);
cin.tie(nullptr);
int M,N; cin >> M >> N;
int Limit = N; //Limit 必要なサイズ fac->x! facinv->1/x! inv->1/x.
vector<mint> fac(Limit+1,1); for(int i=1; i<=Limit; i++) fac.at(i) = fac.at(i-1)*i;
vector<mint> facinv(Limit+1); facinv.at(min((int)mod-1,Limit)) = fac.at(min((int)mod-1,Limit)).inv();
for(int i=min((int)mod-2,Limit-1); i>=0; i--) facinv.at(i) = facinv.at(i+1)*(i+1);
//vector<mint> inv(Limit+1); for(int i=1; i<=Limit; i++) inv.at(i) = fac.at(i-1)*facinv.at(i); //必要なら解放!.
auto nCr = [&](int n, int r) -> mint {
if(n < r || r < 0 || n < 0) return 0;
return fac.at(n)*facinv.at(r)*facinv.at(n-r);
};
auto nPr = [&](int n, int r) -> mint {
if(n < r || r < 0 || n < 0) return 0;
return fac.at(n)*facinv.at(n-r);
};
vector<mint> C(M);
for(int i=0; i<N; i++){
int a; cin >> a; a--;
C.at(a)++;
}
vector<vector<mint>> Bs;
for(auto c : C){
for(int i=0; i<c.v; i++) Bs.emplace_back(vector<mint>{-i,1});
}
vector<mint> A = facinv;
for(int i=1; i<=N; i+=2) A.at(i) = -A.at(i);
auto F = Productmint(Bs);
fps f; for(auto v : F) f.emplace_back(v);
vector<mint> X(N+1);
for(int i=0; i<=N; i++) X.at(i) = i;
auto B = MultiPointEvalution(f,X);
for(int i=0; i<=N; i++) B.at(i) *= facinv.at(i);
C = convolution(A,B);
mint answer = 0;
for(int i=0; i<=N; i++) answer += fac.at(i)*C.at(i);
cout << answer.v << endl;
}