bzoj 1815: [Shoi2006]color 有色图 置换群

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
using namespace std;
#define PROB "graph"
#define MAXE MAXN*MAXN
#define MAXN 66
#define MAXM MAXN*MAXN
typedef long long qword;
int n,m,mod;
struct edge
{
        int x,y;
}e[MAXE];
bool operator < (edge e1,edge e2)
{
        if (e1.x==e2.x)return e1.y<e2.y;
        return e1.x<e2.x;
}
bool vis[MAXN];
int vec[MAXN];
qword sum=0;
int per[MAXM];
int g[MAXN];
qword pow_mod(qword x,qword y)
{
        qword ret=1;
        while (y)
        {
                if (y&1)ret=ret*x%mod;
                x=x*x%mod;
                y>>=1;
        }
        return ret;
}
void init()
{
        int i,j,k;
        int m;
        int x,y;
        for (i=2;i<=n;i++)
        {
                m=0;
                for (j=0;j<i;j++)
                        for (k=j+1;k<i;k++)
                                e[m].x=j,e[m++].y=k;
                for (j=0;j<m;j++)
                {
                        edge et;
                        et.x=min((e[j].x+1)%i,(e[j].y+1)%i);
                        et.y=max((e[j].x+1)%i,(e[j].y+1)%i);
                        per[j]=lower_bound(e,e+m,et)-e;
                }
                for (j=0;j<m;j++)
                {
                        if (~per[j])
                        {
                                g[i]++;
                                x=j;
                                y=per[x];
                                per[x]=-1;
                                while (~per[y]){
                                        x=y;
                                        y=per[x];
                                        per[x]=-1;
                                }
                        }
                }
        }
        return ;
}
int gcd(int x,int y)
{
        return x%y==0?y:gcd(y,x%y);
}
qword fact[MAXN],ufact[MAXN],uval[MAXN];
int gcdv[MAXN][MAXN];
int totv;
void dfs(int s,int b)
{
        if (!s)
        {
                register qword val=1;
                register int res=0,n0=n,i,j;
                for (i=0;i<totv;i++)
                        for (j=i+1;j<totv;j++)
                                res=(res+gcdv[vec[i]][vec[j]])%(mod-1);
                for (i=0;i<totv;i++)
                {
                        res=(res+g[vec[i]])%(mod-1);
                        val=val*fact[n0]%mod*ufact[n0-vec[i]]%mod*ufact[vec[i]]%mod;
                        n0-=vec[i];
                }
                //for (int i=0;i<totv;i++)printf("%d ",vec[i]);printf("\n");
                for (i=0;i<totv;)
                {
                        for (j=i;j<totv && vec[j]==vec[i];j++);
                        val=val*ufact[j-i]%mod;
                        i=j;
                }
                for (i=0;i<totv;i++)
                        val=val*fact[vec[i]]%mod*uval[vec[i]]%mod;
                //printf("val:%d\n",val);
                //printf("res:%d\n",res);
                sum+=pow_mod(m,res)*val%mod;
                sum%=mod;
                return ;
        }
        for (int i=b;i<=s;i++)
        {
                vec[totv++]=i;
                dfs(s-i,i);
                totv--;
        }
}
int main()
{
        freopen("input.txt","r",stdin);
        int i,j,k,x,y,z;
        scanf("%d%d%d",&n,&m,&mod);
        init();
        fact[0]=1;
        for (i=1;i<=n;i++)
                for (j=1;j<=n;j++)
                        gcdv[i][j]=gcd(i,j);
        for (i=1;i<=n;i++)
                fact[i]=fact[i-1]*i%mod;
        for (i=0;i<=n;i++)
                ufact[i]=pow_mod(fact[i],mod-2);
        for (i=1;i<=n;i++)
                uval[i]=pow_mod(i,mod-2);
        dfs(n,1);
        printf("%lld\n",sum*ufact[n]%mod);
}