P5825排列计数【EGF,NTT】

正题

题目链接:https://www.luogu.com.cn/problem/P5825


题目大意

对于每个\(k\),求有多少个长度为\(n\)的排列有\(k\)个位置上升。
\(1\leq n\leq 2\times 10^5\)


解题思路

考虑到同时考虑大于和小于十分麻烦,设\(f_i\)表示钦定\(i\)个上升时的方案

连续的上升段可以视为同一个组,那么整个序列就会被分为\(m=n-k\)段,每个组内都是无序的。

所以可以考虑一下\(\text{EGF}\)来做,因为不能选空段,那么每一段的生成函数就是\(e^x-1\)

也就是\(f_{n-m}=(e^x-1)^m[x^n]\)。二项式定理展开一下

\[f_m=\sum_{i=0}^m\binom{m}{i}(-1)^{m-i}e^{ix} \]

\[=\sum_{i=0}^m\frac{m!}{i!(m-i)!}(-1)^{m-i}\frac{i^n}{n!} \]

\[=\frac{m!}{n!}\sum_{i=0}^m\frac{(-1)^{m-i}}{(m-i)!}\frac{i^n}{i!} \]

\(\text{NTT}\)卷起来就好了。

然后\(g_i\)表示恰好有\(i\)个的话,上二项式反演即可

\[f_i=\sum_{j=0}^i\binom{i}{j}g_j\Rightarrow g_i=\sum_{j=i}(-1)^{j-i}\binom{j}{i}f_j \]

这个也是显然可以卷积快速求得的。

顺带一提的是,这个求得其实就是欧拉数\(\left\langle\begin{matrix} n\\k\end{matrix}\right\rangle\)

联立上面的\(f_i\)\(g_i\)的式子可以得到欧拉数的通式

\[\left\langle\begin{matrix} n\\k\end{matrix}\right\rangle=\sum_{i=0}^{n-k}(-1)^{n-k-i}i^n\binom{n+1}{k+j+1} \]

这个可以一次卷积求得

时间复杂度\(O(n\log n)\)


code

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const ll N=8e5+10,P=998244353;
ll n,m,inv[N],fac[N],f[N],g[N],r[N];
ll power(ll x,ll b){
	ll ans=1;
	while(b){
		if(b&1)ans=ans*x%P;
		x=x*x%P;b>>=1;
	}
	return ans;
}
void NTT(ll *f,ll op){
	for(ll i=0;i<m;i++)
		if(i<r[i])swap(f[i],f[r[i]]);
	for(ll p=2;p<=m;p<<=1){
		ll tmp=power(3,(P-1)/p),len=(p>>1);
		if(op==-1)tmp=power(tmp,P-2);
		for(ll k=0;k<m;k+=p){
			ll buf=1;
			for(ll i=k;i<k+len;i++){
				ll tt=f[i+len]*buf%P;
				f[i+len]=(f[i]-tt+P)%P;
				f[i]=(f[i]+tt)%P;
				buf=buf*tmp%P;
			}
		}
	}
	if(op==-1){
		ll invn=power(m,P-2);
		for(ll i=0;i<m;i++)
			f[i]=f[i]*invn%P;
	}
	return;
}
signed main()
{
	scanf("%lld",&n);
	inv[1]=1;
	for(ll i=2;i<=n;i++)
		inv[i]=P-inv[P%i]*(P/i)%P;
	inv[0]=fac[0]=1;
	for(ll i=1;i<=n;i++)
		inv[i]=inv[i-1]*inv[i]%P,fac[i]=fac[i-1]*i%P;
	for(ll i=0;i<=n;i++)
		f[i]=inv[i]*power(i,n)%P,g[i]=(i&1)?(P-inv[i]):inv[i];
	m=1;while(m<=2*n)m<<=1;
	for(ll i=0;i<m;i++)
		r[i]=(r[i>>1]>>1)|((i&1)?(m>>1):0);
	NTT(f,1);NTT(g,1);
	for(ll i=0;i<m;i++)f[i]=f[i]*g[i]%P;
	NTT(f,-1);
	memset(g,0,sizeof(g));
	for(ll i=0;i<m;i++)
		f[i]=(i<n)?(f[i+1]*fac[i+1]%P):0;
	for(ll i=n;i<m;i++)f[i]=0;
	for(ll i=0;i<n;i++){
		f[i]=f[i]*fac[n-i-1];
		f[i]=(i&1)?(P-f[i]):f[i];
		g[i]=inv[i];
	}
	NTT(f,1);NTT(g,1);
	for(ll i=0;i<m;i++)f[i]=f[i]*g[i]%P;
	NTT(f,-1);
	for(ll i=0;i<n-i-1;i++)swap(f[i],f[n-i-1]);
	for(ll i=0;i<n;i++){
		f[i]=f[i]*inv[i]%P;
		f[i]=((n-i)&1)?f[i]:(P-f[i]);
		printf("%lld ",f[i]%P);
	}
	putchar('0');
	return 0;
}
原文地址:https://www.cnblogs.com/QuantAsk/p/14620442.html