弱化版(其實完全一樣)
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();
}