给一串有n(1<=10^9 && gcd(n,9973))颗珠子的项链染色，已知某些{ci,cj}（i可能等于j）颜色不能出现在相连的珠子中。项链在平面内转动（不能翻转）得到的为等价方案，求不同的染色方案数mod 9973。
由于项链只能绕中心转动，而不能翻转，所以有一个很好的性质：项链转 2*PI*k/n个角度的置换共有 gcd(k,n)个长度均为len=n/gcd(k,n)的循环节，而且可以发现，若gcd(k,n)>1,则从某一个编号的项链开始顺时针（或逆时针）相邻的len个珠子两两处于不同的循环节中。这样k置换（1<=k<=n）的固定构型数等价于：
给长度为gcd(k,n)的序列染色（没写错，确实是序列），其中所有相邻编号（beg和end为相邻的）为合法对（若color i 和 color j 能够分别出现在相邻的两珠子上，那么说(i,j)是合法对）的染色方案。
假设一个长度为 L 的序列的合法表示为 a1,a2,……,aL 那么(a1,a2)(a2,a3),……，(a(L-1),aL),(aL,a1)必定都为1 ，所以令 g[ai,aj]=1表示（ai,aj）为合法对，g[ai,aj]=0表示(ai,aj)为非法对。那么以a1开头的合法序列为k个那么必定存在且仅存在k个不同的序列{b1,b2,b3,b4,……，b(L-1)}满足 g[a1,b1]*g[b1,b2]*……g[b(L-2),b(L-1)]*g[b(L-1),a1]=1; 这样很显然就是矩阵g[ai,aj]的L次幂中g[a1,a1]的的值。
很显然，这题可以用矩阵快速幂取模。
1 #include <iostream>
2 #include <cstring>
3 #include <cstdio>
4 using namespace std;
5
6 const int max_n=15;
7 const int MOD=9973;
8 int N;
9
10 struct Mat{
11 int mat[max_n][max_n];
12 };
13
14 Mat operator*(Mat a,Mat b){
15 Mat c;
16 memset(c.mat,0,sizeof(c.mat));
17 for(int i=0;i<N;i++){
18 for(int j=0;j<N;j++){
19 for(int k=0;k<N;k++){
20 if(a.mat[i][k] && b.mat[k][j]){
21 c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%MOD;
22 }
23 }
24 }
25 }
26 return c;
27 }
28
29 Mat operator^(Mat A,int x){
30 Mat c;
31 memset(c.mat,0,sizeof(c.mat));
32 for(int i=0;i<N;i++) c.mat[i][i]=1;
33 for(;x;x>>=1){
34 if(x&1){
35 c=c*A;
36 }
37 A=A*A;
38 }
39 return c;
40 }
41
42 int mypow(int a,int b){
43 int res=1;
44 for(;b;b>>=1){
45 if(b&1){
46 res*=a; res%=MOD;
47 }
48 a*=a; a%=MOD;
49 }
50 return res;
51 }
52
53 int M,K;
54
55 int PHI(int n){
56 int i,res=n; long long j;
57 for(i=2,j=4LL;j<=(long long)n;i++,j+=i+i-1){
58 if(!(n%i)){
59 res=res/i*(i-1);
60 while(!(n%i)) n/=i;
61 }
62 }
63 if(n>1) res=res/n*(n-1);
64 return res%MOD;
65 }
66
67 Mat A;
68 int solve(int p){
69 int res=0;
70 Mat C=A^p;
71 for(int i=0;i<N;i++){
72 res+=C.mat[i][i]; res%=MOD;
73 }
74 return res;
75 }
76 int main()
77 {
78 int T;
79
80 scanf("%d",&T);
81
82 for(int ncas=1;ncas<=T;ncas++){
83 scanf("%d%d%d",&M,&N,&K);
84 int u,v;
85 for(int i=0;i<N;i++){
86 for(int j=0;j<N;j++){
87 A.mat[i][j]=1;
88 }
89 }
90 for(int i=0;i<K;i++){
91 scanf("%d%d",&u,&v);
92 A.mat[u-1][v-1]=A.mat[v-1][u-1]=0;
93 }
94
95 int ans=0;
96
97 for(int p=1;p*p<=M;p++){
98 if(M%p==0){
99 if(p*p==M){
100 ans = (ans+ PHI(p)*solve(p) )%MOD;
101 }else{
102 ans = (ans+ PHI(p)*solve(M/p)+ PHI(M/p)*solve(p) )%MOD;
103 }
104 }
105 }
106
107 int inv=mypow((M%MOD),MOD-2); //晕，这里没取模导致超出整形，WA了三次
108 ans*=inv; ans%=MOD; //逆元。
109 printf("%d\n",ans);
110 }
111 return 0;
112 }
113