题解 LOJ3284 「USACO 2020 US Open Platinum」Exercise

对于一个排列，它的操作次数是它所有循环圈长度的(operatorname{lcm})。

因为答案是求乘积。我们可以分别计算每个质数的贡献。

对于一个质数(p)，它对答案的贡献是每个排列中，所有循环圈长度（分解质因数后）(p)的次数的最大值之和。即：

[ans=prod_{p}p^{sum_{ ext{一个排列}}max_{ ext{一个循环圈}i}{len_i ext{里}p ext{的次数}}} ]

考虑求指数部分，对(m-1)取模。即：

[sum_{ ext{一个排列}}max_{ ext{一个循环圈}i}{len_i ext{里}p ext{的次数}} ]

考虑min-max容斥。则上式可以转化为：

[sum_{ ext{一个排列}}sum_{ ext{循环圈的一个子集}s}(-1)^{|s|+1}min_{iin s}{len_i ext{里}p ext{的次数}} ]

最小值，可以转化为枚举一个值(x)，若所有数都(geq x)，则对最小值(+1)。于是上式可以写成：

[sum_{ ext{一个排列}}sum_{ ext{循环圈的一个子集}s}(-1)^{|s|+1}sum_{x=1}^{n}prod_{iin s}[len_i ext{里}p ext{的次数}geq x] ]

我们枚举质数(p)，再枚举(x)。问题转化为，求所有排列的，所有循环圈的子集中，有多少个子集使得子集内所有循环圈的大小都是(p^x)的倍数。每个子集要带上容斥系数((-1)^{|s|+1})。

考虑DP。设子集内所有循环圈大小之和是(icdot p^x)的方案数为(dp[i])。转移时，类似于带标号无向连通图计数的方法，我们钦定一个点为一号点。枚举一号点所在的循环圈的大小为(jcdot p^x)。则可以这么转移：

[dp[i]=sum_{j=1}^{i}(-1)cdot dp[i-j]cdot {icdot p^x-1choose jcdot p^x-1}(jcdot p^x-1)! ]

DP部分的时间复杂度是(O(frac{n^2}{p^{2x}}))的。在外层枚举(p)和(x)后，总复杂度是(n^2sum_{p}sum_{x}frac{1}{p^{2x}}leq n^2sum_{i=1}^{n}frac{1}{i^2}leqfrac{n^2pi^2}{6})，也就是(O(n^2))。

参考代码：

在LOJ查看

//LOJ3284
#include <bits/stdc++.h>
using namespace std;

#define pb push_back
#define mk make_pair
#define lob lower_bound
#define upb upper_bound
#define fi first
#define se second
#define SZ(x) ((int)(x).size())

typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;

/*  ------  by:duyi  ------  */ // myt天下第一
const int MAXN=7500;
int n,m,fac[MAXN+5],comb[MAXN+5][MAXN+5],p[MAXN+5],cnt_p,flag[MAXN+5];
bool v[MAXN+5];
int pow_mod(int x,int i,int MOD){
	int y=1;
	while(i){
		if(i&1)y=(ll)y*x%MOD;
		x=(ll)x*x%MOD;
		i>>=1;
	}
	return y;
}
int calc(int s,int MOD){
	static int dp[MAXN+5];
	int lim=n/s;
	dp[0]=MOD-1;
	for(int i=1;i<=lim;++i){
		dp[i]=0;
		for(int j=1;j<=i;++j)if(dp[i-j]){//钦定一个一号点,一号点所在的环长为j*s
			dp[i]=(dp[i]-(ll)dp[i-j]*comb[i*s-1][j*s-1]%MOD*fac[j*s-1]%MOD+MOD)%MOD;
		}
	}
	int res=0;
	for(int i=1;i<=lim;++i)res=(res+(ll)dp[i]*comb[n][i*s]%MOD*fac[n-i*s]%MOD)%MOD;
	return res;
}
int main() {
	cin>>n>>m;
	fac[0]=1;
	for(int i=1;i<=n;++i)fac[i]=(ll)fac[i-1]*i%(m-1);//阶乘
	comb[0][0]=1;
	for(int i=1;i<=n;++i){
		comb[i][0]=1;
		for(int j=1;j<=i;++j){
			int& x=comb[i][j];
			x=comb[i-1][j]+comb[i-1][j-1];
			x=(x<(m-1)?x:x-(m-1));
		}
	}//组合数
	for(int i=2;i<=n;++i){
		if(!v[i])p[++cnt_p]=i;
		for(int j=1;j<=cnt_p&&p[j]*i<=n;++j){
			v[i*p[j]]=1;
			if(i%p[j]==0)break;
		}
	}//质数
	for(int i=1;i<=cnt_p;++i){
		for(int j=p[i];j<=n;j*=p[i]){
			flag[j]=p[i];
		}
	}//质数的幂
	
	int ans=1;
	for(int i=1;i<=n;++i)if(flag[i]){
		ans=(ll)ans*pow_mod(flag[i],calc(i,m-1),m)%m;
	}
	cout<<ans<<endl;
	return 0;
}