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 }
