方便起见给 \(k\) 减 \(1\),考虑答案为:
\[\begin{aligned} &Ans_k=\sum f_{i,k}\binom{n}{i}(n-i)! \end{aligned}\]对于 \(f_{i,k}\) 考虑通过容斥计算,设 \(g_{i,k}\) 表示长度为 \(i\) 的序列存在 \(k\) 处 <
的方案数。
那么就有:
\[g_{i,k}=\sum f_{i,j}\binom{j}{k} \]对于 \(g_{i,k}\),考虑其代表了 \(i-k\) 个连通块,此时相当于将 \(i\) 个带标号球分配给 \(i-k\) 个盒子的方案数,又相当于给 \(1\sim i\) 染上 \(i-k\) 种颜色,且每类颜色非空的方案数,即 \(i!\cdot (e^x-1)^{i-k}[x^i]\)
于是:
\[f_{i,k}=\sum_j \binom{j}{k}(-1)^{j-k}g_{i,j} \]故:
\[Ans_k=\sum_i \binom{n}{i}(n-i)!\sum_j \binom{j}{k}(-1)^{j-k}i!\cdot (e^x-1)^{i-j}[x^i] \] \[Ans_k=n!\sum_i \sum_j \binom{j}{k}(-1)^{j-k}\cdot (e^x-1)^{i-j}[x^i] \] \[Ans_k=n!\sum_j \binom{j}{k}(-1)^{j-k}\sum_{i\ge 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_{i\ge 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)\to x^{-k}F^{-1}(x)\)
- Part2:
由于需要对于每个 \(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)}\to \frac{1}{1-F(x)}\times \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))}\times \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} \]此时答案即:
\[\begin{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)}\to P(W(x))=\frac{W(x)}{H(W(x))}=x\)
于是我们等价于求:
\[\begin{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}\times \bigg(y\frac{1}{(1-H(x))(1-xy)^2}+\frac{H'(x)}{(1-H(x))^2(1-xy)}\bigg)H(x)^{n+1} \\&[x^n]\frac{1}{n+1}\times \Bigg(\frac{1}{1-H(x)}\sum (i+1)x^iy^{i+1}+\frac{H'(x)}{(1-H(x))^2}\sum x^iy^i\Bigg)H(x)^{n+1} \\&[x^ny^{j}]\frac{1}{n+1}\times \Bigg(j\frac{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}\times \Bigg([x^{n-j+1}]j\frac{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}\times \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(n\log 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 ;
}