[loj3284]Exercise

对于一个排列$p_{i}$，假设循环长度依次为$x_{1},x_{2},...,x_{m}$，那么所需步数即${ m lcm}_{i=1}^{m}x_{i}$

由于是乘积，因此可以枚举素数$p$，并统计其的次数（注意这是对$varphi(M)=M-1$取模）

类似于$E(X)=sum_{ige 1}P(Xge i)$，定义$f(k)$表示满足$exists kmid x_{i}$的排列数，次数即$sum_{alphage 1}f(p^{alpha})$

关于$f(k)$，考虑总排列数为$n!$，并去掉$forall 1le ile m,k otmid x_{i}$的排列数即为所求

令$dp_{i}$表示$n=i$且要求$forall 1le jle m,k otmid x_{j}$时的方案数，同样考虑总排列数为$i!$，并去掉$exists kmid x_{j}$的排列数即为所求，后者去枚举$sum_{kmid x_{j}}x_{j}$，显然剩下的部分即为$dp_{i-sum_{kmid x_{j}}x_{j}}$

由此，转移即
$$
dp_{i}=i!-sum_{1le jle lfloorfrac{i}{k} floor}{ichoose jk}g_{j}dp_{i-jk}
$$
其中$g_{i}$表示$n=ik$且要求$forall 1le jle m,kmid x_{j}$时的方案数，枚举$n$​所在循环长度，转移即
$$
g_{i}=sum_{1le jle i}{ik-1choose jk-1}(jk-1)!g_{i-j}
$$
显然计算$g_{i}$的复杂度为$o(frac{n^{2}}{k^{2}})$，同时由于只需要求$dp_{n}$，根据转移只关心于$iequiv n(mod k)$的$dp_{i}$，因此求$dp_{i}$的复杂度同样为$o(frac{n^{2}}{k^{2}})$，也即求$f(k)$的复杂度为$o(frac{n^{2}}{k^{2}})$

而$k$显然$in [1,n]$且不重复，因此复杂度为$o(sum_{i=1}^{n}frac{n^{2}}{i^{2}})$，提出$n^{2}$后该式收敛于$frac{pi^{2}}{6}$

综上，总复杂度即为$o(n^{2})$，可以通过

#include<bits/stdc++.h>
using namespace std;
#define N 7505
#define ll long long
int n,mod,ans,p[N],vis[N],fac[N],C[N][N],g[N],f[N];
int qpow(int n,int m){
    int s=n,ans=1;
    while (m){
        if (m&1)ans=(ll)ans*s%mod;
        s=(ll)s*s%mod;
        m>>=1;
    }
    return ans;
}
int calc(int k){
    memset(g,0,sizeof(g));
    memset(f,0,sizeof(f));
    g[0]=f[0]=1;
    for(int i=1;i<=n/k;i++)
        for(int j=1;j<=i;j++)g[i]=(g[i]+(ll)C[i*k-1][j*k-1]*fac[j*k-1]%(mod-1)*g[i-j])%(mod-1);
    for(int i=0;i<=n;i++)
        if (i%k==n%k){
            f[i]=fac[i];
            for(int j=1;j<=i/k;j++)f[i]=(f[i]-(ll)C[i][j*k]*g[j]%(mod-1)*f[i-j*k]%(mod-1)+mod-1)%(mod-1);
        }
    return (fac[n]-f[n]+mod-1)%(mod-1);
}
int main(){
    for(int i=2;i<N;i++){
        if (!vis[i])p[++p[0]]=i;
        for(int j=1;(j<=p[0])&&(i*p[j]<N);j++){
            vis[i*p[j]]=1;
            if (i%p[j]==0)break;
        }
    }
    scanf("%d%d",&n,&mod);
    fac[0]=1;
    for(int i=1;i<N;i++)fac[i]=(ll)i*fac[i-1]%(mod-1);
    for(int i=0;i<N;i++){
        C[i][0]=C[i][i]=1;
        for(int j=1;j<i;j++)C[i][j]=(C[i-1][j-1]+C[i-1][j])%(mod-1);
    }
    ans=1;
    for(int i=1;i<=p[0];i++)
        for(int j=p[i];j<=n;j*=p[i])ans=(ll)ans*qpow(p[i],calc(j))%mod;
    printf("%d
",ans);
    return 0;
}