[Comet OJ - Contest #6 F]permutation

Description

给出n，求有多少个长度为n的排列没有长度为2~n-1的连续段
n<=10^5

Solution

生成函数学艺不精.jpg
特判掉n<=2，我们考虑一般情况
首先数排列等价于数析合树，我们考虑根节点的形态
若根为合点，则儿子排列为单调上升/下降，且不存在一个儿子为合点，且这个儿子的儿子也为上升/下降
显然上升=下降，设G(x)表示根为合点且儿子单调上升的析合树数量，H(x)表示排列的OGF，考虑儿子数量我们有$G(x)={\sum }_{i\ge 2}(H(x)-G(x){)}^i$，因为每个儿子都不能和自己状态相同
化出来就是$G(x)=\frac{H(x{)}^2}{1+H(x)}$，G可以用多项式求逆求得
考虑设F(x)为答案，即>=4个点形成的无连续段排列的数量的生成函数，那根节点为析点的生成函数就是$$\sum _{i\ge 4}{f}_{i}H(x{)}^i=F(H(x))$$
那么，一棵析合树，根节点要么是合点，要么是析点，要么是单点，于是我们有$$F(H(x))+2G(x)+x=H(x)$$
可以变成$$F(H(x))=P(x)$$
考虑H(x)的复合逆$H(H^{-1}(x))=x$，我们有$$F(x)=P(H^{-1}(x))$$
我们现在要求$[x^n]F(x)$即$[x^n]P(H^{-1}(x))$，根据扩展拉格朗日反演我们有$$[x^n]P(H^{-1}(x))=\frac{1}{n}[x^{n-1}]P^{\prime }(x)(\frac{x}{H(x)}{)}^n$$
直接多项式ln和exp求幂即可做到O(n log n)
OEISA111111，但没有通项

Code

#include <cstdio>
#include <cstring>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
using namespace std;

typedef long long ll;

const int N=(1<<18)+5,Mo=998244353;

ll pwr(ll x,ll y) {
	ll z=1;
	for(;y;y>>=1,x=x*x%Mo)
		if (y&1) z=z*x%Mo;
	return z;
}

int n,len;
ll W[2][N],inv[N],fac[N],h[N],ih[N],g[N],p[N],hn[N];

void pre(int N) {
	fac[0]=1;fo(i,1,N) fac[i]=fac[i-1]*i%Mo;
	inv[0]=inv[1]=1;fo(i,2,N) inv[i]=(ll)-Mo/i*inv[Mo%i]%Mo;
	for(int i=1;i<N;i<<=1) {
		ll wn=pwr(3,(Mo-1)/(i<<1));
		for(int j=0;j<i;j++) W[1][i+j]=j?W[1][i+j-1]*wn%Mo:1;
		wn=pwr(3,Mo-1-(Mo-1)/(i<<1));
		for(int j=0;j<i;j++) W[0][i+j]=j?W[0][i+j-1]*wn%Mo:1;
	}
}

void DFT(ll *a,int len,int flag) {
	for(int i=0,j=0;i<len;++i) {
        	if (i<j) swap(a[i],a[j]);
        	for(int k=len>>1;(j^=k)<k;k>>=1);
        }
	if (flag==-1) flag=0;
	for(int i=1;i<len;i<<=1)
		for(int j=0;j<len;j+=i<<1)
			for(int k=0;k<i;k++) {
				ll u=a[j+k],v=a[j+k+i]*W[flag][i+k];
				a[j+k]=(u+v)%Mo;a[j+k+i]=(u-v)%Mo;
			} 
	if (flag==0) for(int i=0;i<len;i++) a[i]=a[i]*inv[len]%Mo;
}

void get_inv(ll *a,ll *b,int n) {
	if (n==1) {
		b[0]=pwr(a[0],Mo-2);
		return;
	}
	get_inv(a,b,(n+1)>>1);
	static ll c[N];
	int len=1;for(;len<n<<1;len<<=1);
	fo(i,0,len-1) c[i]=0;
	fo(i,0,n-1) c[i]=a[i];
	fo(i,(n+1)>>1,len-1) b[i]=0;
	DFT(b,len,1);DFT(c,len,1);
	fo(i,0,len-1) b[i]=(2*b[i]-b[i]*b[i]%Mo*c[i])%Mo;
	DFT(b,len,-1);
	fo(i,n,len-1) b[i]=0;
}

void get_ln(ll *a,ll *b,int n) {
	static ll c[N],d[N];
	fo(i,0,n-2) c[i]=a[i+1]*(i+1)%Mo;c[n-1]=0;
	get_inv(a,d,n);
	int len=1;for(;len<n<<1;len<<=1);
	fo(i,n,len-1) c[i]=d[i]=0;
	DFT(c,len,1);DFT(d,len,1);
	fo(i,0,len-1) c[i]=c[i]*d[i]%Mo;
	DFT(c,len,-1);
	b[0]=0;fo(i,1,n-1) b[i]=c[i-1]*pwr(i,Mo-2)%Mo;
}

void get_exp(ll *a,ll *b,int n) {
	if (n==1) {b[0]=1;return;}
	get_exp(a,b,(n+1)>>1);
	static ll c[N];
	int len=1;for(;len<n<<1;len<<=1);
	get_ln(b,c,n);fo(i,n,len-1) c[i]=0;
	fo(i,0,n-1) c[i]=(a[i]-c[i])%Mo;c[0]++;
	DFT(c,len,1);DFT(b,len,1);
	fo(i,0,len-1) b[i]=b[i]*c[i]%Mo;
	DFT(b,len,-1);
	fo(i,n,len-1) b[i]=0;
}

void get_pow(ll *a,ll *b,int n,int m) {
	static ll c[N];
	get_ln(a,c,n);
	fo(i,0,n-1) c[i]=c[i]*m%Mo;
	get_exp(c,b,n);
}

int main() {
	pre(1<<18);
	scanf("%d",&n);
	if (n==1) {puts("1");return 0;}
	if (n==2) {puts("2");return 0;}
	for(len=1;len<=n<<1;len<<=1);
	h[0]=1;fo(i,1,n) h[i]=fac[i];get_inv(h,ih,n+1);h[0]=0;
	DFT(h,len,1);fo(i,0,len-1) g[i]=h[i]*h[i]%Mo;DFT(h,len,-1);
	DFT(g,len,-1);fo(i,n+1,len-1) g[i]=0;
	DFT(g,len,1);DFT(ih,len,1);
	fo(i,0,len-1) g[i]=g[i]*ih[i]%Mo;DFT(g,len,-1);
	fo(i,0,n) p[i]=(h[i]-2*g[i])%Mo;p[1]--;
	fo(i,0,n-1) p[i]=p[i+1]*(i+1)%Mo;fo(i,n,len-1) p[i]=0;
	fo(i,0,n-1) h[i]=h[i+1];h[n]=0;fo(i,0,len-1) ih[i]=0;
	get_inv(h,ih,n);get_pow(ih,hn,n,n);
	DFT(hn,len,1);DFT(p,len,1);
	fo(i,0,len-1) p[i]=p[i]*hn[i]%Mo;
	DFT(p,len,-1);
	ll ans=p[n-1]*pwr(n,Mo-2)%Mo;
	printf("%lld\n",(ans+Mo)%Mo);
	return 0;
}