前置知识
\(60\space\text{pts}\) 解法
就像对于一个数，我们经常从素因子之积的角度看待它一样，在这道题中，我们从轮换的角度看待置换。
我们考虑一个轮换变成恒等变换所需次数：一个长度为 \(l\) 轮换，可以看做一个边数为 \(l\) 的环，置换乘法可以看做数字沿着边转一圈。假设某一条边的终点处的数字一开始为 \(k\)，那么这条边的起点的位置编号一定是 \(k\)。以轮换 \((3,1,2)\) 为例，边 \(3\rightarrow1\)（位置编号）为例，起点的位置编号为 \(3\)，终点一开始是 \(3\)。也就是说，一个数字从它开始的位置到达恒等变换的位置（位置编号也为这个数字），需要经过 \(l-1\) 条边，反映到幂次上，贡献为 \(l\)。
而对于整个置换，若变成恒等置换，则贡献为每个轮换长度的最小公倍数。答案即为每一种最小公倍数的平方与置换总数，即 \(n!\) 的商。
我们考虑统计出每一种最小公倍数的出现次数，设 \(dp_{i,j}\) 为将 \(i\) 个位置分配给若干个长度之最小公倍数为 \(j\) 的方案数。
先来考虑局部转移：将 \(dp_{i,j}\) 转移给 \(dp_{i+k\cdot s,\operatorname{lcm}(j,s)}\)，即在已经分配的 \(i\) 个位置中再加入 \(k\) 个长度为 \(s\) 的轮换。对于位置的分配，可以转化为多重集 \(\{i\cdot a_1,s\cdot a_2, \cdots ,s\cdot a_{k+1}\}\) 的排列数，即 \(\frac{(i+k\cdot s)!}{{s!}^k i!}\)。而对于编号的分配，也就是置换环的边“交错情况”的分配，考虑第一个点有 \((n-1)\) 种终点，第二个有 \((n-2)\) 种，以此类推，单个置换环有 \((s-1)!\) 种，共 \(((s-1)!)^k\) 种。这一项与上一项合并得到 \(\frac{(i+k\cdot s)!}{s^k i!}\)。此外，由于每个置换环位置确定后内容也确定，需要消去几个置换环仅仅互换位置的情况，即除以 \(k!\)，得到：
同时看到，在整体的转移过程中，得到 \(dp_{i,j}\) 时必有分子 \(i!\)，这必将在下一步被消去，留下的 \(n!\) 也将在统计答案时因平均被消去，所以每次转移的系数被消为 \(\frac{1}{s!^k k!}\) ，这样最后不必再除 \(n!\)。
部分分代码
#include<bits/stdc++.h>
using namespace std;
unordered_map<long long,long long> dp[300];
unordered_map<long long,bool> vis[300];
long long n,p,t[300],tmaxn[300],inv[100100],invt[300];
long long gcd(long long a,long long b){
long long tem=0;
for(;a;){
tem=a;
a=b%a;
b=tem;
}
return b;
}
long long lcm(long long a,long long b){
return a*b/gcd(a,b);
}
long long qpow(long long x,long long y){
if(!y)
return 1;
long long tem=qpow(x,y>>1);
if(y&1){
return x*tem%p*tem%p;
}else{
return tem*tem%p;
}
}
long long finv(long long x){
if(x<min(p,1ll*100100))
return inv[x];
else
return qpow(x,p-2);
}
void dfs(long long x,long long y){
if(!x){
if(y>=100100)
printf("%lld\n",y);
return;
}
for(int i=1;i<=x;++i){
dfs(x-i,lcm(y,i));
}
}
int main(){
freopen("perm.in","r",stdin);
freopen("perm.out","w",stdout);
scanf("%lld%lld",&n,&p);
inv[0]=1;
inv[1]=1;
for(int i=2;i<min(p,1ll*100100);++i){
inv[i]=(p-p/i)*inv[p%i]%p;
}
t[0]=1;
for(int i=1;i<=n;++i){
t[i]=t[i-1]*i%p;
}
invt[n]=qpow(t[n],p-2);
for(int i=n-1;i;--i){
invt[i]=invt[i+1]*(i+1)%p;
}
invt[0]=1;
long long ans=0;
dp[0][1]=1;
vis[0][1]=1;
for(int s=1;s<=n;++s){
for(int i=n;i>=0;--i){
for(int k=(n-i)/s;k;--k){
for(auto j:vis[i]){
if(!dp[i][j.first])
continue;
vis[i+k*s][lcm(j.first,s)]=1;
long long tem=dp[i][j.first]*finv(qpow(s,k))%p*invt[k]%p;
dp[i+k*s][lcm(j.first,s)]=(dp[i+k*s][lcm(j.first,s)]+tem)%p;
}
}
}
}
for(auto i:vis[n]){
ans=(ans+dp[n][i.first]*(i.first)%p*(i.first)%p)%p;
}
printf("%lld",ans);
return 0;
}
正解
我们发现由于 \(\operatorname{lcm}\) 结果很多，上面做法会超时，考虑如果某个大于 \(\sqrt{n}\) 的因子在 \(\operatorname{lcm}\) 中出现的次数大于 \(1\) 则这个长度必定大于等于这个因子的平方，也大于 \(n\)。因此，每种这样的因子最多只会出现一次。我们把剩下的因子结合起来，将每种 \(\operatorname{lcm}\) 转换为这几个因子的 \(\operatorname{lcm}\) 复合上其它因子的出现情况。每次将前者作为状态，将后者相同的同时转移，枚举复合上的前者，正常转移后将 \(dp\) 的变化量乘上后者的平方(这样后面枚举的值会和前面复合)，统计答案时再复合上前者平方即可。
正解
#include<bits/stdc++.h>
using namespace std;
const long long pr[6]={2,3,5,7,11,13};
struct node{
int pos,lft;
friend bool operator < (node x,node y){
return x.lft<y.lft;
}
}data[300];
long long n,p,base[3000],bcnt,t[300],dp[300][9000],up[300][9000],maxn[300],inv[100100],invt[300];
long long gcd(long long a,long long b){
long long tem=0;
for(;a;){
tem=a;
a=b%a;
b=tem;
}
return b;
}
long long lcm(long long a,long long b){
return a*b/gcd(a,b);
}
long long qpow(long long x,long long y){
if(!y)
return 1;
long long tem=qpow(x,y>>1);
if(y&1){
return x*tem%p*tem%p;
}else{
return tem*tem%p;
}
}
void dfs(int now,long long tim,long long sum){
if(now==6){
++bcnt;
base[bcnt]=tim;
return;
}
dfs(now+1,tim,sum);
long long tem=1;
for(;;){
tem*=pr[now];
if(sum+tem>n)
break;
dfs(now+1,tim*tem,sum+tem);
}
}
int main(){
freopen("perm.in","r",stdin);
freopen("perm.out","w",stdout);
scanf("%lld%lld",&n,&p);
dfs(0,1,0);
sort(base+1,base+bcnt+1);
for(int i=1;i<=n;++i){
data[i].pos=i;
data[i].lft=i;
for(int j=0;j<6;++j){
for(;!(data[i].lft%pr[j]);){
data[i].lft/=pr[j];
}
}
}
sort(data+1,data+n+1);
inv[0]=1;
inv[1]=1;
for(int i=2;i<min(p,1ll*100100);++i){
inv[i]=(p-p/i)*inv[p%i]%p;
}
long long ans=0;
dp[0][1]=1;
for(int lft=1,now=1;lft<=n;lft=now){
for(int i=0;i<=n;++i){
for(int j=1;j<=bcnt;++j){
up[i][j]=dp[i][j];
}
}
long long timr=data[lft].lft*data[lft].lft%p;
for(;data[lft].lft==data[now].lft;++now){
int s=data[now].pos;
for(int j=bcnt;j;--j){
int temlcm=lower_bound(base+1,base+bcnt+1,lcm(base[j],data[now].pos/data[now].lft))-base;
if(temlcm>bcnt)
continue;
for(int i=n-s;i>=0;--i){
if(!up[i][j])
continue;
int temi=i+s;
long long div=inv[s];
for(int k=1;temi<=n;++k,temi+=s,div=div*inv[s]%p*inv[k]%p){
up[temi][temlcm]=(up[temi][temlcm]+up[i][j]*div)%p;
}
}
}
}
for(int i=0;i<=n;++i){
for(int j=1;j<=bcnt;++j){
dp[i][j]=(dp[i][j]+timr*(up[i][j]-dp[i][j]))%p;
}
}
}
for(int i=1;i<=bcnt;++i){
ans=(ans+dp[n][i]*base[i]%p*base[i])%p;
}
ans=(ans+p)%p;
printf("%lld",ans);
return 0;
}