大致题意： 求含(n)个点的互不同构的简单无向图个数。
转化
考虑把一条边是否存在看成黑白染色，这就变成了一个经典的(Polya)计数问题。
应该可以看作是无向图(Polya)定理的板子吧，但确实挺难的。
无向图(Polya)定理
本题中置换群的对象就是这(frac{n(n-1)}2)条边，但由于置换的数量达到(n!)，显然不可以裸暴力。
假设点的置换是(i ightarrow P_i)，则边的置换就是((i,j) ightarrow(P_i,P_j))。
考虑两种置换循环节个数的关系，发现：
- 若点(i,j)属于同一长度为(x)的循环中，((i,j))组成的置换中循环节个数为(lfloorfrac x2 floor)。
- 若点(i,j)分别属于长度为(x,y)的两个循环中，((i,j))组成的置换中循环节个数为(gcd(x,y))。
因此，如果一个点置换长度分别为(L_1,L_2,...,L_t)，则边置换的循环节总数就是：
[sum_{i=1}^tlfloorfrac{L_i}2
floor+sum_{i=1}^tsum_{j=i+1}^tgcd(L_i,L_j)
]
那么就需要求满足循环节长度为(L_1,L_2,...,L_t)的边置换数目，也就是把(1sim n)放入大小分别为(L_1,L_2,...,L_t)的环的方案数。
如果(L)互不相同，方案数就是：
[frac{n!}{prod_{i=1}^tL_i}
]
而若可能有相同的(L)，同样大小的环就会造成重复贡献。
因此我们将(L_i)去重，用(C_i)表示第(i)种值的个数，方案数就是：
[S=frac{n!}{prod_{i=1}^t(L_i)^{C_i} imes C_i!}
]
在这种表示方法下，之前的循环节总数就应该是：
[T=sum_{i=1}^tC_i imes lfloorfrac{L_i}2
floor+sum_{i=1}^t(frac{C_i imes (C_i-1)}2 imes L_i+sum_{j=t+1}^n C_i imes C_j imes gcd(L_i,L_j))
]
最终答案就是：（注意这其实是一个(Polya)定理的形式，只是在(m^{c(g_i)})之前乘上了一个方案数）
[ans=frac1{n!}sum S imes 2^T
]
这道题的核心思想就是把相似的情况放在一起讨论，除去冗余的方法在(Polya)计数问题中是一个很有用的优化。
代码
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 60
#define X 997
using namespace std;
int n,ans,t,l[N+5],c[N+5],pw[X],Inv[N+5],Fac[N+5],IFac[N+5],g[N+5][N+5];
I int QP(RI x,RI y) {RI t=1;W(y) y&1&&((t*=x)%=X),(x*=x)%=X,y>>=1;return t;}
I int gcd(CI x,CI y) {return y?gcd(y,x%y):x;}
I void Calc()//计算贡献
{
RI i,j,S=Fac[n],T=0;for(i=1;i<=t;++i) (T+=c[i]*(l[i]/2))%=X-1;//T统计总长度：同一循环
for(i=1;i<=t;++i) for((T+=c[i]*(c[i]-1)/2*l[i])%=X-1,j=i+1;j<=t;++j) (T+=c[i]*c[j]*g[l[i]][l[j]])%=X-1;//T统计总长度：不同循环
for(i=1;i<=t;++i) (S*=QP(Inv[l[i]],c[i])*IFac[c[i]])%=X;(ans+=S*pw[T])%=X;//S统计方案数，ans统计答案
}
I void dfs(CI x,CI s)//搜索划分方案
{
if(s==n) return Calc();if(x==1) return (void)(l[++t]=1,c[t]=n-s,Calc(),--t);//边界
dfs(x-1,s),++t;for(RI i=1;x*i<=n-s;++i) l[t]=x,c[t]=i,dfs(x-1,s+x*i);--t;//枚举C
}
int main()
{
RI i,j;for(scanf("%d",&n),pw[0]=i=1;i^X;++i) pw[i]=(pw[i-1]<<1)%X;//预处理2的幂
for(Fac[0]=i=1;i<=n;++i) Fac[i]=Fac[i-1]*i%X;//预处理阶乘
for(IFac[n]=QP(Fac[n],X-2),i=n-1;~i;--i) IFac[i]=IFac[i+1]*(i+1)%X;//预处理阶乘逆元
for(i=1;i<=n;++i) for(Inv[i]=QP(i,X-2),j=1;j<=n;++j) g[i][j]=gcd(i,j);//预处理逆元和gcd
return dfs(n,0),printf("%d
",IFac[n]*ans%X),0;//dfs后输出答案
}