POJ-2888 Magic Bracelet（Burnside引理+矩阵优化+欧拉函数+逆元）

Burnside引理经典好题呀！

题解参考 https://blog.csdn.net/maxwei_wzj/article/details/73024349#commentBox 这位大佬的。

这题时间卡得很紧，注意矩阵乘法不能太多次取模，不然会TLE。   因为这题的模数是9973，加完之后再取模也不会炸。

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
using namespace std;
const int N=1e6+10;
const int MOD=9973;
int n,m,k,tot=0;

struct matrix{
    int m[12][12];
    matrix() { memset(m,0,sizeof(m)); }
    friend matrix operator*(matrix a,matrix b) {
        matrix res;
        for (int i=1;i<=10;i++) for (int j=1;j<=10;j++) {
            for (int k=1;k<=10;k++) res.m[i][j]=res.m[i][j]+a.m[i][k]*b.m[k][j];
            res.m[i][j]%=MOD;  //少取模不然会TLE 
        }
        return res;    
    }
};

bool vis[N]; int pri[N];
void prework(int n) {
    memset(vis,false,sizeof(vis));  //vis为0素数  1和数 
    for (int i=2;i<=n;i++){
        if (!vis[i]) pri[++tot]=i;
        for (int j=1;j<=tot&&i*pri[j]<=n;j++){
            vis[i*pri[j]]=1;
            if (i%pri[j]==0) break;
        }
    }
}

int power(int x,int p) {
    int ret=1; x%=MOD;
    for (;p;p>>=1) {
        if (p&1) ret=(ret*x)%MOD;
        x=(x*x)%MOD;
    }
    return ret;
}

int phi(int n) {
    int ret=n;
    for (int i=1;i<=tot && pri[i]<=sqrt((double)n);i++)
        if (n%pri[i]==0) {
            ret=ret/pri[i]*(pri[i]-1);
            while (n%pri[i]==0) n/=pri[i];
        }
    if (n>1) ret=ret/n*(n-1);
    return ret%MOD;    
}

int count(matrix A,int p) {  //计算B=A^p  然后返回sigma(B[i][i]) 
    matrix B;
    for (int i=1;i<=m;i++) B.m[i][i]=1;
    for (;p;p>>=1) {
        if (p&1) B=B*A;
        A=A*A;
    }
    int ret=0;
    for (int i=1;i<=m;i++) ret=(ret+B.m[i][i])%MOD;
    return ret;
}

//POJ-2888 由n(n <= 10^9)个珠子组成的项链，每个珠子共有m(m <= 10)种颜色，再给定k组限制(a, b)表示颜色a和颜色b的珠子不能相邻，问总共有多少种方案满足长度为n的项链。
//题解：Burnside引理 + 矩阵优化 + 欧拉函数 + 逆元。
int main()
{
    prework(1000000);
    int T; cin>>T;
    while (T--) {
        scanf("%d%d%d",&n,&m,&k);
        
        matrix A;  //可达矩阵（类似图论） 
        for (int i=1;i<=k;i++) {
            int x,y; scanf("%d%d",&x,&y);
            A.m[x][y]=A.m[y][x]=1;
        }
        for (int i=1;i<=m;i++) for (int j=1;j<=m;j++) A.m[i][j]=1-A.m[i][j];
        
        int ans=0;
        for (int i=1;i*i<=n;i++)
            if (n%i==0) {
                ans=(ans+count(A,i)*phi(n/i)%MOD)%MOD;
                if (i*i!=n) ans=(ans+count(A,n/i)*phi(i)%MOD)%MOD;
            }
        cout<<ans*power(n,MOD-2)%MOD<<endl;    
    }
    return 0;
}