方便起见给 (k) 减 (1)，考虑答案为：
[egin{aligned}
&Ans_k=sum f_{i,k}inom{n}{i}(n-i)!
end{aligned}]
对于 (f_{i,k}) 考虑通过容斥计算，设 (g_{i,k}) 表示长度为 (i) 的序列存在 (k) 处 < 的方案数。
那么就有：
[g_{i,k}=sum f_{i,j}inom{j}{k}
]
对于 (g_{i,k})，考虑其代表了 (i-k) 个连通块，此时相当于将 (i) 个带标号球分配给 (i-k) 个盒子的方案数，又相当于给 (1sim i) 染上 (i-k) 种颜色，且每类颜色非空的方案数，即 (i!cdot (e^x-1)^{i-k}[x^i])
于是：
[f_{i,k}=sum_j inom{j}{k}(-1)^{j-k}g_{i,j}
]
故：
[Ans_k=sum_i inom{n}{i}(n-i)!sum_j inom{j}{k}(-1)^{j-k}i!cdot (e^x-1)^{i-j}[x^i]
]
[Ans_k=n!sum_i sum_j inom{j}{k}(-1)^{j-k}cdot (e^x-1)^{i-j}[x^i]
]
[Ans_k=n!sum_j inom{j}{k}(-1)^{j-k}sum_{ige j}(e^x-1)^{i-j}[x^i]
]
设后面那个东西为 (f_j)，此时前半部分就是差值卷积。
- (frac{n!}{k!}sum_j j!frac{(-1)^{j-k}}{(j-k)!}f_j)
只需要考虑计算 (f_i)
为此考虑：
[f_j=sum_{ige j}(e^x-1)^{i-j}[x^i]
]
并考虑：
[[x^j]sum_{i=0}^{n-j}(frac{e^x-1}{x})^i
]
设为 (F(x))，即：
[[x^j]frac{1-F(x)^{n-j+1}}{1-F(x)}
]
分成两个 Part：
- Part1，计算 (frac{1}{1-F(x)})，由于常数项为 (0) 无法直接求逆，引入负数幂的定义后表示为 (x^kF^{}(x) o x^{-k}F^{-1}(x))
- Part2：
[[x^j]frac{F(x)^{n-j+1}}{1-F(x)} o [x^{n+1}]frac{(xF(x))^{n-j+1}}{1-F(x)}
]
由于需要对于每个 (j) 进行计算，这里考虑加入一个元以提取信息，我们转为求：
[[x^{n+1}]sum_j frac{(xyF(x))^{n-j+1}}{1-F(x)}
]
[[x^{n+1}][y^{n-j+1}]sum_j frac{(xyF(x))^{j}}{1-F(x)} o frac{1}{1-F(x)} imes frac{1}{1-xyF(x)}
]
我们希望将其表述为多项式复合的形式，这样我们才可以通过拉格朗日反演来描述答案，事实上，我们可以设：
(W(x)=xF(x),H(x)) 为一个多项式满足 (H(W(x))=F(x))，并设 (G(x)=frac{1}{1-H(x)}frac{1}{1-xy})
这样我们就有 (G(W(x))=frac{1}{1-H(W(x))} imes frac{1}{1-W(x)y})
于是我们所求就是：
[[x^{n+1}]G(W(x))=frac{1}{n+1}G(x)'P(x)^{-(n+1)}[x^{-1}]
]
其中 (P(x)) 为 (W(x)) 的复合逆，因此，答案即：
[frac{1}{n+1}G(x)'Big(frac{x}{P(x)}Big)^{n+1}[x^n]
]
我们先求出 (G(x)) 的导数：
[frac{y(1-H(x))+(1-xy)H'(x)}{(1-H(x))^2(1-xy)^2}
]
此时答案即：
[egin{aligned}
&frac{1}{n+1}G(x)'Big(frac{x}{P(x)}Big)^{n+1}[x^n]
\&frac{1}{n+1}frac{y(1-H(x))+(1-xy)H'(x)}{(1-H(x))^2(1-xy)^2}Big(frac{x}{P(x)}Big)^{n+1}[x^n]
end{aligned}]
另一边，求复合逆一直是一个非常困难的问题，但事实上考虑到 (H(W(x))=F(x))，考虑 (P(x)=frac{x}{H(x)} o P(W(x))=frac{W(x)}{H(W(x))}=x)
于是我们等价于求：
[egin{aligned}
\&[x^n]frac{1}{n+1}frac{y(1-H(x))+(1-xy)H'(x)}{(1-H(x))^2(1-xy)^2}Big(H(x)Big)^{n+1}
\&[x^n]frac{1}{n+1} imes igg(yfrac{1}{(1-H(x))(1-xy)^2}+frac{H'(x)}{(1-H(x))^2(1-xy)}igg)H(x)^{n+1}
\&[x^n]frac{1}{n+1} imes Bigg(frac{1}{1-H(x)}sum (i+1)x^iy^{i+1}+frac{H'(x)}{(1-H(x))^2}sum x^iy^iBigg)H(x)^{n+1}
\&[x^ny^{j}]frac{1}{n+1} imes Bigg(jfrac{H(x)^{n+1}x^{j-1}}{1-H(x)}+frac{x^jH'(x)H(x)^{n+1}}{(1-H(x))^2}Bigg)
end{aligned}]
于是答案即：
[frac{1}{n+1} imes Bigg([x^{n-j+1}]jfrac{H(x)^{n+1}}{1-H(x)}+[x^{n-j}]frac{H'(x)H(x)^{n+1}}{(1-H(x))^2}Bigg)
]
注意 (f_j=[x^n][y^{n-j+1}]) 即：
[frac{1}{n+1} imes Bigg([x^{j}](n-j+1)frac{H(x)^{n+1}}{1-H(x)}+[x^{j+1}]frac{H'(x)H(x)^{n+1}}{(1-H(x))^2}Bigg)
]
最后，我们的瓶颈在于求解 (H(x))
注意到 (H(x)=frac{x}{P(x)})，考虑 (P(x)) 为 (W(x)) 即 (e^x-1) 的复合逆。
于是 (P(x)=ln(x+1))，于是 (H(x)=frac{x}{ln(x+1)})
特别的，计算得到的 (H(x)[x^0]) 仍为 (1)，这意味着 (1-H(x)) 是需要进行偏移以保证存在逆元的。
最终复杂度 (mathcal O(nlog n))
注意偏移量。
(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
//视题目调节
const int P = 998244353 ;
const int Gi = 332748118 ;
const int G = 3 ;
// 基本参数
const int N = 8e5 + 5 ;
int gi() {
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 ;
}
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 ;
}
// 求逆
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 ; a[x] = 0 ;
}
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 n, ret, dh[N], h[N], H[N], dH[N], IH[N], IH2[N] ;
int f[N], ans[N], fac[N], inv[N] ;
signed main()
{
n = gi(), ret = n + 10, fac[0] = inv[0] = 1 ;
Poly::_init(ret * 2) ;
rep( i, 1, ret ) fac[i] = fac[i - 1] * i % P ;
rep( i, 0, ret ) inv[i] = fpow(fac[i], P - 2) ;
h[0] = h[1] = 1, Poly::Polyln(h, ret) ;
rep( i, 0, ret ) h[i] = h[i + 1] ;
Poly::PolyInv(h, ret) ;
rep( i, 0, ret ) IH2[i] = IH[i] = dH[i] = H[i] = h[i] ;
Poly::PolyKSM(H, ret, n + 1), Poly::Deriv(dH, ret) ;
rep( i, 0, ret ) IH[i] = IH2[i] = (P - h[i]) % P ;
IH[0] = IH2[0] = 0 ;
Poly::Init(ret + ret), Poly::NTT(IH2, 1) ;
rep( i, 0, Poly::limit ) IH2[i] = IH2[i] * IH2[i] % P ;
Poly::NTT(IH2, 0) ;
rep( i, ret + 1, Poly::limit ) IH2[i] = 0 ;
rep( i, 0, ret ) IH2[i] = IH2[i + 2], IH[i] = IH[i + 1] ;
Poly::PolyInv(IH2, ret), Poly::PolyInv(IH, ret) ;
rep( i, 0, ret ) ans[i] = (P - inv[i + 2]) % P ;
Poly::PolyInv(ans, ret) ;
rep( i, 1, n ) f[i] = ans[i + 1] ;
//----------------------------------------------
Poly::Init(ret * 4) ;
Poly::NTT(IH2, 1), Poly::NTT(H, 1),
Poly::NTT(dH, 1), Poly::NTT(IH, 1) ;
rep( i, 0, Poly::limit )
dh[i] = H[i] * IH[i] % P,
IH2[i] = IH2[i] * dH[i] % P * H[i] % P,
IH[i] = IH[i] * H[i] % P ;
Poly::NTT(IH, 0), Poly::NTT(IH2, 0) ;
//---------------------------------------------
int iv = fpow(n + 1, P - 2) ;
rep( i, 1, n ) f[i] = (f[i] - iv * ((n - i + 1) * IH[i + 1] % P + IH2[i + 1]) % P + P) % P ;
f[0] = n ;
rep( i, 0, ret ) h[i] = 0 ;
rep( i, 0, n ) h[i] = ((n - i) & 1) ? P - inv[n - i] : inv[n - i] ;
rep( i, 0, n ) f[i] = f[i] * fac[i] % P ;
Poly::Init(n * 2 + 2), Poly::NTT(f, 1), Poly::NTT(h, 1) ;
rep( i, 0, Poly::limit ) h[i] = f[i] * h[i] % P ;
Poly::NTT(h, 0) ;
rep( i, 0, n ) f[i] = h[i + n] ;
rep( i, 0, n ) f[i] = f[i] * fac[n] % P * inv[i] % P ;
rep( i, 1, n ) printf("%lld ", f[i - 1] ) ;
return 0 ;
}