给一串有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