POJ 2888 Magic Bracelet （Burnside + 矩阵加速dp）

　　额，这个题也做了好几天了，中间停停断断的……

　　一开始写是超时，想到了正确的思路，用dp[i][j]代表到i个珠子颜色为j的方案数，但是超时了，因为我细节处理的太差，我先枚举了初始状态，然后又枚举了结束状态，多了两个m，因为我想的是判断最后一个和第一个的关系，但其实完全可以让矩阵帮助我多做一次运算，找对应的值就可以，所以结束状态是不需要枚举的，开始状态可以把原先1*N的初始矩阵，改成N*N的矩阵，这里正好是单位阵，每一行都代表一个独立的初始状态，矩阵在计算的时候，就会把答案存储到自身对应的行，最终只要把对角线上的求和即可。这样是很聪明的做法，因为这个题目数据太强，差一点都不行，但是我又不敢减少取模，因为怕溢出。

　　但是今晚改了之后，发现WA了，我一直都以为n的因子不会超过两百个，但是我自己生成数据拍出了288个因子，干脆改成2000过去了，后来发现因子会多于1000个，这也是很符合道理的，算是纠正了我一个误区吧。以后尽量有vector存储吧。但是有意思的是，这个题其实不需要存……是我多此一举了

　　代码：

　　

#include <iostream>
#include <cstdio>
#include <cstring>
#include <vector>
using namespace std;
#define LL long long
const int Mod = 9973;
int n,m,k,ok[12][12];
struct Mar{
    int n,m,a[12][12];
    Mar(int nn,int mm){
        n=nn; m=mm;
        memset(a,0,sizeof(a));
    }
    Mar operator*(const Mar&B)const{
        Mar C(n,B.m);
        for(int i=1;i<=n;i++){
            for(int j=1;j<=m;j++){
                if(a[i][j]){
                    for(int k=1;k<=B.m;k++){
                        C.a[i][k]=(C.a[i][k]+a[i][j]*B.a[j][k]%Mod)%Mod;
                    }

                }
            }
        }
        return C;
    }
};
Mar MarPow(Mar x,int y){
    Mar res(x.n,x.m);
    for(int i=1;i<=x.n;i++){
        res.a[i][i] = 1;
    }
    while(y){
        if(y&1){
           res = res*x;
        }
        x = x*x;
        y>>=1;
    }
    return res;
}
//#include<cmath>
int yz[2000],lyz;
int p[4600],pct;
bool isp[32622];
void init(){
    memset(isp,0,sizeof(isp));
    pct=0; int up = 32600;
    for(int i=2;i<=up;i++){
        if(!isp[i]){
            p[pct++] = i;
        }
        for(int j=0;j<pct;j++){
            LL x = 1LL*p[j]*i;
            if(x > 1LL*up) break;
            isp[x] = 1;
            if(i%p[j]==0) break;
        }
    }
//    cout<<pct<<endl;
}
int Euler(LL x){
    LL res = x;
    for(int i=0;i<pct;i++){
        if(1LL*p[i]*p[i] > x) break;
        if(x == 1) break;
        if(x%p[i]==0){
            res = res/p[i]*(p[i]-1);
            x /= p[i];
            while(x%p[i]==0){
                x/=p[i];
            }
        }
    }
    if(x!=1) res = res/x*(x-1);
    return res%Mod;
}
int ex_gcd(int a,int b,int &x,int &y){
    if(b==0){
        x = 1;  y = 0;
        return a;
    }
    int g = ex_gcd(b,a%b,y,x);
    y -= a/b*x;
    return g;
}
int read(){
    int res = 0; char c;
    while((c=getchar())){
        if(c>='0'&&c<='9'){
            res = res+(c-'0');
            break;
        }
    }
    while((c=getchar())){
        if(c>='0'&&c<='9'){
            res *= 10;
            res = res+(c-'0');
        }else break;
    }
    return res;
}
int main(){
//    freopen("in.cpp","r",stdin);
//    freopen("out1","w",stdout);
    int T;
    init();
    T = read();
    while(T--){
        n = read(); m = read(); k = read();
        lyz = 0;
        for(int i=1;1LL*i*i<=(LL)(n);i++){
            if(n%i==0){
                yz[lyz++]=i;
                int j = n/i;
                if(i!=j) yz[lyz++]=j;
            }
        }
        memset(ok,0,sizeof(ok));
        while(k--){
            int u,v;
            u = read(); v = read();
            ok[u][v] = ok[v][u] = 1;
        }
        Mar B(m,m);
        for(int i=1;i<=m;i++){
            for(int j=1;j<=m;j++){
                if(!ok[i][j]) {
                    B.a[j][i] = 1;
                }
            }
        }
        int ans = 0;
        for(int i=0;i<lyz;i++){
            int all = yz[i]; int sum = Euler(n/yz[i]);
            Mar Now = MarPow(B,all);
            int tmp = 0;
            for(int j=1;j<=m;j++){
                tmp = (tmp+Now.a[j][j])%Mod;
            }
            ans = (ans+(tmp*sum)%Mod)%Mod;
        }
//        printf("%d
",ans);
//        printf("ans1 = %d
",ans);
//         printf("now ans = "); cout<<ans/n<<endl;
        int xx,yy;
        ex_gcd(n,Mod,xx,yy);
        xx = (xx%Mod+Mod)%Mod;
//        xx = qpow(n,Mod-2);
        ans = (ans*xx)%Mod;
        printf("%d
",ans);
    }
    return 0;
}
/**
5

2 4 1
1 2

1000000000 10 0

1000000000 10 3
1 2
3 4
5 5

9972 10 1
1 10

1234 1 0
*/