POJ 2888 Magic Bracelet （有限制的Burnside引理）

POJ 2154的加强,对相邻元素的颜色做出了限制，但是仍然用到了 循环个数 ： gcd(i,n)  循环长度：n/gcd(i,n)的结论，另外需要利用一个特性为对于含gcd(i,n)个长度为n/gcd(i,n)的置换,从置换的某一位置起的gcd（i,n）个元素是分别属于不同的循环的，同一个循环内的颜色是相同的，因此我们对于一个这样的置换求解的方案数便为：用m种颜色染色 a1,a2,...agcd(i,n)使得满足规定的条件，在这里我们用到了一个小技巧，利用邻接矩阵乘幂来求解方案数，具体方法见程序代码，总之做完这道题感觉还是很有收获的。。。

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
using namespace std;
const int MOD = 9973;
const int MAXN = 65540;
const int MAT_SIZE = 15;
bool is_prime[MAXN];
int prime[MAXN], len;
int T, N, M, K;

struct Matrix
{
    public:
        int n;
        int mat[MAT_SIZE][MAT_SIZE];
    public:

        void clear()
        {
            for(int i = 1; i <= n; ++i)
            {
                for(int j = 1; j <= n; ++j)
                {
                    mat[i][j] = 0;
                }
            }
        }

        void Init()
        {
            for(int i = 1; i <= n; ++i)
            {
                for(int j = 1; j <= n; ++j)
                {
                    if(i == j)
                        mat[i][j] = 1;
                    else
                        mat[i][j] = 0;
                }
            }
        }

        friend Matrix operator* (const Matrix &a, const Matrix &b)
        {
            Matrix c;
            c.n = a.n;
            c.clear();

            for(int i = 1; i <= a.n; ++i)
            {
                for(int j = 1; j <= a.n; ++j)
                {
                    for(int k = 1; k <= a.n; ++k)
                    {
                        if(a.mat[i][k] && b.mat[k][j])
                            c.mat[i][j] = (c.mat[i][j] + a.mat[i][k]*b.mat[k][j]);
                    }
                    c.mat[i][j] %= MOD;
                }
            }
            return c;
        }

        friend Matrix operator^ (Matrix a, int x)
        {
            Matrix c;
            c.n = a.n;
            c.Init();

            while(x)
            {
                if(x & 1)
                    c = c * a;
                a = a * a;
                x >>= 1;
            }
            return c;
        }

}A;


void Init()
{
    int i, j;
    len = 0;
    memset(is_prime, true, sizeof(is_prime));
    is_prime[0] = is_prime[1] = false;
    prime[len++] = 2;
    for(i = 4; i < MAXN; i += 2)
        is_prime[i] = false;
    for(i = 3; i * i <= MAXN; i += 2)
    {
        if(is_prime[i])
        {
            prime[len++] = i;
            for(j = i * i; j < MAXN; j += i)
                is_prime[j] = false;
        }
    }

    for( ; i < MAXN; ++i)
        if(is_prime[i])
            prime[len++] = i;
    return ;
}

int euler(int n)
{
    int ret = n;
    for(int i = 0; i < len && (long long)prime[i]*prime[i] <= n; ++i)
    {
        if(n % prime[i] == 0)
        {
            ret = ret / prime[i] * (prime[i] - 1);
            while(n % prime[i] == 0)
            {
                n /= prime[i];
            }
        }
    }
    if(n > 1)
        ret = ret / n * (n - 1);

    return ret % MOD;   //注意取余
}

int cal(int x)
{
    Matrix a; a.n = M;
    int ans = 0;
    a = A^x;
    for(int i = 1; i <= a.n; ++i)
    {
        ans = (ans + a.mat[i][i]) % MOD;
    }
    return ans;
}

int mul_mod(int a, int b)
{
    int ret = 0;
    while(b)
    {
        if(b & 1)
            ret = (ret + a) % MOD;
        a = (a + a) % MOD;
        b >>= 1;
    }
    return ret;
}

int power_mod(int a, int b)
{
    int ret = 1;
    while(b)
    {
        if(b & 1)
            ret = mul_mod(ret, a);
        a = mul_mod(a, a);
        b >>= 1;
    }
    return ret;
}

int main()
{
    //freopen("aa.in", "r", stdin);
    //freopen("bb.out", "w", stdout);

    int u, v;
    int ans;

    Init();
    scanf("%d", &T);
    while(T--)
    {
        scanf("%d %d %d", &N, &M, &K);
        ans = 0;
        A.n = M;
        for(int i = 1; i <= A.n; ++i)
        {
            for(int j = 1; j <= A.n; ++j)
            {
                A.mat[i][j] = 1;
            }
        }
        for(int i = 1; i <= K; ++i)
        {
            scanf("%d %d", &u, &v);
            A.mat[u][v] = A.mat[v][u] = 0;
        }

        for(int i = 1; i * i <= N; ++i)
        {
            if(N % i == 0)
            {
                if(i * i == N)
                {
                    ans = (ans + cal(i)*euler(i)) % MOD;
                }
                else
                {
                    ans = (ans + (cal(i)*euler(N/i)%MOD) + (cal(N/i)*euler(i))) % MOD;
                }
            }
        }
        ans = (ans * power_mod(N%MOD, MOD-2)) % MOD;
        cout << ans << endl;
    }
    return 0;
}