2021 ICPC上海 B. Strange Permutations [NTT+启发式合并]

传送门:CF

[前题提要]:无

---

发现需要求出满足对于每一个$i,Q_{i+1}!=P_{Q_i}$的排列个数,显然的,会想到应该容斥,赛时因为没啥时间了,所以想了个容斥感觉不是很能做就换题了,事实上本题容斥是可做的.

对于本题的容斥,显然我们应该枚举当前满足$Q_{i+1}!=P_{Q_i}$的$i$的个数.你会发现事实上我们并不需要知道对于枚举的$i$来说,具体是那几个条件,只需要求出个数即可.所以假设$f(i)$表示存在$i$个不合法的下标的排列的个数.那么我们最终的答案就是${\sum }_{i=0}^{n-1}(-1{)}^{i}f(i)$.

那么现在我们的问题是如何求出$f(i)$.为了考虑方便起见,我们不妨将原题使用图意义来描述,考虑建出一个完全图,那么对于一个排列来说,等价于完全图中的一个哈密顿路径.考虑设限制边为$<a,b>$,不难发现限制相当于$p_a=b$,也就是限制边即为$<a,p_a>$.那么考虑我们现在选出了$i$条限制边(注意此时我们需要满足选出的限制边不能成环,不然这种情况是不合法的,因为哈密顿路径不存在环).那么对于$i$条限制边来说,我们设其形成了$x$个联通块,那么此时相当于影响了$i+x$个点,那么此时选边的一个组合数贡献为$(n-x-i+x)!=(n-i)!$,我们发现其只和$i$有关.所以我们现在只需要求出选出限制边的的方案数即可.注意此时并不是直接$C_n^i$,因为我们选出的限制边得需要保证不能成环!!!

那么我们此时观察限制边,我们会发现所有的边都是$<i,p_i>$,这是一个经典的连边,这意味着限制边形成的图应该是一个基环树森林.并且特殊的,我们似乎还忽略了一个条件,对于所有的$i,p_i$都是排列,这意味着什么,这说明事实上我们所有点的入度和出度都是1,那么此时我们的基环树森林就是一个个环.到这里,本题就变得好做起来了.对于每一个环,假设其的大小为$Size$,那么我们显然可以挑选$[0,Size-1]$条边,然后我们需要的是选边和为$i$的方案数.此时是一个经典的背包问题,但是我们的环的个数是$O(n)$级别的,我们的环的大小总和是$O(n)$级别的,我们直接背包复杂度是$O(n^2)$级别的,所以我们考虑使用卷积(事实上对于本题,使用卷积依旧需要启发式合并,然而我们的背包事实上也是可以使用启发式合并优化的,感觉强行使用dp也不是不行,但是细节有点麻烦,不如直接上NTT,本题到这里实在是太卷积了).不难发现对于每一个环的边的贡献多项式是${\sum }_{i=0}^{Size-1}{C}_{Size}^i{x}^i$,所以考虑直接使用NTT合并.当然因为NTT合并两个多项式的复杂度是$Sizelog(Size)$的,直接合并复杂度依旧是$n^2$级别.所以需要用启发式合并,也就是每次选择两个最大系数最小的两个多项式进行合并.最后复杂度是$nlo{g}^{2}n$的.

