[POJ 2888]Magic Bracelet[Polya Burnside 置换 矩阵]

也许更好的阅读体验
\(\mathcal{Description}\)

大意：给一条长度为\(n\)的项链，有\(m\)种颜色，另有\(k\)条限制，每条限制为不允许\(x,y\)颜色连在一起。要求有多少种本质不同的染色方式，本质不同的两种染色方式必须旋转不能互相得到。
输入方式：
第一行 \(t,\)表示t组数据
接下来\(t\)组数据：
每组数据第一行为\(n,m,k\)
接下来\(k\)行，每行两个数\(x,y\)表示不允许\(x,y\)颜色连在一起。
答案对9973取模
\((1 ≤ n ≤ 10^9, gcd(n, 9973) = 1), m (1 ≤ m ≤ 10), k (1 ≤ k ≤ m(m − 1) /2)\)
\(\mathcal{Solution}\)

本篇题解假设大家都会没有\(k\)条限制的版本

由于染色有限制，所以\(polya\)定理就不好用了
用\(burnside\)定理解决(发现大家标题都打得polya,于是也这么打标题了)

首先枚举不同的置换，即枚举循环长度，这一块和用\(polya\)定理有点像
由于\(n\)很大，要用\(\varphi\)来优化
一个置换中，循环的元素的所有颜色必须相同，所以我们要计算有多少个循环，这些循环有多少种染色方法

计算染色方法
考虑在该置换内一个循环一个循环的染色，我们可以看做从不同颜色节点走向另一节点
用邻接矩阵\(f[i][j]\)表示从\(i\)能否走向\(j\)
除去那\(k\)条限制，所有其他的\(i,j\) , \(f[i][j]=1\),即在染为第\(i\)种颜色后是可以染为第\(j\)种颜色的

这时候想到邻接矩阵的妙用，用矩乘计算\(n\)次之后得到的\(f[i][j]\)的意义发生改变
\(f[i][j]\)表示从\(i\)出发，走\(n\)步，最后结束在\(j\)有多少种走法。
由于是一个走一圈，所以 从\(i\)出发应在\(i\)结束，答案贡献为\(f[i][i]\)
这样就可以计算长度为\(n\)的循环的染色方法了：求出原矩阵的\(n\)次方后的矩阵后\(\sum_{i=1}^mf[i][i]\)

注意取模
代码

/*******************************
Author:Morning_Glory
LANG:C++
Created Time:2019年07月04日 星期四 14时25分13秒
*******************************/
#include <cstdio>
#include <fstream>
#include <cstring>
using namespace std;
const int maxn = 25;
const int mod = 9973;
//{{{cin 读入优化
struct IO{
    template<typename T>
    IO & operator>>(T&res){
        res=0;
        bool flag=false;
        char ch;
        while((ch=getchar())>'9'||ch<'0')    flag|=ch=='-';
        while(ch>='0'&&ch<='9') res=(res<<1)+(res<<3)+(ch^'0'),ch=getchar();
        if (flag)    res=~res+1;
        return *this;
    }
}cin;
//}}}
int t,n,m,k,x,y,ni,ans;
//{{{Matrix
struct Matrix{
    int rec[maxn][maxn];
    Matrix(){   memset(rec,0,sizeof(rec));}
    Matrix operator = (const Matrix &y){
        for (int i=1;i<=m;++i)
            for (int j=1;j<=m;++j)
                rec[i][j]=y.rec[i][j];
        return *this;
    }
    friend Matrix operator * (const Matrix &a,const Matrix &b){
        Matrix t;
        for (int i=1;i<=m;++i)
            for (int j=1;j<=m;++j)
                for (int k=1;k<=m;++k)
                    t.rec[i][j]=(t.rec[i][j]+a.rec[i][k]*b.rec[k][j])%mod;
        return t;
    }
    Matrix operator ^ (int b){
        Matrix s,a;
        a=*this;
        for (int i=1;i<=m;++i)  s.rec[i][i]=1;
        for (;b;b>>=1,a=a*a)
            if (b&1)    s=s*a;
        return s;
    }
}a,b;
//}}}Martix
//{{{get_phi
int get_phi (int x)
{
    int res=x;
    for (int i=2;i*i<=x;++i)
        if (x%i==0){
            res=res/i*(i-1);
            while (x%i==0)  x/=i;
        }
    if (x>1)    res=res/x*(x-1);
    return res%mod;//此处要取模
}
//}}}
//{{{ksm
int ksm (int a,int b)
{
    int s=1;
    a%=mod;
    for (;b;a=1ll*a*a%mod,b>>=1)
        if (b&1)    s=1ll*s*a%mod;
    return s;
}
//}}}
//{{{calc
int calc (int x)
{
    b=a^x;
    int res=0;
    for (int i=1;i<=m;++i)  res=(res+b.rec[i][i])%mod;
    return res;
}
//}}}
int main()
{
    cin>>t;
    while (t--){
        cin>>n>>m>>k;
        ans=0,ni=ksm(n,mod-2);//ni n的逆元
        for (int i=1;i<=m;++i)
            for (int j=1;j<=m;++j)
                a.rec[i][j]=1;
        while (k--){
            cin>>x>>y;
            a.rec[x][y]=a.rec[y][x]=0;
        }
        for (int i=1;i*i<=n;++i)
            if (n%i==0){
                ans=(ans+get_phi(i)*calc(n/i)%mod)%mod;
                if (i*i!=n) ans=(ans+get_phi(n/i)*calc(i)%mod)%mod;
            }
        ans=ans*ni%mod;
        printf("%d\n",ans);
    }
    return 0;
}