思路挺清晰的。不过，我就是WA。不清楚为什么，很多数据都过了。
其实，一个置换后若有循环节个数为K，则N必定可以除以尽K。而K正好可以看成一个环。为什么呢？看前K个珠子，就是一个环，而后面的若干个K个珠子，不过就是不停的重复而已。这样，循环节的个数可以由最大公约数求得。那么，这个K个珠子构成的环符合题意的有多少种呢？很巧妙的一个方法是，用矩阵表示，若颜色相邻则I,J可以为1，否则为0。矩阵相乘有一个应用就是求的路径数啊。
最后，求逆元即可。可我的就是不过，求大神路过时指点。
#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
",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;
}
以下是别人的代码：http://blog.csdn.net/tmeteorj/article/details/8654330
#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
",ans*getni()%mod);
}
return 0;
}