poj2888 Magic Bracelet
给你一个正n（<10^9）边形和m（<10）种色料，要求给正n边形顶点染色并且规定k组颜色对不能相邻，
输入保证n与mod互质，计数染色总方案数（绕图形中心旋转后相同的方案算一种）对mod取模。
问题的关键在于保证给定的颜色对不相邻，而如果我们用传统的容斥来排除不合法的染色方案在题目给定的n的规模下是不具有可操作性的。
考虑n边形的σi旋转置换，显然有循环：
0 : 0 -> i % n -> 2 *i % n ->... -> 0
1 : 1 -> (i + 1) % n -> (i * 2 + 1) % n ->... -> 1
...
d - 1 : d - 1 -> (i + d - 1) % n -> (i * 2 + d - 1) % n ->... d - 1
即存在d个循环，考虑到每个点的位置是等价的且恰处于一个循环节中,有
各循环节长度相等，其值为l = n / d。
第j + 1个循环中的点p满足 p Ξ j (mod n), 因此若(j + d1 * i ) % n = j,
有d1 * i Ξ 0 (mod n), d1 的最小值为lcm(i, n) / i = n * gcd(i, n)即循环节长度l。
解出d = gcd(i, n)。
显然每个循环节内染色相同，假设第i个循环节染色c(i)。
则n边形染色:c(0) -> c(1) ->...-> c(d - 1) -> c(0) -> c(1) ->....c(d - 1)。
这等价于我们用m种颜色对d 边形染色，同时保证染色合法。
也就是说保证染色：(c(0)c(1)) (c(1)c(2))...(c(d-1)c(0))中成对颜色是合法的。（*）
即便到此用容斥解决也是不可行的，需要将其实现到可快速计算的载体中。
令f(p, q)=
1，若p, q色料可以相邻。
0，其他。
则（*）等价于f(c(0), c(1)) * ... * f(c(d-1), c(0)) = 1。
也就是说每种可行方案为该连乘式贡献1。
考虑布尔矩阵F：
F^n（i，j）表示从i到j路径长度为n的路径总数（这里假设任两点之间距离为1）。
如此所求方案数即∑φ(n / d) * F^d * Inversion(n, mod) % mod (d | n)。
http://poj.org/problem?id=2888
1 #include <cstdio> 2 #include <cstring> 3 #include <algorithm> 4 #include <cmath> 5 using namespace std; 6 const int mod = 9973; 7 const int maxn = 12; 8 int n, m, k, n1; 9 int prime[50], k1; 10 int p2[50], k2; 11 int cnt[50]; 12 int ans; 13 14 struct Matrix{ 15 int matrix[maxn][maxn]; 16 }B, C; 17 18 Matrix operator * (Matrix A, Matrix B){ 19 Matrix C; 20 memset(C.matrix, 0, sizeof C.matrix); 21 for(int i = 0; i < m; i++){ 22 for(int j = 0; j < m; j++){ 23 for(int u = 0; u < m; u++){ 24 int p = A.matrix[i][u], q = B.matrix[u][j]; 25 if(p && q) C.matrix[i][j] += p * q; 26 } 27 C.matrix[i][j] %= mod; 28 } 29 } 30 return C; 31 } 32 33 Matrix operator ^ (Matrix A, int p){ 34 Matrix C; 35 memset(C.matrix, 0, sizeof C.matrix); 36 for(int i = 0; i < m; i++) C.matrix[i][i] = 1; 37 while(p){ 38 if(p & 1) C = C * A; 39 p >>= 1; 40 A = A * A; 41 } 42 return C; 43 } 44 45 int phi(int num){ 46 k2 = 0; 47 for(int i = 0; i < k1; i++) if(num % prime[i] == 0) num /= prime[i], p2[k2++] = prime[i]; 48 for(int i = 0; i < k2; i++) num *= p2[i] - 1; 49 return num % mod; 50 } 51 52 int power(int a, int p){ 53 a %= mod; 54 int ans = 1; 55 while(p){ 56 if(p & 1) ans *= a, ans %= mod; 57 p >>= 1; 58 a *= a; 59 a %= mod; 60 } 61 return ans; 62 } 63 64 void dfs(int pos, int num){ 65 if(pos == k1){ 66 C = B ^ num; 67 int ans1 = 0; 68 for(int i = 0; i < m; i++) ans1 += C.matrix[i][i], ans1 %= mod; 69 ans1 = (ans1 * phi(n / num)) % mod; 70 ans = (ans + ans1) % mod; 71 return; 72 } 73 int tem = num; 74 for(int i = 0; i <= cnt[pos]; i++){ 75 dfs(pos + 1, tem); 76 tem *= prime[pos]; 77 } 78 } 79 80 void solve(){ 81 n1 = n; 82 k1 = 0; 83 memset(cnt, 0, sizeof cnt); 84 if(!(n & 1)){ 85 prime[k1++] = 2; 86 while(!(n & 1)) n >>= 1, ++cnt[k1 - 1]; 87 } 88 int mid = (int)sqrt(n); 89 for(int i = 3; i <= mid; i += 2){ 90 if(n % i == 0){ 91 prime[k1++] = i; 92 while(n % i == 0) n /= i, ++cnt[k1 - 1]; 93 mid = (int)sqrt(n); 94 } 95 } 96 if(n != 1) prime[k1++] = n, ++cnt[k1 - 1]; 97 n = n1; 98 ans = 0; 99 dfs(0, 1); 100 ans = ans * power(n, mod - 2) % mod; 101 printf("%d\n", ans); 102 } 103 104 int main(){ 105 //freopen("in.txt", "r", stdin); 106 int T; 107 scanf("%d", &T); 108 while(T--){ 109 scanf("%d%d%d", &n, &m, &k); 110 for(int i = 0; i < m; i++) for(int j = 0; j < m; j++) B.matrix[i][j] = 1; 111 for(int i = 0, p, q; i < k; i++){ 112 scanf("%d%d", &p, &q); 113 --p, --q; 114 B.matrix[p][q] = B.matrix[q][p] = 0; 115 } 116 solve(); 117 } 118 return 0; 119 }
