对于一个排列，它的操作次数是它所有循环圈长度的(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))。
参考代码：
//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;
}