PS:简单口胡证明一下启发式的复杂度,我们考虑证明对于一个多项式,它乘的次数是log次的.这样最后的复杂度就是$nlo{g}^{2}n$了.考虑如果想要一个多项式乘的次数最多,我们每次都选择尽量的选择它,也就是它每次都是较小的两个之一,不难发现其贡献的最大项的次数是指数级别增长的,因为每过几次,另外一个多项式的最大系数必然比它前几次的要大,所以这里贡献了一个倍数.显然的我们会感觉这个证明事实上是不太充分的(逃,但是由于博主也不是很会证明,所以就当贡献一个证明的思考方向好了

那么至此,本题结束.

---

下面是具体的代码部分:

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define root 1,n,1
#define ls (rt<<1)
#define rs (rt<<1|1)
#define lson l,mid,rt<<1
#define rson mid+1,r,rt<<1|1
#define pd __gnu_pbds
inline ll read() {
	ll x=0,w=1;char ch=getchar();
	for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') w=-1;
	for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
	return x*w;
}
inline void print(__int128 x){
	if(x<0) {putchar('-');x=-x;}
	if(x>9) print(x/10);
	putchar(x%10+'0');
}
#define maxn 1000010
#define int long long
const int mod=998244353;
const double eps=1e-8;
#define	int_INF 0x3f3f3f3f
#define ll_INF 0x3f3f3f3f3f3f3f3f
int f[maxn],g[maxn];int rev[maxn];
int qpow(int a,int b) {
	int ans=1;
	while(b) {
		if(b&1) ans=ans*a%mod;
		b>>=1;
		a=a*a%mod;
	}
	return ans;
}
int fac[maxn];int in_fac[maxn];
void init(int limit=1e5) {
	fac[0]=1;
	for(int i=1;i<=limit;i++) {
		fac[i]=fac[i-1]*i%mod;
	}
	in_fac[limit]=qpow(fac[limit],mod-2);
	for(int i=limit-1;i>=0;i--) {
		in_fac[i]=in_fac[i+1]*(i+1)%mod;
	}
}
int C(int a,int b) {
	return fac[a]*in_fac[b]%mod*in_fac[a-b]%mod;
}
void NTT(int *a,int n,int inv) {
	for(int i=0;i<=n;i++) {
		if(i<rev[i]) swap(a[i],a[rev[i]]);
	}
	for(int len=1;len<=(n>>1);len<<=1) {
		int gn=qpow(inv==1?3:qpow(3,mod-2),(mod-1)/(len<<1));
		for(int i=0;i<=n;i+=(len<<1)) {
			int g0=1;
			for(int j=0;j<=len-1;j++) {
				int x=a[i+j],y=g0*a[i+j+len]%mod;
				a[i+j]=(x+y)%mod;a[i+j+len]=(x-y+mod)%mod;
				g0=g0*gn%mod;
			}
		}
	}
	return ;
}
int p[maxn];vector<int>edge[maxn];
int Circle[maxn];int tot=0;int vis[maxn];int cnt=0;
void dfs(int u,int pre_u) {
	cnt++;vis[u]=1;
	for(auto v:edge[u]) {
		if(v==pre_u||vis[v]) continue;
		dfs(v,u);
	}
}
struct heapnode{
	int Size;vector<int>F;
	bool operator < (const heapnode &rhs) const {
		return Size>rhs.Size;
	}
};
int temp1[maxn],temp2[maxn];int N,M;int n;int ans=0;
void solve() {
	priority_queue<heapnode>q;
	for(int i=1;i<=tot;i++) {
		vector<int>temp;
		for(int j=0;j<Circle[i];j++) {
			int num=C(Circle[i],j);
			temp.push_back(num);
		}
		q.push({(int)temp.size(),temp});
	}
	while(1) {
		if(q.size()==1) break;
		auto f1=q.top();q.pop();auto f2=q.top();q.pop();
		N=0;M=0;
		for(auto i:f1.F) {
			temp1[N++]=i;
		}
		for(auto i:f2.F) {
			temp2[M++]=i;
		}
		N--;M--;
		int limit=1,len=0;
		while(limit<=N+M) len++,limit<<=1;
		for(int i=1;i<=limit;i++) {
			rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
		}
		NTT(temp1,limit,1);NTT(temp2,limit,1);
		for(int i=0;i<=limit;i++) {
			temp1[i]=temp1[i]*temp2[i]%mod;
		}
		NTT(temp1,limit,-1);
		int inv=qpow(limit,mod-2);
		vector<int>temp;
		for(int i=0;i<=N+M;i++) {
			temp.push_back(temp1[i]*inv%mod);
		}
		q.push({(int)temp.size(),temp});
		for(int i=0;i<=limit;i++) {
			temp1[i]=temp2[i]=0;
		}
	}
	auto f=q.top();
	for(int i=0;i<f.F.size();i++) {
		int state=(i&1)?-1:1;
		ans+=state*fac[n-i]%mod*f.F[i]%mod;ans%=mod;
	}
	ans=(ans%mod+mod)%mod;
}
signed main() {
	init();
	n=read();
	for(int i=1;i<=n;i++) {
		p[i]=read();
		edge[i].push_back(p[i]);
		edge[p[i]].push_back(i);
	}
	for(int i=1;i<=n;i++) {
		if(!vis[i]) {
			cnt=0;
			dfs(i,0);
			Circle[++tot]=cnt;
		}
	}
	solve();
	cout<<ans<<endl;
	return 0;
}