noip模拟 *置换

暴力就是枚举每种情况，算 \(lcm\)。
然后仔细想想（虽然考场没想）就会发现，最后图会形成若干个环
我们要求的其实就是环大小的 \(lcm\)
根据这个性质，可以得到一个 \(60\) 分的 \(dp\)
令 \(dp_{i,j}\) 表示一共选了 \(i\) 个数，\(lcm\) 为 \(j\) 的方案数
我们枚举环的大小 \(s\) 和数量 \(a\)，可以转移

\[dp_{i+s*a,lcm(j,s)}=dp_{i,j}*\prod_{z=1}^{a}\binom{i+s*z}{s}(a!)^{-1}*((s-1)!)^{a} \]

转移的第二项是在求选数的方案，最好自己理解一下
第三项是求单位环连边的方案

正解

\(lcm\) 过大（状态数太大），在写 \(60\) 分的时候还要写 \(vec\)
考虑怎么优化
把一个 \(lcm\) 拆成若干质因子，发现大于 \(\sqrt{n}\) 的质因子出现次数不会大于 \(1\)
那我们把所有小于 \(\sqrt{n}\) 的质因子可以构成的 \(lcm\) 处理出来，发现数量不超过 \(2000\) 个
那这道题就可以做了
处理出每个数被小于 \(\sqrt{n}\) 的质因子处理尽的值（设其为 \(t\) ），把 \(t*t\) 乘进方案就好了
因为 \((ab)^2=a^2b^2\)

Code

#include <bits/stdc++.h>
#define re register
#define db double
#define int long long
#define ll long long
#define pir make_pair
#define memset(a) { for(re int i=1,lim=2*n;i<=lim;i++) for(re int j=1,lim=2*n;j<=lim;j++) a[i][j]=0; }
using namespace std;
const int maxn=5e3+10;
const int INF=1e9+10;
// const int mol=998244353;
inline int qpow(int a,int b,int mol) { int ans=1; while(b) { if(b&1) (ans*=a)%=mol; (a*=a)%=mol; b>>=1; } return ans; }
inline int read(){
    int x=0,w=1;char ch=getchar();
    while(ch<'0'||ch>'9') w=(ch=='-')?-1:1,ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    return x*w;
}

int prim[10]={0,2,3,5,7,11,13,17,19};
int n,mol,a[maxn],g[520][maxn],f[520][maxn],fac[maxn],inv[maxn];
inline int C(int n,int m) { return fac[n]*inv[n-m]%mol*inv[m]%mol; }
inline int lcm(int a,int b) { return a*b/__gcd(a,b); }
int ks[maxn],tot; int pre[maxn],dik[maxn];
inline void dfs(int num,int sum,int x) { 
    if(num>n) { return; }
    if(x>6) { ks[++tot]=sum; return; }
    int tmp=1;
    dfs(num,sum,x+1);
    for(re int i=1;i<=n;i++) {
        tmp*=prim[x];
        if(num+tmp>n) break;
        dfs(num+tmp,sum*tmp,x+1);
    }
}
signed main(void) {
    // freopen("erp.in","r",stdin);
    freopen("perm.in","r",stdin); freopen("perm.out","w",stdout); 
    n=read(); mol=read();
    fac[0]=1; for(re int i=1;i<=n;i++) fac[i]=fac[i-1]*i%mol;
    inv[n]=qpow(fac[n],mol-2,mol); for(re int i=n;i>=1;i--) inv[i-1]=inv[i]*i%mol;
    dfs(0,1,1); sort(ks+1,ks+tot+1);
    for(re int i=1;i<=n;i++) {
        pre[i]=dik[i]=i;
        for(re int j=1;j<=6;j++) while(dik[i]%prim[j]==0) dik[i]/=prim[j];
    }
    sort(pre+1,pre+1+n,[](int a,int b){return dik[a]<dik[b];});
    f[0][1]=1;
    for(re int i=1,d=i;i<=n;i=d) {
        for(re int j=0;j<=n;j++) for(re int k=1;k<=tot;k++) g[j][k]=f[j][k];
        while(dik[pre[i]]==dik[pre[d]]) {
            int ls=pre[d]; d++;
            for(re int j=tot;j>=1;j--) {
                int lc=lcm(ks[j],ls/dik[ls]);
                if(lc>ks[tot]) continue;
                lc=lower_bound(ks+1,ks+1+tot,lc)-ks;
                for(re int a=n-ls;a>=0;a--) if(g[a][j]) {
                    int ds=1;
                        for(re int z=1;z*ls+a<=n;z++) {
                            (ds*=C(z*ls+a,ls))%=mol;
                            (g[z*ls+a][lc]+=g[a][j]*ds%mol*qpow(fac[ls-1],z,mol)%mol*inv[z]%mol)%=mol;
                        }
                }
            }
        }
        for(re int j=0;j<=n;j++) for(re int k=1;k<=tot;k++) (f[j][k]+=(g[j][k]-f[j][k])*dik[pre[i]]%mol*dik[pre[i]]%mol)%=mol;
    }
    int ans=0;
    for(re int i=1;i<=tot;i++) (ans+=f[n][i]*ks[i]%mol*ks[i]%mol)%=mol;
    (ans*=inv[n])%=mol;
    printf("%lld\n",ans);
}