POJ 2888 Magic Bracelet [Polya 矩阵乘法]

传送门

题意：竟然扯到哈利波特了....

和上一题差不多，但颜色数很少，给出不能相邻的颜色对

---

可以相邻的连边建图矩阵乘法求回路个数就得到$f(i)$了....

感觉这样的环上有限制问题挺套路的...旋转的等价循环个数$t$我们很清楚了，并且环上每$t$个元素各属于不同的循环，我们只要求出$t$个元素满足限制的方案数就能得到$C(f)$了

然后再加上$gcd$取值讨论就降到根号了

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N=1e5+5,P=9973;
typedef long long ll;
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1; c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0'; c=getchar();}
    return x*f;
}
int n,m,k,u,v;
int p[N];
bool notp[N];
void sieve(int n){
    for(int i=2;i<=n;i++){
        if(!notp[i]) p[++p[0]]=i;
        for(int j=1;j<=p[0]&&i*p[j]<=n;j++){
            notp[i*p[j]]=1;
            if(i%p[j]==0) break;
        }
    }
}
inline int phi(int n){
    int re=n,m=sqrt(n);
    for(int i=1;i<=p[0]&&p[i]<=m&&p[i]<=n;i++) if(n%p[i]==0){
        re=re/p[i]*(p[i]-1);
        while(n%p[i]==0) n/=p[i];
    }
    if(n>1) re=re/n*(n-1);
    return re%P;
}
struct Matrix{
    int a[11][11];
    int* operator [](int x){return a[x];}
    Matrix(){memset(a,0,sizeof(a));}
    void ini(){for(int i=1;i<=10;i++) a[i][i]=1;}
}a;
Matrix operator *(Matrix a,Matrix b){
    Matrix c;
    for(int k=1;k<=m;k++)
        for(int i=1;i<=m;i++) if(a[i][k])
            for(int j=1;j<=m;j++) if(b[k][j])
                (c[i][j]+=a[i][k]*b[k][j])%=P;
    return c;
}
Matrix operator ^(Matrix a,int b){
    Matrix re;re.ini();
    for(;b;b>>=1,a=a*a)
        if(b&1) re=re*a;
    return re;
}
inline void mod(int &x){if(x>=P) x-=P;}
int f(int x){
    Matrix b=a^x;
    int re=0;
    for(int i=1;i<=m;i++) mod(re+=b[i][i]);
    return re;
}
inline int Pow(int a,int b){
    int re=1;
    a%=P;
    for(;b;b>>=1,a=a*a%P)
        if(b&1) re=re*a%P;
    return re;
}
inline int Inv(int a){return Pow(a,P-2);}
void solve(){
    int m=sqrt(n),ans=0;
    for(int i=1;i<=m;i++) if(n%i==0){
        mod(ans+= f(i)*phi(n/i)%P);
        if(i*i!=n) mod(ans+= f(n/i)*phi(i)%P);
    }
    printf("%d
",ans*Inv(n)%P);
}
int main(){
    freopen("in","r",stdin);
    sieve(32000);
    int T=read();
    while(T--){
        n=read();m=read();k=read();
        for(int i=1;i<=m;i++) for(int j=1;j<=m;j++) a[i][j]=1;
        for(int i=1;i<=k;i++){
            u=read();v=read();
            a[u][v]=a[v][u]=0;
        }
        solve();
    }
}