POJ - 2888 Magic Bracelet Ploya+矩阵

题目：有n个珠子的环，有m种颜色，旋转视为相同。但是有一些颜色限制，规定某两种种颜色不能相邻。

思路：运用矩阵，a->b->c->a，表示一个长度为3的一种方案，而且最终是回到起点的。

这样颜色限制解决了，然后用Burnside引理去做。

注意：要尽量少取模，不然会超时。

代码：

#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<iostream>
#include<algorithm>
#include<ctime>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<string>
#include<vector>
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<list>
#include<numeric>
using namespace std;
#define LL long long
#define ULL unsigned long long
#define INF 0x3f3f3f3f3f3f3f3f
#define mm(a,b) memset(a,b,sizeof(a))
#define PP puts("*********************");
template<class T> T f_abs(T a){ return a > 0 ? a : -a; }
template<class T> T gcd(T a, T b){ return b ? gcd(b, a%b) : a; }
template<class T> T lcm(T a,T b){return a/gcd(a,b)*b;}
// 0x3f3f3f3f3f3f3f3f

const int MOD=9973;
const int maxn=4e4+50;
bool isprime[maxn];
int prime[maxn],tol;
void make_prime(int n){
    for(int i=0;i<=n;i++)
        isprime[i]=true;
    isprime[0]=isprime[1]=false;
    tol=0;
    for(int i=2;i<=n;i++){
        if(isprime[i])
            prime[tol++]=i;
        for(int j=0;j<tol;j++){
            if(i*prime[j]<=n)
                isprime[i*prime[j]]=false;
            else
                break;
            if(i%prime[j]==0)
                break;
        }
    }
}
const int MAXN = 13;
int sz;
struct Matrix {
	int a[MAXN][MAXN];
	void toZero() {
		for(int i = 0; i < sz; ++i)
			for(int j = 0; j < sz; ++j)
				a[i][j] = 0;
	}
	void toOne() {
		for(int i = 0; i < sz; ++i) {
			for(int j = 0; j < sz; ++j)
				a[i][j] = (int)(i == j);
		}
	}
	Matrix operator + (const Matrix &B) const {
		Matrix A(*this);
		int i, j;
		for(i = 0; i < sz; ++i) {
			for(j = 0; j < sz; ++j)
				//A.a[i][j] += B.a[i][j];
				A.a[i][j] = (A.a[i][j] + B.a[i][j]) % MOD;
		} return A;
	}
	Matrix operator * (const Matrix &B) const {
		Matrix A(*this), res;
		res.toZero();
		int i, j, k;
		for(i = 0; i < sz; ++i) {
			for(j = 0; j < sz; ++j) {
				if(A.a[i][j] == 0) continue;
				for(k = 0; k < sz; ++k){
					//res.a[i][k] += A.a[i][j] * B.a[j][k];
					res.a[i][k] = (res.a[i][k] + (A.a[i][j] * B.a[j][k]));
					res.a[i][k]%=MOD;
				}
			}
		} return res;
	}
	Matrix operator ^ (int p) const {
		Matrix res, X(*this);
		res.toOne();
		for(; p; p >>= 1) {
			if(p & 1) res = res * X;
			X = X * X;
		} return res;
	}
	void Disp() {
		printf("---------------------\n");
		for(int i = 0; i < sz; ++i) {
			for(int j = 0; j < sz; ++j)
				printf("%4d", a[i][j]);
			printf("\n");
		}
		printf("---------------------\n");
	}
}A;
int phi(int n){
    int res=n;
    for(int i=0;i<tol&&prime[i]*prime[i]<=n;i++)
        if(n%prime[i]==0){
            res=res-res/prime[i];
            while(n%prime[i]==0){
                n/=prime[i];
            }
        }
    if(n>1)
        res=res-res/n;
    return res%MOD;
}
int solve(int k){
    Matrix B=A^k;
    int sum=0;
    for(int i=0;i<sz;i++)
        sum=(sum+B.a[i][i])%MOD;
    return sum;
}
int pow_mod(int a,int b){
    int res=1;
    a=a%MOD;
    while(b){
        if(b&1) res=res*a%MOD;
        a=a*a%MOD;
        b/=2;
    }
    return res;
}
int main(){

    int T,n,m,k,u,v;
    make_prime(maxn-50);
    scanf("%d",&T);
    while(T--){
        scanf("%d%d%d",&n,&m,&k);
        sz=m;
        for(int i=0;i<sz;i++)
            for(int j=0;j<sz;j++)
                A.a[i][j]=1;
        while(k--){
            scanf("%d%d",&u,&v);u--;v--;
            A.a[u][v]=A.a[v][u]=0;
        }
        int ans=0,i;
        for(i=1;i*i<n;i++){
            if(n%i)
                continue;
            ans=(ans+phi(i)*solve(n/i)%MOD)%MOD;
            ans=(ans+phi(n/i)*solve(i)%MOD)%MOD;
        }
        if(i*i==n)
            ans=(ans+phi(i)*solve(i)%MOD)%MOD;
        ans=ans*pow_mod(n,MOD-2)%MOD;
        printf("%d\n",ans);
    }
    return 0;
}