POJ 2888 Magic Bracelet（burnside引理+矩阵）

题意：一个长度为n的项链，m种颜色染色每个珠子。一些限制给出有些颜色珠子不能相邻。旋转后相同视为相同。有多少种不同的项链？

思路：这题有点综合，首先，我们对于每个n的因数i，都考虑这个因数i下的不变置换个数，然后乘以(n/i)的欧拉函数加到ans上面，然后再让ans乘以n在模p下的逆元。至于怎么求因数i下的不变置换个数，相信大家都做过没有限制的，至于有限制的，大家可以考虑一下这样：初始数组a[m][m]都为1，对于每个限制x,y,都令a[x][y]=a[y][x]=0,我们有一个数列：b1,b2,b3,b4...bm，那么对于每个i∈[1,m-1]，都有a[bi][bi+1]=1，这样想到了什么？就是矩阵乘法，对于刚刚的矩阵a我们让s=a^i，然后把每个a[i][i]加起来，这样就代表：每个i出发又回到i的方法数，这样就变成因数i下的不变置换个数了！

  1 #include<algorithm>
  2 #include<cstdio>
  3 #include<cmath>
  4 #include<cstring>
  5 #include<iostream>
  6 const int Mod=9973;
  7 int a[105][105],b[105][105],c[105][105],s[105][105],n,m,K;
  8 int read(){
  9     char ch=getchar();
 10     int t=0,f=1;
 11     while (ch<‘0‘||ch>‘9‘){
 12         if (ch==‘-‘) f=-1;
 13         ch=getchar();
 14     }
 15     while (‘0‘<=ch&&ch<=‘9‘){
 16         t=t*10+ch-‘0‘;
 17         ch=getchar();
 18     }
 19     return t*f;
 20 }
 21 void exgcd(int a,int b,int &x,int &y){
 22     if (b==0){
 23         x=1;
 24         y=0;
 25         return;
 26     }
 27     exgcd(b,a%b,x,y);
 28     int t=x;
 29     x=y;
 30     y=(t-a/b*y);
 31 }
 32 int inv(int n){
 33     int A,B,x,y;
 34     A=n;B=Mod;
 35     exgcd(A,B,x,y);
 36     return (x+Mod)%Mod;
 37 }
 38 int euler(int n){
 39     int res=n,a=n;
 40     for (int i=2;i*i<=n;i++)
 41         if (a%i==0){
 42             res=res/i*(i-1);
 43             while (a%i==0) a/=i;
 44         }
 45     if (a>1) res=res/a*(a-1);
 46     return res%Mod;
 47 }
 48 void mult_matrix1(){
 49     for (int i=1;i<=m;i++)
 50      for (int j=1;j<=m;j++)
 51       c[i][j]=0;
 52     for (int i=1;i<=m;i++)
 53      for (int j=1;j<=m;j++)
 54       for (int k=1;k<=m;k++)
 55        c[i][j]=(c[i][j]+s[i][k]*b[k][j])%Mod;
 56     for (int i=1;i<=m;i++)
 57      for (int j=1;j<=m;j++)
 58        s[i][j]=c[i][j];
 59 }
 60 void mult_matrix2(){
 61     for (int i=1;i<=m;i++)
 62      for (int j=1;j<=m;j++)
 63       c[i][j]=0;
 64     for (int i=1;i<=m;i++)
 65      for (int j=1;j<=m;j++)
 66       for (int k=1;k<=m;k++)
 67        c[i][j]=(c[i][j]+b[i][k]*b[k][j])%Mod;
 68     for (int i=1;i<=m;i++)
 69      for (int j=1;j<=m;j++)
 70        b[i][j]=c[i][j];
 71 }
 72 int work(int x){
 73     for (int i=1;i<=m;i++)
 74      for (int j=1;j<=m;j++)
 75       b[i][j]=a[i][j],s[i][j]=0;
 76     int ans=0;
 77     for (int i=1;i<=m;i++)
 78       s[i][i]=1;
 79     while (x){
 80         if (x%2) mult_matrix1();
 81         mult_matrix2();
 82         x/=2;
 83     }
 84     for (int i=1;i<=m;i++) ans=(ans+s[i][i])%Mod;
 85     return ans;
 86 }
 87 int main(){
 88     int T=read();
 89     while (T--){
 90        n=read();m=read();K=read();
 91        for (int i=1;i<=m;i++)
 92          for (int j=1;j<=m;j++)
 93            a[i][j]=1;
 94        for (int i=1;i<=K;i++){
 95             int x=read(),y=read();
 96             a[x][y]=a[y][x]=0;
 97        }
 98        int ans=0;
 99        for (int i=1;i*i<=n;i++)
100         if (n%i==0){
101             ans=(ans+work(i)*euler(n/i))%Mod;
102             if (i*i!=n)
103             ans=(ans+work(n/i)*euler(i))%Mod;
104         }
105        ans=(ans*inv(n%Mod))%Mod;
106        printf("%d\n",ans);
107     }
108 }