题目描述
定义一个正整数序列\(\texttt{p}\)，称其是合法的当且仅当对于所有在\(\texttt{p}\)中出现过且\(> 1\)的正整数\(k\)，存在\(i<j\)，满足\(p_i=k-1,p_j=k\)。
定义\(f(k)\)为\(k\)在所有长度为\(n\)的好序列中的出现次数之和。
\(\forall k\in[1,n]\)求\(f(k)\)。
Easy Version(\(n\leq 5000\))
可证明每个合法序列和一个排列一一对应：
- 1.考虑怎么把排列对应成一个合法序列：
考虑把排列划分成若干极长的连续下降段，若第\(i\)段里有数\(x\)，就让\(p_x=i\)。
这样构造的\(p\)一定是合法的，否则就说明存在一个\(k\)，使得所有\(k\)的出现位置都在\(k-1\)的出现位置前面，那么这两个连续段应该会变成一个，就矛盾了。
- 2.一个合法序列对应一个排列
这是显然的，把上述过程逆过来就好了。
因此这是一个双射。
那么对于一个位置\(i\)，它对\(f(k)\)的贡献就是它被划分在第\(k\)段的排列数。
对于位置\(i\)，若\(p_i<p_{i+1}\)，我们称之为一个上升。
则\(i\)被划分在第\(k\)段的排列数等价与前\(i\)个位置恰好有\(k-1\)个上升。
设长度为\(n\)且有\(m\)个上升的排列个数为\(dp_{n,m}\)，则
\[f(k)=\sum_{i=1}^ndp_{i,k-1}C_n^i(n-i)! \]\(dp_{n,m}\)的递推是好求的，考虑加入最小的数：
-
加入后上升个数不变，有\(m+1\)种选择，因为可以填到末尾
-
加入后上升个数加一，有\(n-m\)种选择，因为可以填到开头
总复杂度\(O(n^2)\).
#include<bits/stdc++.h>
using namespace std;
const int N = 5050;
int dp[N][N];
int n;
const int mod = 998244353;
inline int myplus(int a,int b){return (a+b>=mod?a+b-mod:a+b);}
int ans[N],C[N][N],fac[N];
int main()
{
cin>>n;
dp[1][0]=1;
for(int i=1;i<=n;i++)
for(int j=0;j<=i-1;j++)
{
dp[i+1][j]=myplus(dp[i+1][j],1ll*dp[i][j]*(j+1)%mod);
dp[i+1][j+1]=myplus(dp[i+1][j+1],1ll*dp[i][j]*(i-j)%mod);
}
C[0][0]=1;fac[0]=1;
for(int i=1;i<=n;i++)
{
fac[i]=1ll*fac[i-1]*i%mod;
C[i][0]=1;
for(int j=1;j<=i;j++)
C[i][j]=myplus(C[i-1][j-1],C[i-1][j]);
}
for(int k=1;k<=n;k++)
for(int i=1;i<=n;i++)
ans[k]=myplus(ans[k],1ll*dp[i][k-1]*C[n][i]%mod*fac[n-i]%mod);
for(int i=1;i<=n;i++)printf("%d ",ans[i]);
return 0;
}
Hard Version(\(n\leq 10^5\))
上文中的\(f_{n,m}\)其实就是欧拉数，记为\(E_{n,m}\),其含义为有m个上升的n阶排列数量。
\(dp\)的求法显然不能拓展，考虑更直接的刻画。
二项式反演，记\(D_{n,m}\)表示钦定\(m\)个上升，剩下位置随便的方案数。
则
\[D_{n,m}=\sum_{k=m}^n \binom{k}{m}E_{n,k} \]\[E_{n,m}=\sum_{k=m}^n(-1)^{k-m}\binom{k}{m}D_{n,k} \]考虑\(D_{n,m}\)，把连续的上升看成一个连续段，则每个段的方案数唯一，不同段之间是一个二项卷积，一共有\(n-m\)个连续段，故：
连续段的EGF：\(\lbrace0,1,1,……\rbrace\),即\(e^x-1\)。
总的EGF：\((e^x-1)^{n-m}\)
即：\(D_{n,m}=n^{n-m}\)
\[E_{n,m}=n!\sum_{k=m}^n(-1)^{k-m}\binom{k}{m}[x^n](e^x-1)^{n-k} \]\(Easy Version\)的答案式：
\[ans_k=\sum_{i=1}^nE_{i,k-1}\binom{n}{i}(n-i)! \]方便起见，我们改成\(ans_{k+1}\)，这样右边会整齐一点。
\[\begin{aligned} ans_{k+1} &=\sum_{i=1}^nE_{i,k}\binom{n}{i}(n-i)!\\ &=\sum_{i=1}^n\binom{n}{i}(n-i)i!\sum_{j=k}^i(-1)^{j-k}\binom{j}{k}[x^i](e^x-1)^{i-j}!\\ &=n!\sum_{i=1}^n\sum_{j=k}^i(-1)^{j-k}\binom{j}{k}[x^i](e^x-1)^{i-j}!\\ &=\frac{n!}{k!}\sum_{j=k}^n\frac{j!(-1)^{j-k}}{(j-k)!}\sum_{i=j}^n[x^i](e^x-1)^{i-j}!\\ \end{aligned} \]设\(S_j=\sum_{i=j}^n[x^i](e^x-1)^{i-j}!\)
则
\[ans_{k+1}=\frac{n!}{k!}\sum_{j=k}^n\frac{j!(-1)^{j-k}}{(j-k)!}S_{j} \]是一个差卷积的形式，可以直接算。
问题就是求出\(S_j\)
\[\begin{aligned} S_j&=\sum_{i=j}^n[x^i](e^x-1)^{i-j}!\\ &=\sum_{i=j}^n[x^j](\frac{e^x-1}{x})^{i-j}!\\ &=\sum_{i=0}^{n-j}[x^j](\frac{e^x-1}{x})^i!\\ \end{aligned} \]设\(F(x)=\frac{e^x-1}{x}\)
\[S_j=[x^j]\frac{1-F(x)^{n-j+1}}{1-F(x)} \]\[=[x^j](\frac{1}{1-F(x)}-\frac{F(x)^{n-j+1}}{1-F(x)}) \]左边几乎可以直接算，唯一需要注意的是\(1-F(x)\)的常数项为0不能直接算，但是根据某些知识
\[\frac{1}{F(x)\times x^{-1}}\times x^{-1}=\frac{1}{F} \]因此可以先平移一下，再平移回来。
重点看右边，设为\(R_j\)：
\[R_j=[x^j]\frac{F(x)^{n-j+1}}{1-F(x)}=[x^n]\frac{(xF(x))^{n-j+1}}{x-xF(x)}=[x^n]\frac{(e^x-1)^{n-j+1}}{x-e^x+1} \]设\(G(x)=e^x-1,P(x)=\ln(x+1)\)，则\(P(G(x))=x\)，即\(P,G\)互为复合逆。
\[[x^n]\frac{(e^x-1)^{n-j+1}}{x-e^x+1}=[x^n]\frac{G(x)^{n-j+1}}{x-G(x)} \]设\(H_j(x)=\frac{x^{n-j+1}}{P(x)-x}\)，则\(R_j=[x^n]H_j(G(x))\)
由扩展另类拉格朗日反演得：
\[\begin{aligned} R_j &=[x^n]H_j(G(x))\\ &=[x^n]H^j(x)P'(x)(\frac{x}{P(x)})^{n+1}\\ &=[x^n]\frac{x^{n-j+1}}{P(x)-x}P'(x)(\frac{x}{P(x)})^{n+1}\\ &=[x^{k-1}]\frac{1}{P(x)-x}P'(x)(\frac{x}{P(x)})^{n+1} \end{aligned} \]然后就做完了。
注意求逆的时候要先平移一下。
复杂度\(O(n\log n)\)
#include<bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
#define pb push_back
typedef std::vector<int> Poly;
typedef std::pair<int,int> pii;
typedef std::pair<Poly,Poly> ppp;
Poly Pre(Poly A,int L){A.resize(L);return A;}
void Prt(Poly A){for(int x:A)printf("%d ",x);puts("");}
Poly Rev(Poly A){return std::reverse(A.begin(),A.end()),A;}
Poly Shift(Poly A,int L=1){for(int i=L;i<(int)A.size();++i)A[i-L] = A[i];for(int i=max((int)A.size()-L,0);i<(int)A.size();++i)A[i] = 0;return A;}
const int mdig = 18,mod = 998244353,g = 3,maxn = (1<<mdig);
int Pow(int a,int x){int ans = 1,bas = a;while(x){if(x&1)ans = 1ll*ans*bas%mod;bas = 1ll*bas*bas%mod,x >>= 1;}return ans;}
#define getinv(x) Pow(x,mod-2)
const int ivg = Pow(g,mod-2),I = Pow(g,(mod-1)/4);
int gs[2][mdig][maxn>>1];
void NTT(std::vector<int>&A,int N,int op){
static int rev[maxn],flag;static ull tmp[maxn];
if(!flag){flag = 1;for(int op=0;op<2;++op)for(int i=0;i<mdig;++i){int l=(1<<i),w=Pow(op?ivg:g,(mod-1)/(2*l));for(int p=0,c=1;p<l;++p,c=1ll*c*w%mod)gs[op][i][p] = c;}}
for(int i=0;i<N;++i)rev[i] = (rev[i>>1]>>1)|((i&1)*N/2);
for(int i=0;i<N;tmp[i]=A[i],++i)if(i<rev[i])std::swap(A[i],A[rev[i]]);
for(int l=1,d=0;l<N;l*=2,++d){for(int r=0;r<N;r+=(2*l)){int *GS=gs[op][d];ull *tp1=&tmp[l+r],*tp2=&tmp[r];for(int p=0,t;p<l;++p,++GS,++tp1,++tp2)t = 1ll*(*GS)*(*tp1)%mod,(*tp1) = (*tp2)+mod-t,(*tp2) += t;}if(d == 17)for(int i=0;i<N;++i)tmp[i] %= mod;}
for(int i=0,inv=Pow(N,mod-2);i<N;++i)A[i] = tmp[i]*(op?inv:1)%mod;}
Poly Val(int x){Poly ret;ret.pb(x);return ret;}
Poly Mul(Poly A,int x){x = (x%mod+mod)%mod;for(int &v:A)v = 1ll*v*x%mod;return A;}
Poly Der(Poly A){for(int i=0;i<A.size()-1;++i)A[i] = 1ll*A[i+1]*(i+1)%mod;*A.rbegin() = 0;return A;}
Poly Int(Poly A){static int ifac[maxn]={0,1},now=1;if(now < A.size())for(int i=now+1;i<=A.size();++i)ifac[i] = mod-1ll*(mod/i)*ifac[mod%i]%mod;A.pb(0);for(int i=A.size()-1;i;--i)A[i] = 1ll*A[i-1]*ifac[i]%mod;A[0] = 0;return A;}
Poly Pls(Poly A,Poly B){A.resize(max(A.size(),B.size()));for(int i=0;i<B.size();++i)A[i] = (A[i]+B[i])%mod;return A;}
Poly Mns(Poly A,Poly B){A.resize(max(A.size(),B.size()));for(int i=0;i<B.size();++i)A[i] = (A[i]-B[i]+mod)%mod;return A;}
Poly operator +(Poly A,Poly B){return Pls(A,B);}
Poly operator -(Poly A,Poly B){return Mns(A,B);}
Poly operator *(Poly A,int x){return Mul(A,x);}
Poly operator *(int x,Poly A){return Mul(A,x);}
Poly Times(Poly A,Poly B,int Len=-1){int L = 1;Len = ~Len?Len:A.size()+B.size()-1;while(L<A.size()+B.size()-1)L *= 2;A.resize(L),B.resize(L);NTT(A,L,0),NTT(B,L,0);for(int i=0;i<L;++i)A[i] = 1ll*A[i]*B[i]%mod;NTT(A,L,1); A.resize(Len);return A; }
Poly Inv(Poly A,int Len=-1){Len = ~Len?Len:A.size();Poly now;now.pb(getinv(A[0]));int dig=1;while(dig<Len){dig *= 4;Poly t1 = Pre(now,dig),t2 = Pre(Pre(A,dig/2),dig);NTT(t1,dig,0),NTT(t2,dig,0);for(int i=0;i<dig;++i)t1[i] = 1ll*t1[i]*t1[i]%mod*t2[i]%mod;NTT(t1,dig,1),now = Mns(Mul(now,2),Pre(t1,dig/=2));}now.resize(Len);return now;}
Poly Ln(Poly A,int Len=-1){Len = ~Len?Len:A.size();Poly ret = Int(Times(Der(A),Inv(A,Len),Len));ret.resize(Len);return ret;}
Poly Exp(Poly A,int Len=-1){Len = ~Len?Len:A.size();Poly now;now.pb(1);int dig=1;while(dig<Len)dig *= 2,now = Times(now,Mns(Pls(Pre(A,dig),Val(1)),Ln(now,dig)),dig);now.resize(Len);return now;}
Poly Pow(Poly A,int K){return Exp(Mul(Ln(A),K));}
int n,fac[maxn],ifac[maxn],inv[maxn],N;
Poly F,P,G,H;
int main()
{
scanf("%d",&n),N = n+3;
fac[0] = ifac[0] = 1;for(int i=1;i<=N;++i)fac[i] = 1ll*fac[i-1]*i%mod;
ifac[N] = getinv(fac[N]);for(int i=N;i>=2;--i)ifac[i-1] = 1ll*ifac[i]*i%mod;
inv[1] = 1;for(int i=2;i<=N;++i)inv[i] = mod-1ll*(mod/i)*inv[mod%i]%mod;
F.resize(N),P.resize(N);for(int i=1;i<N;++i)F[i] = ifac[i],P[i] = 1ll*(i&1?1:mod-1)*inv[i]%mod;
H = Times(Der(P),Pow(Inv(Shift(P,1)),n+1),N);P[1] = (P[1]-1+mod)%mod;H = Times(Inv(Shift(P,2)),H,N);
for(int i=0;i<n;++i)H[i] = H[i+1];F = Shift(F,1);G = Times(F,Inv(Shift(Val(1)-F,1)),N);
for(int i=0;i<n;++i)H[i] = (G[i+1]+mod-H[i])%mod;G.resize(n),H.resize(n);for(int i=0;i<n;++i)G[i] = 1ll*H[i]*fac[i]%mod*(i&1?mod-1:1)%mod;for(int i=0;i<n;++i)H[i] = ifac[i];
G = Times(Rev(G),H,n);for(int i=0;i<n;++i)printf("%lld ",1ll*(i&1?mod-1:1)*ifac[i]%mod*fac[n]%mod*G[n-i-1]%mod);
return 0;
}