弱化版(其实完全一样)
u1s1，洛谷上这题的第一个题解写得很不错，可以参考
直接边讲Polya定理边做这题
问题引入：n颗珠子组成的手串，每颗珠子有两种不同的颜色， 如果两个手串能够在旋转或翻转之后完全一样，就称它们是等价的，对手串的等价类计数。我们先把手串破环为链，两个长度为n的01序列等价当且仅当能够在循环移位或翻转后完全一样，求等价类数量。
注意两个序列"相等"指的是外表完全相同，"等价"指的是能够通过转化后外表完全相同。
Polya定理说这个问题的答案是：\(\frac 1{|G|}\sum_{g\in G}c(g,X)\)。其中X表示所有外表不同的序列的集合；G表示把所有变换看成置换之后的置换群(置换群的概念请自行搜索)；\(c(g,X)\)表示X中满足"将置换g施加到序列x上之后x的外表仍然不变"的元素x的个数，专业点说是不动点的个数。对于任意置换与等价类计数的问题，都可以用这个式子求。
具体证明可以看这里
对于这道题也是一样，这题中X是所有点有编号、点有权值(1-k)、边有权值(1/0表示存在或不存在)的n个点的无向完全图的集合；G是置换群(由于这题中的变换是给节点重新编号，所以G是所有1-n的排列的集合)；\(c(g,X)\)表示X中满足"将置换g施加到无向图x上之后x的外表仍然不变"的元素x的个数。
假设现在我们枚举g，考虑求出\(c(g,X)\)。对于g，如果我们把\(i\to g_i\)连有向边，是可以得到若干环的，假设这些环的大小从小到大排列为\(b_1,b_2\cdots b_m\)。对于其中的一条边\(u\to v\)，原图中的节点u在经过变换之后它的编号就要变成v了。考虑如果想让变换前后的图外表完全一样，这个图需要满足什么条件。首先g连出的每一个环内的点权值都必须相同，且所有的环必须覆盖[1,k]中的全部颜色(题目要求，用dp预处理方案数即可)。接下来的限制其实跟上面那个弱化版一样，每个长度为\(b_i\)的环内的所有边被分成了\(\lfloor \frac {b_i}2 \rfloor\)类，每一类的权值都必须相等；两个长度为\(b_i,b_j\)的环之间的\(b_ib_j\)条边被分成了\(gcd(b_i,b_j)\)类，每一类的权值都必须相等。
对于一个序列\(b_1,b_2\cdots b_m\)，我们已经能快速地用上面的方法算出它对答案的贡献，现在还要知道有多少个g对应这个序列，其实就是个简单的排列组合问题。令\(c_i\)表示b序列中值i出现的次数。对应这个b序列的g的个数为：\(\frac{n!}{\prod_{i=1}^m, b_i!}\cdot (\prod_{i=1}^m (b_i-1)!)\cdot \frac 1{\prod c_i! }\)，其中第一部分为多重组合数，用来选出每个环的元素；第二部分是把每个环中的所有元素排成有序环的方案数；第三部分是除掉相同的\(b_i\)算重的次数。
列一下最后答案的式子：\(ans=\sum_{b_1\cdots b_m} \frac 1{\prod b_i\prod c_i!}\cdot dp_{m,k}\cdot 2^{(\sum \lfloor \frac {b_i}2 \rfloor )+\sum_{i,j}gcd(b_i,b_j)}\)，其中\(dp_{i,j}\)表示1-i一共i个元素染色，且占了\([1,j]\)中的所有颜色的方案数。
我们直接暴搜枚举\(b\)数组所有可能的情况，然后用上面的方法暴力计算就行。时间复杂度\(O(能过)\)。
点击查看代码
#include <bits/stdc++.h>
#define rep(i,n) for(int i=0;i<n;++i)
#define repn(i,n) for(int i=1;i<=n;++i)
#define LL long long
#define pii pair <int,int>
#define pdd pair <double,double>
#define fi first
#define se second
#define mpr make_pair
#define pb push_back
void fileio()
{
#ifdef LGS
freopen("in.txt","r",stdin);
freopen("out.txt","w",stdout);
#endif
}
void termin()
{
#ifdef LGS
std::cout<<"\n\nEXECUTION TERMINATED";
#endif
exit(0);
}
using namespace std;
LL n,k,MOD,ans=0,fac[40],inv[40],rinv[40],dp[40][40];
vector <LL> d;
LL qpow(LL x,LL a)
{
LL res=x,ret=1;
while(a>0)
{
if(a&1) (ret*=res)%=MOD;
a>>=1;
(res*=res)%=MOD;
}
return ret;
}
void dfs(LL sum,LL mx)
{
if(sum==n)
{
LL res=1;
rep(i,d.size()) (res*=rinv[d[i]])%=MOD;
map <LL,LL> mp;rep(i,d.size()) ++mp[d[i]];
for(auto it:mp) (res*=inv[it.se])%=MOD;
(res*=dp[d.size()][k])%=MOD;
LL tot=0;
rep(i,d.size()) tot+=d[i]/2;
rep(i,d.size()) for(int j=i+1;j<d.size();++j) tot+=__gcd(d[i],d[j]);
(res*=qpow(2,tot))%=MOD;
(ans+=res)%=MOD;
return;
}
for(LL nxt=max(mx,1LL);nxt<=n-sum;++nxt)
{
d.pb(nxt);
dfs(sum+nxt,nxt);
d.pop_back();
}
}
int main()
{
fileio();
cin>>n>>k>>MOD;
dp[0][0]=1;
rep(i,n+3) rep(j,k+1) if(dp[i][j])
{
(dp[i+1][j]+=dp[i][j]*j)%=MOD;
(dp[i+1][j+1]+=dp[i][j]*(j+1))%=MOD;
}
fac[0]=1;repn(i,35) fac[i]=fac[i-1]*i%MOD;
rep(i,34) inv[i]=qpow(fac[i],MOD-2),rinv[i]=qpow(i,MOD-2);
dfs(0,0);
cout<<ans<<endl;
termin();
}