题解

由双射可以想到置换，置换一类的问题通常用polya定理解决，若对边进行置换，推不动，便换一个思路，对点进行置换，而对边进行染色，边分为存在和不存在两种。

\[ p\acute{o}lya定理：L=\frac{1}{\left| G \right|}(m^{c(g_1)}+m^{c(g_2)}+...+m^{c(g_n)}) \]

在此题中：
\[ \left| G \right|=n!,m=2,c(g_i)为循环节个数 \]
对于一个点的置换
\[ a_1,a_2...a_n \]
可以分成m个循环节(b为每个循环节的大小)
\[ b_1,b_2...b_m \]
对于bi，循环节内的边数为
\[ {b*(b-1)\over 2} \]
对于边(i,j)，i在b中循环一周，j也在b中循环一周，一个循环的长度为b，故循环节个数为
\[ \left \lceil \frac{b-1}{2} \right \rceil=\left \lfloor \frac{b}{2} \right \rfloor \]
对于每个循环节之间的边，若两个循环节大小分别为b1,b2，则总边数为b1b2，边(i,j)，i循环一周要b1次，j循环一周要b2次，故循环长度为lcm(b1,b2)，循环节个数为
\[ \frac{b_1b_2}{lcm(b_1,b_2)}=gcd(b_1,b_2) \]
无向图总个数为：\(n^{n-2}\)，直接枚举复杂度爆炸，则对于b相同的置换，统计此类置换的个数，对于b，若大小为b的循环节个数有N个，总数为sum(=bN)，处理到b时，总共还有tot个a未确定，置换个数为：
\[ {tot \choose sum}\frac{sum!}{b^N N!}\qquad（/b因为环形排列中相同循环个数，N!因为N个循环块之间的排列顺序与答案无关） \]
枚举b时从大到小或从小到大枚举，即减少枚举次数，又避免排序带来的复杂度。

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#define mod 100000007
using namespace std;
int n,d[60];
long long f[60],inv[60],g[60][60],p[60],c[60][60],ans=0;
int gcd(int a,int b){return (!b)?a:gcd(b,a%b);}
long long Pow(long long A,int T)
{
    long long S=1ll;
    while(T)
    {
        if(T&1)S=(S*A)%mod;
        A=(A*A)%mod,T>>=1;
    }
    return S;
}
void dfs(int x,int mx,int sum)
{
    if(sum==0)
    {
        long long tmp=1ll,T=1ll;
        d[x+1]=0;
        for(int i=1;i<=x;i++)tmp=(tmp*p[d[i]])%mod;
        for(int i=1;i<=x;i++)for(int j=i+1;j<=x;j++)tmp=(tmp*g[d[i]][d[j]])%mod;
        for(int i=1,tot=n;i<=x;)
        {
            int pos=i,S=0;
            while(d[pos]==d[pos+1]&&pos<=x)pos++;
            S=(pos-i+1)*d[i],tmp=(tmp*c[tot][S])%mod,tot-=S;
            long long a1=Pow(1ll*d[i],pos-i+1),a2=a1*f[pos-i+1]%mod,a3=f[S]*Pow(a2,mod-2)%mod;
            tmp=(tmp*a3)%mod,i=pos+1;
        }
        ans=(ans+tmp)%mod;
        return;
    }
    for(int i=min(mx,sum);i>=1;i--)d[x+1]=i,dfs(x+1,i,sum-i);
}
int main()
{
    scanf("%d",&n);
    f[0]=1ll;
    for(int i=1;i<=n;i++)f[i]=f[i-1]*1ll*i%mod;
    inv[n]=Pow(f[n],mod-2);
    for(int i=n-1;i>=0;i--)inv[i]=inv[i+1]*1ll*(i+1)%mod;
    for(int i=1;i<=n;i++)for(int j=1,x;j<=n;j++)g[i][j]=Pow(2ll,gcd(i,j));
    for(int i=1;i<=n;i++)p[i]=Pow(2ll,i/2);
    c[1][0]=c[1][1]=1ll;
    for(int i=2;i<=n;i++)
    {
        c[i][0]=1ll;
        for(int j=1;j<=i;j++)c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
    }
    dfs(0,n,n);
    printf("%lld\n",(ans*inv[n])%mod);
    return 0;
}