POJ 2888

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#define MOD 9973
using namespace std;

struct Matrix{
	int m[12][12];
}mat[50];
int m;
bool isprime[35000];
int prime[35000],np;
int dive[100],dn;

Matrix operator * (Matrix a,Matrix b){
	Matrix ret;
	for(int i=1;i<=m;i++){
		for(int j=1;j<=m;j++){
			ret.m[i][j]=0;
			for(int k=1;k<=m;k++){
				ret.m[i][j]=(ret.m[i][j]+(a.m[i][k]*b.m[k][j])%MOD)%MOD;
			}
		}
	}
	return ret;
}
void getprime(){
	memset(isprime,true,sizeof(isprime));
	np=0;
	for(int i=2;i<35000;i++){
		if(isprime[i]){
			prime[np++]=i;
			for(int j=i*i;j<35000;j+=i){
				isprime[j]=false;
			}
		}
	}
}

void divn(int n){
	dn=0;
	int L =(int)sqrt(n*1.0);
	for(int i=1;i<=L;i++){
		if(n%i==0){
			dive[dn++]=i;
			if(i!=n/i)
			dive[dn++]=n/i;
		}
	}
}

void getinit(){
	for(int i=1;i<50;i++){
		mat[i]=mat[i-1]*mat[i-1];
	}
}

int phi(int p){
	int n=p;
	int res=p;
	for(int i=0;i<np&&prime[i]*prime[i]<=n;i++){
		if(p%prime[i]==0){
			res=res-res/prime[i];
			while(p%prime[i]==0)
			p/=prime[i];
		}
	}
	if(p>1)
	res=res-res/p;
	return res%MOD;
}

int quick(int b){
	 Matrix ans;
	 memset(ans.m,0,sizeof(ans.m));
	 for(int i=1;i<=m;i++)
	 ans.m[i][i]=1;
	 int k=0;
	 while(b){
	 	if(b&1) ans=ans*mat[k];
	 	b>>=1;
	 	k++;
	 }
	 int res=0;
	 for(int i=1;i<=m;i++)
	 res=(res+ans.m[i][i])%MOD;
	 return res;
}

void exgcd(int a,int b,int &x,int &y){
	if(b==0){
       x=1;
       y=0;
       return ;
    }
	exgcd(b,a%b,x,y);
    int t=x;
    x=y;
    y=t-a/b*y;
}

void slove(int n,int m){
	int x,y;
	exgcd(n,MOD,x,y);
	int ans=0;
	for(int i=0;i<dn;i++){
		ans=(ans+(phi(n/dive[i])%MOD)*quick(dive[i]))%MOD;
	}
	x=(x%MOD+MOD)%MOD;
	ans=(ans*x)%MOD;
	printf("%d\n",ans);
}

int main(){
	int n,k,p,q,T;
	getprime();
	scanf("%d",&T);
	while(T--){
		scanf("%d%d%d",&n,&m,&k);
		for(int i=1;i<=m;i++){
			for(int j=1;j<=m;j++)
			mat[0].m[i][j]=1;
		}
		for(int i=1;i<=k;i++){
			scanf("%d%d",&p,&q);
			mat[0].m[p][q]=mat[0].m[q][p]=0;
		}
		getinit();
		divn(n);
		slove(n,m);
	}
	return 0;
}

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long LL;
const int mr=100000;
const LL mod=9973;
bool notp[mr];
int pr[mr],fac[102],num[102];
int pn,top,n,m;
LL ans;
struct MAT
{
    LL bas[13][13];
    void init()
    {
        memset(bas,0,sizeof(bas));
    }
} mat[50];
MAT mul(MAT a,MAT b)
{
    MAT c;
    c.init();
    for(int i=1; i<=m; i++)
        for(int k=1; k<=m; k++)
        {
            if(a.bas[i][k])
            {
                for(int j=1; j<=m; j++)
                {
                    c.bas[i][j]+=a.bas[i][k]*b.bas[k][j];
                    if(c.bas[i][j]>=mod)
                        c.bas[i][j]%=mod;
                }
            }
        }
    return c;
}
void getpri()//筛素数
{
    pn=0;
    memset(notp,0,sizeof(notp));
    for(int i=2; i<mr; i++)
    {
        if(!notp[i])
        {
            pr[pn++]=i;
        }
        for(int j=0; j<pn && i*pr[j]<mr; j++)
        {
            int k=i*pr[j];
            notp[k]=1;
            if(i%pr[j]==0)break;
        }
    }
}
void divn()
{
    int nn=n;
    top=0;
    int lim=(int)sqrt((double(nn)))+1;
    for(int i=0; pr[i]<=lim; i++)
    {
        if(nn%pr[i]==0)
        {
            fac[top]=pr[i];
            num[top]=0;
            while(nn%pr[i]==0)
                num[top]++,nn/=pr[i];
            top++;
        }
    }
    if(nn>1)
        fac[top]=nn,num[top++]=1;
}
int phi(int x)
{
    int i, res=x;
    for (i=0;pr[i]<(int)sqrt((double)x)+1;i++)
    if(x%pr[i]==0)
    {
        res=res/pr[i]*(pr[i]-1);
        while(x%pr[i]==0)x/=pr[i];
    }
    if(x>1)res=res/x*(x-1);
    return res;
}
void solve(int r)
{
    int res=phi(n/r);
    MAT mt;
    mt.init();
    for(int i=1;i<=m;i++)
        mt.bas[i][i]=1;
    for(int i=1,tp=r;tp;i++,tp>>=1)
    if(tp&1)mt=mul(mt,mat[i]);
    for(int i=1;i<=m;i++)
    {
        ans+=mt.bas[i][i]*res;
        if(ans>=mod)ans%=mod;
    }
}
void dfs(int id,int sum)
{
    if(id==top)
    {
        solve(sum);
        return;
    }
    else
    {
        dfs(id+1,sum);
        for(int ct=0; ct<num[id]; ct++)
            dfs(id+1,sum=sum*fac[id]);
    }
}
void init()
{
    for(int i=2; i<50; i++)
        mat[i]=mul(mat[i-1],mat[i-1]);
}
int Egcd (int a,int b, int &x, int &y)
{
    if (b==0)
    {
        x=1,y=0;
        return a;
    }
    LL d, tp;
    d = Egcd (b, a%b, x, y);
    tp = x;
    x = y;
    y = tp - a/b*y;
    return d;
}
int getni()
{
    int x,y;
    Egcd(n,mod,x,y);
    return (x%mod+mod)%mod;
}
int main()
{
    getpri();
    int T;
    for(scanf("%d",&T); T; T--)
    {
        int k;
        scanf("%d%d%d",&n,&m,&k);
        ans=0;
        for(int i=1; i<=m; i++)
            for(int j=1; j<=m; j++)
                mat[1].bas[i][j]=1;
        for(int a,b,i=0; i<k; i++)
        {
            scanf("%d%d",&a,&b);
            mat[1].bas[a][b]=mat[1].bas[b][a]=0;
        }
        init();
        divn();
        dfs(0,1);
        printf("%d\n",ans*getni()%mod);
    }
    return 0;
}