第一次出现合法状态的期望时间可以转换为所有非法状态的出现概率乘以其在此处期望停留的时间。
设此非法状态抽了 \(r\) 张卡，那么其有 \(\binom{m}{r}^{-1}\) 的概率出现，经过 \(\frac{m}{m-r}\) 时间后会到达下一个状态。
只需要计数所有非法状态的数量即可。
对每个极长连续段考虑，设其长度为 \(n\)，考虑其中选出 \(r\) 张牌的方案数应为：
\[[x^{n+1-r}]\bigg(\sum x^i[i<k]\bigg)^{r} \]这一部分基于这样的一个考虑，我们添加一个辅助元 \(n+1\) 并钦定其不被选中，这样的话我们枚举没有被选择的元素的个数，并为其前面添增 \(0,1...k-1\) 个被选中的元素，这样只需要将计算得到的答案反转系数即可得到选中了 \(r\) 个元素的方案数。
因此，我们只需要考虑：
\[[x^{n+1-r}]\bigg(\sum x^i[i<k]\bigg)^{r} \]接下来我们将结尾段的 \(0\) 加入计数的段落中，并添加辅助元 \(u\) 表示被选择的 \(0\) 的数量，设 \(G(x)=\sum_{i\ge 1}x^i[i\le k]=x\frac{1-x^k}{1-x}\)
则我们希望求解的就是：
\[\sum (uG(x))^i=\frac{1}{1-uG(x)}[x^{n+1}] \]设 \(H(x)=\frac{1}{1-ux}\)，则所求为 \(H(G(x))[x^{n+1}]\)，设 \(F(x)\) 为 \(G(x)\) 的复合逆。
我们直接上拓展拉格朗日反演：
\[\frac{1}{1-uG(x)}[x^{n+1}]=\frac{1}{n+1}\frac{u}{(1-ux)^2}F(x)^{-(n+1)}[x^{-1}] \]也即：
\[\frac{1}{n+1}\frac{u}{(1-ux)^2}\Big(\frac{x}{F(x)}\Big)^{(n+1)}[x^{n}] \]事实上，由于 \(G(x)=\frac{x-x^{k+1}}{1-x}\)，所以 \(F(x)\) 会满足：
\[\frac{F(x)^{k+1}-F(x)}{F(x)-1}=x\to F(x)^{k+1}-(1+x)F(x)+x=0 \]牛迭部分的方程如下：
\[F(x)=F_0(x)-\frac{x^{k}F_0(x)^{k+1}-(1+x)F_0(x)+1}{x^{k}(k+1)F_0(x)^{k}-(1+x)} \]上面的 \(F(x)\) 是 \(\frac{F(x)}{x}\)。
\(Code:\)
#include<bits/stdc++.h>
using namespace std ;
#define drep( i, s, t ) for( register int i = (t); i >= (s); -- i )
#define rep( i, s, t ) for( register int i = (s); i <= (t); ++ i )
#define re register
#define int long long
#define vi vector<int>
//视题目调节
const int P = 998244353 ;
const int Gi = 332748118 ;
const int G = 3 ;
// 基本参数
const int N = 2e6 + 5 ;
int n, m, len, B[N], A[N], Dx[N] ;
vector<int> f[N] ;
int read() {
char cc = getchar() ; int cn = 0, flus = 1 ;
while( cc < '0' || cc > '9' ) { if( cc == '-' ) flus = - flus ; cc = getchar() ; }
while( cc >= '0' && cc <= '9' ) cn = ( cn * 10 + cc - '0' ) % P, cc = getchar() ;
return cn * flus ;
}
int fpow( int x, int k ) {
int ans = 1, base = x ;
while( k ) {
if( k & 1 ) ans = ans * base % P ;
base = base * base % P, k >>= 1 ;
} return ans % P ;
}
void inc(int &x, int y) {
((x += y) >= P) && (x -= P) ;
}
void dec(int &x, int y) {
((x -= y) < 0) && (x += P) ;
}
namespace Poly {
int limit, L, Inv, R[N], inv[N], gi[N << 1], igi[N << 1] ;
void Init( int x ) {
limit = 1, L = 0 ; while( limit < x ) limit <<= 1, ++ L ;
rep( i, 0, limit ) R[i] = ( R[i >> 1] >> 1 ) | ( ( i & 1 ) << ( L - 1 ) ) ;
Inv = fpow( limit, P - 2 ) ;
}
void _init( int x ) { // max x range
for( re int i = 0; i <= x; ++ i ) inv[i] = fpow( i, P - 2 ) ;
Init( x + x ) ; int num = 0 ;
for( re int k = 1; k < limit; k <<= 1 ) {
int d = fpow( G, ( P - 1 ) / ( k << 1 ) ) ;
for( re int j = 0, g = 1; j < k; ++ j, g = g * d % P ) gi[++ num] = g ;
} num = 0 ;
for( re int k = 1; k < limit; k <<= 1 ) {
int d = fpow( Gi, ( P - 1 ) / ( k << 1 ) ) ;
for( re int j = 0, g = 1; j < k; ++ j, g = g * d % P ) igi[++ num] = g ;
}
}
void NTT( int *a, int type ) { //NTT多项式乘法
for( re int i = 0; i < limit; ++ i ) if( i < R[i] ) swap( a[i], a[R[i]] ) ;
int num = 0 ;
for( re int k = 1; k < limit; k <<= 1 ) {
for( re int i = 0; i < limit; i += ( k << 1 ) )
for( re int j = i, d = num + 1; j < i + k; ++ j, ++ d ) {
int Nx = a[j], Ny = a[j + k] * ( ( type == 1 ) ? gi[d] : igi[d] ) % P ;
a[j] = ( Nx + Ny ) % P, a[j + k] = ( Nx - Ny + P ) % P ;
} num += k ;
}
if( type != 1 ) rep( i, 0, limit ) a[i] = a[i] * Inv % P ;
}
void NTT( vi& a, int type ) { //NTT多项式乘法
for( re int i = 0; i < limit; ++ i ) if( i < R[i] ) swap( a[i], a[R[i]] ) ;
int num = 0 ;
for( re int k = 1; k < limit; k <<= 1 ) {
for( re int i = 0; i < limit; i += ( k << 1 ) )
for( re int j = i, d = num + 1; j < i + k; ++ j, ++ d ) {
int Nx = a[j], Ny = a[j + k] * ( ( type == 1 ) ? gi[d] : igi[d] ) % P ;
a[j] = ( Nx + Ny ) % P, a[j + k] = ( Nx - Ny + P ) % P ;
} num += k ;
}
if( type != 1 ) rep( i, 0, limit ) a[i] = a[i] * Inv % P ;
}
void Mul(vi& a, vi& b) {
int l1 = a.size(), l2 = b.size() ;
Init(l1 + l2), a.resize(limit + 1), b.resize(limit + 1) ;
NTT(a, 1), NTT(b, 1) ;
rep( i, 0, limit ) a[i] = a[i] * b[i] % P ;
NTT(a, 0) ; int len = l1 + l2 ;
rep( i, len, limit ) a[i] = 0 ;
a.resize(len) ;
}
// 求逆
int H[N], Gx[N] ;
void PolyInv( int *a, int x ) {
Gx[0] = fpow( a[0], P - 2 ) ; int lim = 1 ;
while( lim < x ) {
lim <<= 1, Init( lim << 1 ) ;
for( re int i = 0; i < lim; ++ i ) H[i] = a[i] ;
NTT( H, 1 ), NTT( Gx, 1 ) ;
for( re int i = 0; i < limit; ++ i ) Gx[i] = ( 2ll - Gx[i] * H[i] % P + P ) % P * Gx[i] % P ;
NTT( Gx, 0 ) ; for( re int i = lim; i < limit; ++ i ) Gx[i] = 0 ;
}
for( re int i = 0; i <= x; ++ i ) a[i] = Gx[i] ;
rep( i, 0, limit ) Gx[i] = H[i] = 0 ;
}
void Deriv( int *a, int x ) { //求导
rep( i, 1, x ) a[i - 1] = a[i] * i % P ;
}
void Integ( int *a, int x ) { //积分
drep( i, 1, x ) a[i] = a[i - 1] * inv[i] % P ;
a[0] = 0 ;
}
int Gt[N], F[N] ;
void Polyln( int *a, int x ) {
rep( i, 0, x ) F[i] = Gt[i] = a[i] ;
PolyInv( Gt, x ), Deriv( F, x ) ;
Init( x << 1 ), NTT( Gt, 1 ), NTT( F, 1 ) ;
rep( i, 0, limit ) F[i] = F[i] * Gt[i] % P ;
NTT( F, 0 ), Integ( F, x ) ;
rep( i, 0, x - 1 ) a[i] = F[i] ;
rep( i, 0, limit ) F[i] = Gt[i] = 0 ;
}
int _Gt[N], _F[N], _H[N] ;
void Polyexp( int *a, int x ) {
_Gt[0] = 1 ; int lim = 1 ;
while( lim < x ) {
lim <<= 1 ; rep( i, 0, lim - 1 ) _F[i] = a[i], _H[i] = _Gt[i] ;
Polyln( _H, lim ), Init( lim << 1 ) ;
NTT( _Gt, 1 ), NTT( _H, 1 ), NTT( _F, 1 ) ;
rep( i, 0, limit ) _Gt[i] = ( 1ll - _H[i] + _F[i] + P ) % P * _Gt[i] % P ;
NTT( _Gt, 0 ) ; rep( i, lim, limit ) _Gt[i] = 0 ;
}
rep( i, 0, x - 1 ) a[i] = _Gt[i] ;
rep( i, 0, limit ) _F[i] = _H[i] = _Gt[i] = 0 ;
}
void PolyKSM( int *a, int x, int k ) {
Polyln( a, x ) ;
rep( i, 0, x - 1 ) a[i] = a[i] * k % P ;
Polyexp( a, x ) ;
}
int _Nw[N], _Gx[N], _Gy[N], _Gz[N], _Go[N], _GG[N] ;
void Newton(int *a, int x) {
_Nw[0] = 1 ; int lim = 1 ;
while( lim < x ) {
lim <<= 1 ; rep( i, 0, lim - 1 ) _Gz[i] = _Gx[i] = _Gy[i] = _Nw[i] ;
PolyKSM(_Gx, lim, m) ;
rep( i, 0, lim ) _GG[i] = _Gx[i] ;
Init(lim << 1) ;
NTT(_GG, 1), NTT(_Gy, 1) ;
rep( i, 0, limit ) _Gy[i] = _Gy[i] * _GG[i] % P ;
NTT(_Gy, 0) ;
rep( i, lim, limit ) _Gy[i] = 0 ;
rep( i, 0, lim ) _Gx[i] = _Gx[i] * (m + 1) % P ;
rep( i, 0, lim ) _Go[i] = _Gx[i], _Gx[i] = 0 ;
rep( i, m, lim - 1 ) _Gx[i] = _Go[i - m] ;
dec(_Gx[0], 1), dec(_Gx[1], 1), PolyInv(_Gx, lim) ;
for(re int i = lim - 1; i; -- i)
inc( _Gz[i], _Gz[i - 1] ) ;
rep( i, 0, lim ) _Go[i] = _Gy[i], _Gy[i] = 0 ;
rep( i, m, lim - 1 ) _Gy[i] = _Go[i - m] ;
rep( i, 0, lim ) dec(_Gy[i], _Gz[i]) ;
inc(_Gy[0], 1), Init(lim << 1) ;
NTT(_Gx, 1), NTT(_Gy, 1) ;
rep( i, 0, limit ) _Gx[i] = _Gx[i] * _Gy[i] % P ;
NTT(_Gx, 0) ;
rep( i, lim, limit ) _Gx[i] = 0 ;
rep( i, 0, lim - 1 ) dec(_Nw[i], _Gx[i]) ;
}
rep( i, 0, x - 1 ) a[i] = _Nw[i] ;
}
}
int top, st[N], fac2[N], inv2[N] ;
void solve(int l, int r) {
if(l == r) return ;
int mid = (l + r) >> 1 ;
solve(l, mid), solve(mid + 1, r) ;
Poly::Mul(f[l], f[mid + 1]) ;
}
int C(int x, int y) {
return 1ll * fac2[x] * inv2[y] % P * inv2[x - y] % P ;
}
signed main()
{
n = read(), m = read(), fac2[0] = inv2[0] = 1 ;
rep( i, 1, n ) B[i] = read() ;
rep( i, 1, n ) fac2[i] = fac2[i - 1] * i % P ;
rep( i, 0, n ) inv2[i] = fpow(fac2[i], P - 2) ;
Poly::_init(n * 2), Poly::Newton(A, n) ;
Poly::PolyInv(A, n) ;
sort(B + 1, B + n + 1), B[0] = -1 ; int l = 0 ;
rep( i, 1, n ) {
if(B[i] != B[i - 1] + 1) st[++ top] = l, l = 0 ; ++ l ;
} st[++ top] = l ;
for(re int i = 1; i <= top; ++ i) {
int u = st[i], iv = fpow(u + 1, P - 2) ;
f[i].resize(u + 1) ;
rep( j, 0, u ) Dx[j] = A[j] ;
Poly::PolyKSM(Dx, u + 1, u + 1) ;
rep( j, 1, u ) f[i][j] = j * Dx[u - j + 1] % P * iv % P ;
rep( j, 1, u ) if((u - j + 1) > j) swap(f[i][j], f[i][u - j + 1]) ;
f[i][0] = 1 ;
}
solve(1, top) ;
int Ans = 0 ;
rep( i, 1, n - 1 ) Ans = (Ans + f[1][i] * fpow(C(n, i), P - 2) % P * n % P * fpow(n - i, P - 2)) % P ;
cout << (Ans + 1) % P << endl ;
return 0 ;
}