「USACO 2020 US Open Platinum」Exercise

对于一个排列，答案是每个轮换环长的 (operatorname{lcm})

考虑对每个质数算贡献。

[ans=prod_{p}p^{sum_{egeq 1}sum_{Tsubseteq{1,...,n},sum_{c_iin T}=n}e[max(ord_p(c_i))==e]} ]

只要算 (p) 的指数上的式子 (mod p-1) 的结果。

[egin{aligned} sum_{egeq 1}sum_{Tsubseteq{1,...,n},sum_{c_iin T}=n}e[max(ord_p(c_i))==e]&=sum_{egeq}sum_{Tsubseteq{1,...,n},sum_{c_iin T}=n}[max(ord_p(c_i))>=e]\ &=sum_{egeq 1}sum_{Tsubseteq{1,...,n},sum_{c_iin T}=n}sum_{Ssubseteq T}(-1)^{mid S mid +1}prod_{c_i in S}[p^e|c_i]\ &=sum_{egeq 1}sum_{Ssubseteq {1,...,n}}(-1)^{|S| +1}{n choose {|S|}}(n-|S|)!prod_{c_i in S}[p^e|c_i] end{aligned} ]

枚举(e)，对后面的式子dp，(dp_i) 表示有 (i) 个点的答案。设 (x=p^e)
枚举 (i) 号点所在环的大小。(dp_i=-sum_{j=1}^i{{i-1} choose {j-1}}(j-1)!dp_{i-j}[x mid j])

发现有用的dp值只有 (n/x) 个。

枚举(p,e)，总复杂度是 (O(sum_psum_e(frac{n}{p^e})^2)leq O(sum_{x=1}^{n}(frac{n}{x})^2)=O(n^2))

// Author -- Frame

#include<bits/stdc++.h>

#define lowbit(x) ((x)&(-(x)))
#define Finline __inline__ __attribute__ ((always_inline))
#define DEBUG fprintf(stderr,"Running on Line %d in Function %s
",__LINE__,__FUNCTION__)
#define SZ(x) ((int)x.size())
#define mkpr std::make_pair
#define pb push_back

typedef long long ll;
typedef unsigned int uint;
typedef unsigned long long ull;
typedef std::pair<int,int> pi;
typedef std::pair<ll,ll> pl;

const int inf=0x3f3f3f3f,Inf=0x7fffffff;
const ll INF=0x3f3f3f3f3f3f3f3f;

Finline uint rnd()
{
	static uint seed=416;
	return seed^=seed>>5,seed^=seed<<17,seed^=seed>>13;
}
template <typename _Tp>_Tp gcd(const _Tp &a,const _Tp &b){return (!b)?a:gcd(b,a%b);}
template <typename _Tp>Finline _Tp abs(const _Tp &a){return a>=0?a:-a;}
template <typename _Tp>Finline _Tp max(const _Tp &a,const _Tp &b){return a<b?b:a;}
template <typename _Tp>Finline _Tp min(const _Tp &a,const _Tp &b){return a<b?a:b;}
template <typename _Tp>Finline void chmax(_Tp &a,const _Tp &b){(a<b)&&(a=b);}
template <typename _Tp>Finline void chmin(_Tp &a,const _Tp &b){(b<a)&&(a=b);}
template <typename _Tp>Finline void read(_Tp &x)
{
	char ch(getchar());
	bool f(false);
	while(ch<48||ch>57) f|=ch==45,ch=getchar();
	x=ch&15,ch=getchar();
	while(ch>=48&&ch<=57) x=(((x<<2)+x)<<1)+(ch&15),ch=getchar();
	if(f) x=-x;
}
template <typename _Tp,typename... Args>Finline void read(_Tp &t,Args &...args){read(t);read(args...);}
template <typename _Tp,typename... Args>Finline _Tp min(const _Tp &a,const _Tp &b,const Args &...args){return a<b?min(a,args...):min(b,args...);}
template <typename _Tp,typename... Args>Finline _Tp max(const _Tp &a,const _Tp &b,const Args &...args){return a<b?max(b,args...):max(a,args...);}
Finline int read_str(char *s)
{
	char ch(getchar());
	while(ch==' '||ch=='
'||ch=='
') ch=getchar();
	char *tar=s;
	*tar=ch,ch=getchar();
	while(ch!=' '&&ch!='
'&&ch!='
'&&ch!=EOF) *(++tar)=ch,ch=getchar();
	return tar-s+1;
}

const int N=7505;
int mod;
template<typename _Tp1,typename _Tp2>Finline void add(_Tp1 &a,_Tp2 b){(a+=b)>=mod&&(a-=mod);}
template<typename _Tp1,typename _Tp2>Finline void sub(_Tp1 &a,_Tp2 b){(a-=b)<0&&(a+=mod);}
struct FastMod{
	typedef __uint128_t LLL;
	ull b,m;
	void init(ull b) {this->b=b,m=ull((LLL(1)<<64)/b);}
	Finline ull operator()(ull a){
		ull r=a-(ull)((LLL(m)*a)>>64)*b;
		return r>=b?r-b:r;
	}
}M;
int fac[N],dp[N],C[N][N];
bool pr[N];
int p[N],pos;
void sieve()
{
	for(int i=2;i<N;++i)
	{
		if(!pr[i]) p[++pos]=i;
		for(int j=1;j<=pos&&i*p[j]<N;++j)
		{
			pr[i*p[j]]=true;
			if(!(i%p[j])) break;
		}
	}
}
ll ksm(ll a,ll b,ll mod)
{
	ll res=1;
	while(b)
	{
		if(b&1) res=res*a%mod;
		a=a*a%mod,b>>=1;
	}
	return res;
}
int main()
{
	sieve();
	int n,P;read(n,P);mod=P-1;M.init(mod);
	fac[0]=1;for(int i=1;i<=n;++i) fac[i]=M(1ll*fac[i-1]*i);
	for(int i=0;i<=n;++i) for(int j=C[i][0]=1;j<=i;++j) C[i][j]=C[i-1][j-1],add(C[i][j],C[i-1][j]);
	int ans=1;
	for(int _=1;_<=pos&&p[_]<=n;++_)
	{
		int tmp=0;
		for(int cur=p[_];cur<=n;cur*=p[_])
		{
			int cnt=n/cur;
			memset(dp,0,(cnt+3)<<2);dp[0]=mod-1;
			for(int i=1;i<=cnt;++i) for(int j=1;j<=i;++j) sub(dp[i],M(M(1ll*C[i*cur-1][j*cur-1]*fac[j*cur-1])*dp[i-j]));
			for(int i=1;i<=cnt;++i) add(tmp,M(M(1ll*C[n][i*cur]*fac[n-i*cur])*dp[i]));
		}
		if(tmp) ans=1ull*ans*ksm(p[_],tmp,P)%P;
	}
	printf("%d
",ans);
	return 0;
}