[SHOI2006]有色图

感觉polya定理还没搞太清,写点东西助于理解.
首先注意到这题求的是本质不同的染色方案,这个好像很符合polya定理求的东西.
\[l=\frac{1}{|G|}\sum_{i\in G}m^{w(i)}\]
其中\(w(i)\)表示的是i置换下的循环节个数
要用这个定理首先必须要找置换.发现边的置换很难找,于是考虑寻找点的置换和边的置换的联系.
如果我们暴力枚举点的置换.再把置换写成循环的形式,这样这些点就被我们分成一组一组的了.每个组就是一个循环.
一条边是连在两个点之间的,我们再来考虑在点的置换确定的条件下,边的置换个数的循环个数.
这个个数怎么求呢?
分两种情况讨论:

case1 两个点在一个循环内

现在想象l个点它们围成一个环,它们之间连满了边,边数是\(\left(\begin{matrix}l \\2\end{matrix}\right)\),如果环的元素是奇数个的话,你转l次就可以转回原来的样子,点一一对应,边也一一对应.那么边的循环长度显然就为l了,然后又有\(\frac{l(l-1)}{2}\)条边,那么循环节个数就是\(\frac{l-1}2\)了.
如果l为偶数,那么就会有一些边,你转180度,就可以转回来,它们的循环长度为\(\frac l 2\).
这些边有多少个呢?显然是\(\frac{l}{2}\)个,我们愉快地减掉它,就好了.循环个数:\[\frac{\frac{l(l-1)}{2}-\frac{l}2}{l}+\frac{\frac{l}{2}}{\frac{l}2}=\frac l 2\]

case2 两个点在不同循环内

同样用上面的方法,根据循环长度,计算循环个数,假设这两个点所在的点循环的长度分别为\(l_1,l_2\).那么它们合成的循环长度就为\(lcm(l_1,l_2)\),然后边数显然为\(l_1*l_2\)这就很美滋滋了.循环个数就为\(\frac{l_1*l_2}{lcm(l_1,_2)}=gcd(l_1,l_2)\)
所以这样以来,已知点的置换,得到的边置换的循环个数就为:
\[c=\sum_{i=1}^k\lfloor\frac{l_i}{2}\rfloor+\sum_{i=1}^k\sum_{j=i+1}^ngcd(l_i,l_j)\]

但是这样显然是不行的.
很容易发现对称群个数\(|S_n|=n!\)这个显然是不可以接受的.那我们想办法简化状态,之前的状态里包含了每个循环里具体有哪些点的信息,但是这些信息是多余的,我们只要知道每个循环里有多少个元素就可以了,那么我们原来枚举的东西,就可以代替成枚举整数的拆分.
但是光知道一种状态下的答案是不够的,我们还得计算有多少种点的置换对应这种状态.
设被当前状态为:被分成了k个循环,每个循环的元素个数为\(l_i\).
我们算的是状态对应的点的置换数,想一想发现,一个循环里点的顺序不同代表的置换是不同的,如果固定第一个点的话,后面的不同排列个数为\((l_i-1)!\)
所以如果每个l都不相同的话,那么\(S=\frac{n!}{\prod_{i=1}^k l_i!}*\prod_{i=1}^k (l_i-1)!=\frac{n!}{\prod_{i=1}^k l_i}\)
如果有相同的,那么设\(q_i\)为第i种l的个数.
所以\[S=\frac{n!}{(\prod_{i=1}^kl_i)*(\prod_{i=1}^{kind}q_i)}\]
这里可以用费马小定理求逆元.
\[\therefore ans=\sum S*m^c\]

// luogu-judger-enable-o2
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#define maxn 105
#define ll long long
using namespace std;
ll mod,f[maxn],ans;
int n,m,rec[maxn];
ll gcd(ll a,ll b)
{if(b==0)return a;return gcd(b,a%b);
}
ll ksm(ll a,ll b)
{ll res=1;while(b){if(b&1)res=(res*a)%mod;a=a*a%mod;b>>=1;}return res;
}
void calc(int x)
{ll sum=0,mul=1;int now=1;for(int i=1;i<=x;i++)sum+=rec[i]/2;for(int i=1;i<=x;i++)for(int j=i+1;j<=x;j++)sum+=gcd(rec[i],rec[j]);for(int i=1;i<=x;i++)(mul*=rec[i])%=mod;for(int i=2;i<=x;i++){if(rec[i]!=rec[i-1])mul=mul*f[now]%mod,now=0;now++;}mul=mul*f[now]%mod;mul=f[n]*ksm(mul,mod-2)%mod;//cout<<n<<" "<<mul<<endl;ans=(ans+mul*ksm(m,sum)%mod)%mod;
}
void dfs(int k,int x,int s)
{if(!x)calc(k-1);if(x<s)return;for(int i=s;i<=x;i++){rec[k]=i;dfs(k+1,x-i,i);}
}
int main()
{cin>>n>>m>>mod;f[0]=1;for(int i=1;i<=n;i++)f[i]=f[i-1]*i%mod;dfs(1,n,1);ans=ans*ksm(f[n],mod-2)%mod;printf("%lld",ans);return 0;
}