如果这篇博客帮助到你，可以请我喝一杯咖啡~
CC BY-NC-SA 4.0 （除特别声明或转载文章外）
好神仙啊。
对于任何一个好序列 $a$，可以通过这样的方式把它转化为一个排列 $p$：先倒着列出 1 的出现位置，然后倒着列出 2 的出现位置，如此类推。
我们发现，一个排列和一个好的序列一一对应，第 $p_i$ 个位置的值为 $p_1,\cdots,p_i$ 间的升高数加一。
设答案为 $f_{0,\cdots,n-1}$，则有：
\[f_k=\sum\limits_{i=1}^n\left<\begin{matrix}i\\k\end{matrix}\right>\begin{pmatrix}n\\i\end{pmatrix}(n-i)!\]这样 F1 就做完了。
我们注意到复杂度瓶颈在欧拉数的计算，至少有 $k$ 次升高的方案数为 $n![x^n] (e^x-1)^{n-k}$（钦定 $k$ 个位置填 $<$，小于号连接的部分称为一段，则总共有 $n-k$ 段，一段的 $\text{EGF}$ 为 $e^x-1$，即一个非空集合），那么欧拉数就有
\[\left<\begin{matrix}n\\k\end{matrix}\right>=\sum\limits_{i=k}^n\begin{pmatrix}i\\k\end{pmatrix}(-1)^{i-k}n^{n-i}\]则有
\[\begin{aligned} f_k=&\sum\limits_{i=1}^n\sum\limits_{j=k}^i\begin{pmatrix}j\\k\end{pmatrix}(-1)^{j-k}i^{i-j}\begin{pmatrix}n\\i\end{pmatrix}(n-i)!\\ =&\frac{n!}{k!}\sum\limits_{i=1}^n\sum\limits_{j=k}^i\frac{j!(-1)^{j-k}}{(j-k)!}[x^i](e^x-1)^{i-j}\\ =&\frac{n!}{k!}\sum\limits_{j=k}^n\frac{j!(-1)^{j-k}}{(j-k)!}\sum\limits_{i=j}^n[x^i](e^x-1)^{i-j}\\ \end{aligned}\]注意到如果设 $g_j=\sum\limits_{i=j}^n[x^i] (e^x-1)^{i-j}$，上面这个东西就很好计算了。
这样复杂度瓶颈在求出 $g$，继续推式子。
\[\begin{aligned} g_j&=\sum\limits_{i=j}^n[x^i](e^x-1)^{i-j}\\ &=[x^j]\sum\limits_{i=j}^n(\frac{e^x-1}{x})^{i-j}\\ &=[x^j]\sum\limits_{i=0}^{n-j}(\frac{e^x-1}{x})^i\\ &=[x^j]\frac{1-(\frac{e^x-1}{x})^{n-j+1}}{1-\frac{e^x-1}{x}}\\ &=[x^j]\frac{1}{1-\frac{e^x-1}{x}}-[x^j]\frac{(\frac{e^x-1}{x})^{n-j+1}}{1-\frac{e^x-1}{x}}\\ \end{aligned}\]下面使用 $F(x)=\frac{e^x-1}{x}$。
第一部分直接多项式求逆，$[x^j]\frac{1}{1-F(x)}=[x^j]\frac{x^{-1}}{x^{-1}(1-F(x))}=[x^{j+1}]\frac{1}{x^{-1}(1-F(x))}$。
后面部分 $[x^j]\frac{F(X)^{n-j+1}}{1-F(x)}=[x^{n+1}]\frac{(xF(x))^{n-j+1}}{1-F(x)}$。
我们用一种神仙做法——加一个元以区分信息，那么这个东西就等价于 $[x^{n+1}y^{n-j+1}]\sum\limits_{i=0}^\infin\frac{(xF(x)y)^i}{1-F(x)}=[x^{n+1}y^{n-j+1}] (\frac{1}{1-F(x)}\frac{1}{1-xF(x)y})$。
设 $W(x)=xF(x)=e^x-1$，$H(x)$ 满足 $\frac{W(x)}{H(W(x))}=x$，$\frac{F(x)}{H(W(x))}=1$，即 $F(x)=H(W(x))$。
设 $G(x)=\frac{1}{1-H(x)}\frac{1}{1-xy}$，那么我们求的就是 $[x^{n+1}]G(W(x))$。
设 $P(x)$ 为 $W(x)$ 的复合逆，根据扩展拉格朗日反演得，$[x^{n+1}]G(W(x))=\frac{1}{n+1}[x^n]G’(x)(\frac{x}{P(x)})^{n+1}$。
$P(x)=\ln(x+1)$，$G’(x)=\frac{y+H’(x)-yH(x)-xyH’(x)}{(1-H(x))^2(1-xy)^2}$，$W(x)=e^x-1$，$H(x)=\frac{x}{\ln(x+1)}$。
所以我们要求的实际上是 $\frac{1}{n+1}[x^ny^{n-j+1}]\frac{y(1-H(x))+H’(x)(1-xy)}{(1-H(x))^2(1-xy)^2}H(x)^{n+1}$。
\[\frac{1}{n+1}[x^ny^{n-j+1}]\frac{y(1-H(x))+H'(x)(1-xy)}{(1-H(x))^2(1-xy)^2}H(x)^{n+1}\\ =\frac{1}{n+1}[x^ny^{n-j+1}](\frac{1}{(1-H(x))}\sum\limits_{i=0}^\infin (i+1)x^iy^{i+1}+\frac{H'(x)}{(1-H(x))^2}\sum\limits_{i=0}^\infin x^iy^i)H(x)^{n+1}\\ =\frac{1}{n+1}[x^n](\frac{1}{(1-H(x))}(n-j+1)x^{n-j}+\frac{H'(x)}{(1-H(x))^2}x^{n-j+1})H(x)^{n+1}\\ =\frac{1}{n+1}((n-j+1)[x^j]\frac{H(x)^{n+1}}{(1-H(x))}+[x^{j-1}]\frac{H'(x)H(x)^{n+1}}{(1-H(x))^2})\\ =\frac{1}{n+1}((n-j+1)[x^{j+1}]\frac{H(x)^{n+1}}{x^{-1}(1-H(x))}+[x^{j+1}]\frac{H'(x)H(x)^{n+1}}{x^{-2}(1-H(x))^2})\\\]就可以直接算了。
#include<cstdio>
#include<algorithm>
int const mod=998244353;
int pow(int x,int y){
int res=1;
while(y){
if(y&1)res=1ll*res*x%mod;
x=1ll*x*x%mod;
y>>=1;
}
return res;
}
void dft(int *a,int lim){
for(int mid=lim>>1;mid;mid>>=1)
for(int rt=pow(3,(mod-1)/(mid<<1)),r=mid<<1,j=0;j<lim;j+=r)
for(int p=1,k=0;k<mid;k++,p=1ll*p*rt%mod){
int x=a[j+k],y=a[j+mid+k];
a[j+k]=(x+y)%mod,a[j+mid+k]=1ll*(x-y+mod)*p%mod;
}
}
void idft(int *a,int lim){
for(int mid=1;mid<lim;mid<<=1)
for(int rt=pow(332748118,(mod-1)/(mid<<1)),r=mid<<1,j=0;j<lim;j+=r)
for(int p=1,k=0;k<mid;k++,p=1ll*p*rt%mod){
int x=a[j+k],y=1ll*p*a[j+mid+k]%mod;
a[j+k]=(x+y)%mod,a[j+mid+k]=(x-y+mod)%mod;
}
for(int i=0,p=pow(lim,mod-2);i<lim;i++)a[i]=1ll*a[i]*p%mod;
}
void inv(int const *a,int *ans,int lim){
static int b[400010];
for(int i=0;i<lim<<1;i++)ans[i]=b[i]=0;
ans[0]=pow(a[0],mod-2);
for(int n=2;n<=lim;n<<=1){
for(int i=0;i<n;i++)b[i]=a[i];
dft(ans,n<<1),dft(b,n<<1);
for(int i=0;i<n<<1;i++)ans[i]=ans[i]*(2-1ll*ans[i]*b[i]%mod+mod)%mod,b[i]=0;
idft(ans,n<<1);
for(int i=n;i<n<<1;i++)ans[i]=0;
}
}
void der(int const *a,int *ans,int lim){
for(int i=1;i<lim;i++)ans[i-1]=1ll*i*a[i]%mod;
ans[lim-1]=0;
}
void inte(int const *a,int *ans,int lim){
for(int i=1;i<lim;i++)ans[i]=1ll*pow(i,mod-2)*a[i-1]%mod;
ans[0]=0;
}
void ln(int const *a,int *ans,int lim){
static int b[400010];
for(int i=0;i<lim<<1;i++)b[i]=0;
der(a,b,lim);
inv(a,ans,lim);
dft(b,lim<<1),dft(ans,lim<<1);
for(int i=0;i<lim<<1;i++)b[i]=1ll*b[i]*ans[i]%mod,ans[i]=0;
idft(b,lim<<1);
inte(b,ans,lim);
}
void exp(int const *a,int *ans,int lim){
static int b[400010];
for(int i=0;i<lim<<1;i++)b[i]=ans[i]=0;
ans[0]=1;
for(int n=2;n<=lim;n<<=1){
ln(ans,b,n);
b[0]=(a[0]-b[0]+1+mod)%mod;
for(int i=1;i<n;i++)b[i]=(a[i]-b[i]+mod)%mod;
dft(b,n<<1),dft(ans,n<<1);
for(int i=0;i<n<<1;i++)ans[i]=1ll*ans[i]*b[i]%mod;
idft(ans,n<<1);
for(int i=n;i<n<<1;i++)ans[i]=0;
}
}
void pow(int const *a,int k,int *ans,int lim){
static int b[400010];
ln(a,b,lim);
for(int i=0;i<lim;i++)b[i]=1ll*b[i]*k%mod;
exp(b,ans,lim);
}
int fac[400010],iv[400010],f[400010],g[400010],h[400010],n,lim=1<<17,w[400010],q[400010],e[400010];
int main(){
scanf("%d",&n);
fac[0]=1;
for(int i=1;i<=lim+10;i++)fac[i]=1ll*fac[i-1]*i%mod;
iv[lim+10]=pow(fac[lim+10],mod-2);
for(int i=lim+10;i;i--)iv[i-1]=1ll*iv[i]*i%mod;
for(int i=0;i<lim;i++)f[i]=(mod-iv[i+2])%mod;
inv(f,g,lim);
for(int i=0;i<lim;i++)g[i]=g[i+1],f[i]=(i&1?mod-1ll:1ll)*pow(i+1,mod-2)%mod;
inv(f,h,lim);
pow(h,n+1,w,lim);
der(h,q,lim);
for(int i=0;i<lim;i++)f[i]=0;
dft(w,lim<<1),dft(q,lim<<1);
for(int i=0;i<lim<<1;i++)q[i]=1ll*q[i]*w[i]%mod;
idft(q,lim<<1);
for(int i=lim;i<lim<<1;i++)q[i]=0;
for(int i=0;i<lim;i++) h[i]=(mod-h[i+1])%mod;
inv(h,e,lim);
dft(q,lim<<1),dft(e,lim<<1);
for(int i=0;i<lim<<1;i++)q[i]=1ll*q[i]*e[i]%mod,w[i]=1ll*w[i]*e[i]%mod;
idft(q,lim<<1),idft(w,lim<<1);
for(int i=lim;i<lim<<1;i++)q[i]=0;
dft(q,lim<<1);
for(int i=0;i<lim<<1;i++)q[i]=1ll*q[i]*e[i]%mod;
idft(q,lim<<1);
for(int i=1;i<=n;i++)g[i]=(g[i]-1ll*pow(n+1,mod-2)*(1ll*(n-i+1)*w[i+1]%mod+q[i+1])%mod+mod)%mod;
g[0]=n;
for(int i=n+1;i<lim;i++)g[i]=0;
for(int i=0;i<=n;i++)g[i]=1ll*fac[i]*g[i]%mod;
for(int i=0;i<=n;i++)f[i]=(i&1?mod-1ll:1ll)*iv[i]%mod;
for(int i=0;i<=n;i++)if(i<n-i)std::swap(g[i],g[n-i]);
dft(f,lim<<1),dft(g,lim<<1);
for(int i=0;i<lim<<1;i++)f[i]=1ll*f[i]*g[i]%mod;
idft(f,lim<<1);
for(int i=0;i<n;i++)printf("%lld ",1ll*fac[n]*iv[i]%mod*f[n-i]%mod);
return 0;
}