【题解】POJ 2888 魔法手镯[Burnside+莫反+矩乘]
CSDN博客

#include<cstdio>
#include<bits/stdc++.h>
using namespace std;
template<typename T> void qread(T &x){
    x=0; int f=1; char c=getchar();
    for(; !isdigit(c); c=getchar()) if(c=='-') f=-1;
    for(; isdigit(c); c=getchar()) x=x*10+(c-'0');
    x*=f;
}
typedef long long LL;
const int N=5e6+10, M=15;
const LL P=9973;
struct matx{
    LL a[M][M];
    matx(){
        for(int i=1; i<=10; i++)
            for(int j=1; j<=10; j++)
                a[i][j]=0;
    }
} ;
int m;
matx operator*(matx A, matx B){
    matx C;
    for(int i=1; i<=m; i++)
        for(int j=1; j<=m; j++)
            for(int k=1; k<=m; k++)
                C.a[i][j]=(C.a[i][j]+A.a[i][k]*B.a[k][j])%P;
    return C;
}
matx mq_pow(matx a, LL b){
    matx res;
    for(int i=1; i<=m; i++) res.a[i][i]=1;
    while(b){
        if(b&1) res=res*a;
        a=a*a; b/=2;
    }
    return res;
}
LL calc(LL x){
    LL res=x;
    for(LL i=2; i<=sqrt(x+1); i++) if(x%i==0){
        res=res/i*(i-1);
        while(x%i==0) x/=i;
    }
    if(x>1) res=res/x*(x-1);
    return res;
}
matx mt;
LL work(LL d){
    LL res=0;
    matx sum=mq_pow(mt, d);
    for(int i=1; i<=m; i++)
        res=(res+sum.a[i][i])%P;

【题解】POJ 2888 魔法手镯[Burnside+莫反+矩乘]-CSDN博客

#include<bits/stdc++.h>
using namespace std;
template<typename T> void qread(T &x){
    x=0; int f=1; char c=getchar();
    for(; !isdigit(c); c=getchar()) if(c=='-') f=-1;
    for(; isdigit(c); c=getchar()) x=x*10+(c-'0');
    x*=f;
}
typedef long long LL;
const int N=5e6+10, M=15;
const LL P=9973;
struct matx{
    LL a[M][M];
    matx(){
        for(int i=1; i<=10; i++)
            for(int j=1; j<=10; j++)
                a[i][j]=0;
    }
} ;
int m;
matx operator*(matx A, matx B){
    matx C;
    for(int i=1; i<=m; i++)
        for(int j=1; j<=m; j++)
            for(int k=1; k<=m; k++)
                C.a[i][j]=(C.a[i][j]+A.a[i][k]*B.a[k][j])%P;
    return C;
}
matx mq_pow(matx a, LL b){
    matx res;
    for(int i=1; i<=m; i++) res.a[i][i]=1;
    while(b){
        if(b&1) res=res*a;
        a=a*a; b/=2;
    }
    return res;
}
LL calc(LL x){
    LL res=x;
    for(LL i=2; i<=sqrt(x+1); i++) if(x%i==0){
        res=res/i*(i-1);
        while(x%i==0) x/=i;
    }
    if(x>1) res=res/x*(x-1);
    return res;
}
matx mt;
LL work(LL d){
    LL res=0;
    matx sum=mq_pow(mt, d);
    for(int i=1; i<=m; i++)
        res=(res+sum.a[i][i])%P;
    return res;
}
LL q_pow(LL a, LL b){
    LL res=1;
    while(b){
        if(b&1) res=res*a%P;
        a=a*a%P; b/=2;
    }
    return res;
}
int main(){
    int T; qread(T);
    while(T--){
	    for(int i=1; i<=10; i++)
	        for(int j=1; j<=10; j++)
	            mt.a[i][j]=1;
        LL n; int K; 
        qread(n); qread(m); qread(K);
    LL inv=q_pow(n, P-2);

    for(int i=1; i&lt;=K; i++){
        int a, b; 
        qread(a); qread(b);
        mt.a[a][b]=0;
        mt.a[b][a]=0;
    }
    if(n==1){
		printf("%d\n", m);
		continue;
	}

    LL ans=0;
    for(LL d=1; d&lt;=sqrt(n+1); d++) if(n%d==0){
		ans=(ans+calc(n/d)*work(d)%P)%P;
        if(d*d!=n) ans=(ans+calc(d)*work(n/d)%P)%P;
    }
    printf("%lld\n", ans*inv%P);
}
return 0;

}