POJ 2888 Magic Bracelet

//
//  main.cpp
//  Richard
//
//  Created by 邵金杰 on 16/8/3.
//  Copyright © 2016年 邵金杰. All rights reserved.
//



#include<cstdio>
#include<cstring>
using namespace std;
const int maxn=50000;
int primes[maxn+1],is_prime[maxn+1];
const int mod=9973;
int n,m,k;
int u,v;
int cnt=0;
struct Matrix{
    int Mat[15][15];
};
Matrix operator * (Matrix a,Matrix b)
{
    Matrix c;
    for(int i=0;i<m;i++)
    {
        for(int j=0;j<m;j++)
        {
            c.Mat[i][j]=0;
            for(int k=0;k<m;k++)
            {
                c.Mat[i][j]=(c.Mat[i][j]+a.Mat[i][k]*b.Mat[k][j])%mod;
            }
        }
    }
    return c;
}
Matrix operator ^ (Matrix a,int b)
{
    Matrix temp;
    Matrix base=a;
    memset(temp.Mat,0,sizeof(temp.Mat));
    for(int i=0;i<m;i++) temp.Mat[i][i]=1;
    while(b)
    {
        if(b&1) temp=temp*base;
        base=base*base;
        b>>=1;
    }
    return temp;
}
void prime()
{
    memset(is_prime,1,sizeof(is_prime));
    for(int i=2;i<=maxn;i++)
    {
        if(is_prime[i]) primes[cnt++]=i;
        for(int j=0;j<i&&i*primes[j]<=maxn;j++)
        {
            is_prime[i*primes[j]]=0;
            if(i%primes[j]==0) break;
        }
    }
}
int solve(Matrix a,int b)
{
    int ret=0;
    Matrix temp=a^b;
    for(int i=0;i<m;i++)
        ret=(ret+temp.Mat[i][i])%mod;
    return ret;
}

int euler(int n)
{
    int ret=1;
    for(int i=0;i<cnt&&primes[i]*primes[i]<=n;i++)
    {
        if(n%primes[i]==0)
        {
            n/=primes[i];ret=(ret*(primes[i]-1));
            while(n%primes[i]==0) {n/=primes[i];ret=(ret*primes[i])%mod;}
        }
    }
    if(n>1) ret=(ret*(n-1));
    return ret%mod;
}
int Pow(int a,int b)
{
    int result=1;
    int base=a%mod;
    while(b)
    {
        if(b&1) result=(result*base)%mod;
        base=(base*base)%mod;
        b>>=1;
    }
    return result;
}
int polya(Matrix a)
{
    int result=0,i;
    for(i=1;i*i<n;i++)
    {
        if(n%i==0){
            result=(result+euler(i)*solve(a,n/i)+euler(n/i)*solve(a,i))%mod;
        }
    }
    if(i*i==n) result=(result+euler(i)*solve(a,i))%mod;
    return (result*Pow(n%mod,mod-2))%mod;
}
int main()
{
    prime();
    int t;
    scanf("%d",&t);
    while(t--)
    {
        Matrix a;
        scanf("%d%d%d",&n,&m,&k);
        for(int i=0;i<m;i++)
            for(int j=0;j<m;j++)
                a.Mat[i][j]=1;
        for(int i=0;i<k;i++){
            scanf("%d%d",&u,&v);
            u--,v--;
            a.Mat[u][v]=a.Mat[v][u]=0;
        }
        printf("%d\n",polya(a));
    }
    return 0;
}