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：

#include <iostream>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#define MXN 9973;
using namespace std;

int n , m, T, k;
int b[50][10][10], prime[20000], bo[50001], v[10][10], u[10][10],  tot = 0, ans;

void find_prime(){
     for (int i = 2; i <= 50000; ++i)
       if (!bo[i]){
           for (int j = i + i; j <= 50000; j += i)
             bo[j] = 1;
             prime[++tot] = i;
       }
}

void multi_m(int a[][10], int  b[][10] , int c[][10] ){
      for (int i = 0; i < m; ++i)
         for (int j = 0; j < m; ++j){
           c[i][j] = 0;
           for (int k = 0; k < m; ++k)
             c[i][j] =(c[i][j] + a[i][k]* b[k][j]) % MXN;
         }
             
    /* for (int i = 0; i < m; ++i){
       for (int j = 0; j < m; ++j)
          printf("%d  ", c[i][j]);
       puts("");
     }
     puts("");*/
}

void init(){
     int x, y;
     scanf("%d%d%d",&n, &m, &k);
     for (int i = 0; i < m; ++i)
       for (int j = 0; j < m; ++j)
           b[1][i][j] = 1;
     for (int i = 1; i <= k; ++i){
          scanf("%d%d", &x, &y);
          --x; --y;
          b[1][x][y] = b[1][y][x] =  0;
     }
     for (int i = 2; i <= 30; ++i)
         multi_m(b[i-1],  b[i-1],  b[i]);
         
  /*  for (int k = 1; k <= 30; ++k){
     for (int i = 0; i < m; ++i){
       for (int j = 0; j < m; ++j)
          printf("%d ", b[k][i][j]);
          puts("");
     }
    puts("\n");
   }*/
}

int phi(int now){
    int rea = now;
    for (int i = 1; prime[i] * prime[i] <= now; ++i)
       if (now % prime[i] == 0){
         rea -= rea / prime[i];
         while (now % prime[i] == 0) now /= prime[i];
       }
    if (now > 1)  rea -= rea / now;
    return rea % MXN;
}

int cal(int now){
    for (int i = 0; i < m; ++i)
       for (int j = 0; j < m; ++j)
          if (i == j) v[i][j] = 1;
             else v[i][j] = 0;
    int i = 1, ret = 0;
    
    while (now){
        if (now & 1) {
        
           memcpy(u, v, sizeof(v));   
           //memset(v, 0, sizeof(v));
           multi_m(u, b[i], v);   
        } 
        now >>= 1;
        ++i;
    }
    for (int i = 0; i < m; ++i)
          ret = (ret + v[i][i]) % MXN;
    return ret;    
}

void work(int now){
    int l = phi(now); 
    ans = (ans + l * cal(n / now)) %  MXN;
}

void solve(){
    ans = 0;
    for (int i = 1; i * i <= n; ++i)
        if (n % i == 0){
              work(i);
              if (i*i != n) work(n / i);
        }
    n %= MXN;
    int num = MXN; 
    for (int i = 1; i <= num - 2; ++i) ans = (ans * n) % MXN; 
    printf("%d\n", ans);
}

int main(){
      freopen("poj2888.in", "r", stdin);
      freopen("poj2888.out", "w", stdout);
      scanf("%d", &T);
      find_prime();
      while (T--){
          init();
          solve();
      }
      fclose(stdin); fclose(stdout);
}