经典的有限制条件的Burnside计数+矩阵乘法！！！
对于这种限制条件的情况我们可以通过矩阵连乘得到，先初始化矩阵array[i][j]为1.如果颜色a和颜色b不能涂在相邻的珠子，
那么array[a][b] = array[b][a] = 0;
对于具有n/L个循环节的置换(L为每个循环节的长度)。先求出array[][]的n/L次幂，然后将这个矩阵的array[i][i]
(1<=i<=m)全部加起来即为这种置换下涂色不变的方法数。
代码如下：
1 #include<iostream> 2 #include<stdio.h> 3 #include<algorithm> 4 #include<iomanip> 5 #include<cmath> 6 #include<cstring> 7 #include<vector> 8 #define ll __int64 9 #define pi acos(-1.0) 10 #define MAX 50000 11 using namespace std; 12 int m,mod=9973,cnt,prime[50002]; 13 bool f[50004]; 14 struct Matrix 15 { 16 int a[15][15]; 17 }; 18 Matrix Mul(Matrix A,Matrix B) 19 { 20 Matrix ans; 21 for(int i=0;i<m;i++) 22 for(int j=0;j<m;j++){ 23 ans.a[i][j] = 0; 24 for(int k=0;k<m;k++) 25 { 26 if (A.a[i][k] && B.a[k][j]) 27 ans.a[i][j] = ans.a[i][j]+A.a[i][k]*B.a[k][j]; 28 } 29 ans.a[i][j]%=mod;//一定要在这里去摸，否则会超时的 30 } 31 return ans; 32 } 33 int pows(Matrix A,int b) 34 { 35 int i,j,re=0; 36 Matrix ans; 37 for(i=0;i<m;i++) 38 for(j=0;j<m;j++) 39 ans.a[i][j]=(i==j); 40 while(b){ 41 if (b&1) ans = Mul(A,ans); 42 b>>=1; 43 A = Mul(A,A); 44 } 45 for (i=0;i<m;i++) 46 re = (re+ans.a[i][i])%mod; 47 return re; 48 } 49 int Euler(int n) 50 { 51 int ans = n; 52 for (int i=2;i*i<=n;i++){ 53 if (n%i==0){ 54 ans=(ans/i)*(i-1); 55 n/=i; 56 while(n%i==0) 57 n/=i; 58 } 59 } 60 if(n>1) ans=(ans/n)*(n-1); 61 return ans%mod; 62 } 63 int mypow(int a,int b) 64 { 65 int ans = 1; 66 while(b){ 67 if (b&1) ans = (ans*a)%mod; 68 b>>=1; 69 a=(a*a)%mod; 70 } 71 return ans; 72 } 73 int main(){ 74 int i,n,sum,t,aa,bb,k,j; 75 scanf("%d",&t); 76 while(t--){ 77 scanf("%d%d%d",&n,&m,&k); 78 Matrix p; 79 for(i=0;i<m;i++) 80 for(j=0;j<m;j++) 81 p.a[i][j]=1; 82 for(i=0;i<k;i++){ 83 scanf("%d%d",&aa,&bb); 84 p.a[aa-1][bb-1]=p.a[bb-1][aa-1]=0; 85 } 86 sum = 0; 87 for(i=1;i*i<=n;i++){ 88 if(n%i==0){ 89 sum = (sum+pows(p,i)*Euler(n/i))%mod; 90 if(i*i!=n) 91 sum = (sum+pows(p,n/i)*Euler(i))%mod; 92 } 93 } 94 int x=mypow(n%mod,mod-2); 95 sum = (sum*x)%mod; 96 printf("%d ",sum); 97 } 98 return 0; 99 }