poj2888

题意：给定N（10^9）种珠子及颜色总数M，并给出一些限制条件，如颜色i与j不能相邻，求总共有多少种方案数，结果Mod 9973

思路：假设如果原题没有颜色限制的话，那就是一题polya的模板题，再加上phi（欧拉函数）的优化即可（雷同poj2154）

但是因为有了所谓的限制，那就不能直接套，我么该怎么办呢？

这我们很容易想到先前讲的那个以burnside的思想，polya的置换数的方法求，那怎么做？

我们把每个循环节数计算出来，那么在当前下，处于同一个环上的必须颜色相同，那么同一个环就可以看成一个点，然后对这些点进行染色，满足相邻的限制即可，这样的方案数就是当前置换下的方案数。举个例子：

比如：n = 6， M = 3 ，限制为1与2不相邻

那么在旋转2个位置情况下，可得到置换：（1 3 5）（2 4 6）

那么就等价于用3种颜色给2个位置涂色，并满足颜色1与颜色2不相邻

这样的话就可以用矩阵来解决方案数的问题，把每种颜色间连一条边，如果颜色不能涂在相邻处便为0，否者为1，

设这个矩阵为A，那么A^m上的aij便表示i通过M个位置最后一个位置为J的方案数，然后再用快速幂优化即可。

PS：这里还要涉及到形如 a/b mod c的问题

处理方法要用到费马小定理，因为 b与c互质，则 b ^(c -1)= 1 (mod c)

所以 a /b mod c= a /b *（b ^ (c - 1)） mod c = a* (b ^ (c - 2)) mod c

这样就可以化除为乘。。甚是神奇

code：

  1 #include <iostream>
  2 #include <cstring>
  3 #include <cmath>
  4 #include <algorithm>
  5 #include <cstdio>
  6 #include <cstdlib>
  7 #define MXN 9973;
  8 using namespace std;
  9 
 10 int n , m, T, k;
 11 int b[50][10][10], prime[20000], bo[50001], v[10][10], u[10][10],  tot = 0, ans;
 12 
 13 void find_prime(){
 14      for (int i = 2; i <= 50000; ++i)
 15        if (!bo[i]){
 16            for (int j = i + i; j <= 50000; j += i)
 17              bo[j] = 1;
 18              prime[++tot] = i;
 19        }
 20 }
 21 
 22 void multi_m(int a[][10], int  b[][10] , int c[][10] ){
 23       for (int i = 0; i < m; ++i)
 24          for (int j = 0; j < m; ++j){
 25            c[i][j] = 0;
 26            for (int k = 0; k < m; ++k)
 27              c[i][j] =(c[i][j] + a[i][k]* b[k][j]) % MXN;
 28          }
 29              
 30     /* for (int i = 0; i < m; ++i){
 31        for (int j = 0; j < m; ++j)
 32           printf("%d  ", c[i][j]);
 33        puts("");
 34      }
 35      puts("");*/
 36 }
 37 
 38 void init(){
 39      int x, y;
 40      scanf("%d%d%d",&n, &m, &k);
 41      for (int i = 0; i < m; ++i)
 42        for (int j = 0; j < m; ++j)
 43            b[1][i][j] = 1;
 44      for (int i = 1; i <= k; ++i){
 45           scanf("%d%d", &x, &y);
 46           --x; --y;
 47           b[1][x][y] = b[1][y][x] =  0;
 48      }
 49      for (int i = 2; i <= 30; ++i)
 50          multi_m(b[i-1],  b[i-1],  b[i]);
 51          
 52   /*  for (int k = 1; k <= 30; ++k){
 53      for (int i = 0; i < m; ++i){
 54        for (int j = 0; j < m; ++j)
 55           printf("%d ", b[k][i][j]);
 56           puts("");
 57      }
 58     puts("\n");
 59    }*/
 60 }
 61 
 62 int phi(int now){
 63     int rea = now;
 64     for (int i = 1; prime[i] * prime[i] <= now; ++i)
 65        if (now % prime[i] == 0){
 66          rea -= rea / prime[i];
 67          while (now % prime[i] == 0) now /= prime[i];
 68        }
 69     if (now > 1)  rea -= rea / now;
 70     return rea % MXN;
 71 }
 72 
 73 int cal(int now){
 74     for (int i = 0; i < m; ++i)
 75        for (int j = 0; j < m; ++j)
 76           if (i == j) v[i][j] = 1;
 77              else v[i][j] = 0;
 78     int i = 1, ret = 0;
 79     
 80     while (now){
 81         if (now & 1) {
 82         
 83            memcpy(u, v, sizeof(v));   
 84            //memset(v, 0, sizeof(v));
 85            multi_m(u, b[i], v);   
 86         } 
 87         now >>= 1;
 88         ++i;
 89     }
 90     for (int i = 0; i < m; ++i)
 91           ret = (ret + v[i][i]) % MXN;
 92     return ret;    
 93 }
 94 
 95 void work(int now){
 96     int l = phi(now); 
 97     ans = (ans + l * cal(n / now)) %  MXN;
 98 }
 99 
100 void solve(){
101     ans = 0;
102     for (int i = 1; i * i <= n; ++i)
103         if (n % i == 0){
104               work(i);
105               if (i*i != n) work(n / i);
106         }
107     n %= MXN;
108     int num = MXN; 
109     for (int i = 1; i <= num - 2; ++i) ans = (ans * n) % MXN; 
110     printf("%d\n", ans);
111 }
112 
113 int main(){
114       freopen("poj2888.in", "r", stdin);
115       freopen("poj2888.out", "w", stdout);
116       scanf("%d", &T);
117       find_prime();
118       while (T--){
119           init();
120           solve();
121       }
122       fclose(stdin); fclose(stdout);
123 }