Description
用(m(mleq10))种颜色给一个(n(nleq10^9))个点的环染色，其中有(k)对颜色不能相邻。求在旋转同构下本质不同的方案数，对(9973)取模。
Solution
设(f(x))表示用(m)种颜色给一个(x)个点的环染色的在满足约束下的方案数。那么
[egin{align*}
ncdot ans&=sum_{i=1}^n f(gcd(i,n)) \
&= sum_{d|n} f(d) sum_{d|i}^n [gcd(i,n)=d] \
&= sum_{d|n} f(d) sum_{i=1}^{frac{n}{d}} [gcd(i,frac{n}{d})=1] \
&= sum_{d|n} f(d) varphi(frac{n}{d})
end{align*}$$ $varphi(x)$可以在对$n$分解质因数后简单求解，那么接下来只需要考虑如何求$f(x)$。
设$dp[i][j]$表示一个长度为$i$的序列以颜色$j$结尾时的合法方案数。构造转移矩阵$M$，满足颜色$i,j$不能相邻时$M_{i,j}=0$，其余的$M_{i,j}=1$。那么有$dp[i] imes M=dp[i+1]$，$dp[x]=dp[1] imes M^{x-1}$。
$f(x)$等于初始颜色为$cin[1,m]$，即$dp[1][c]=1$时$dp[x+1][c]$的和。而$dp[x+1][c]=sum_{i=1}^m dp[1][i] imes M^x[i][c]=M^x[c][c]$，所以$f(x)=sum_{c=1}^m M^x[c][c]$。
> 时间复杂度$O(sigma_0(n)m^3logn)$。
##Code
```cpp
//Magic Bracelet
#include <cstdio>
#include <cstring>
inline char gc()
{
static char now[1<<16],*s,*t;
if(s==t) {t=(s=now)+fread(now,1,1<<16,stdin); if(s==t) return EOF;}
return *s++;
}
inline int read()
{
int x=0; char ch=gc();
while(ch<'0'||'9'<ch) ch=gc();
while('0'<=ch&&ch<='9') x=x*10+ch-'0',ch=gc();
return x;
}
const int N=2e6+10;
int prCnt,pr[N]; bool prNot[N];
void init(int n)
{
for(int i=2;i<=n;i++)
{
if(!prNot[i]) pr[++prCnt]=i;
for(int j=1;j<=prCnt;j++)
{
if(i*pr[j]>n) break;
prNot[i*pr[j]]=true;
if(i%pr[j]==0) break;
}
}
}
const int P=9973;
int n,m,k;
int dCnt,dCnt1,d[2000],d1[50];
struct mtx
{
int c,r,x[20][20];
mtx(int _r,int _c) {r=_r,c=_c,memset(x,0,sizeof x);}
};
mtx operator *(mtx A,mtx B)
{
mtx C=mtx(A.r,B.c);
for(int i=1;i<=A.r;i++)
for(int k=1;k<=A.c;k++)
for(int j=1;j<=B.c;j++)
C.x[i][j]=(C.x[i][j]+A.x[i][k]*B.x[k][j])%P;
return C;
}
mtx I(int n)
{
mtx A=mtx(n,n);
for(int i=1;i<=n;i++) A.x[i][i]=1;
return A;
}
mtx pow(mtx A,int y)
{
mtx r=I(A.r),t=A;
for(int i=y;i;i>>=1,t=t*t) if(i&1) r=r*t;
return r;
}
int pow(int x,int y)
{
int r=1,t=x%P;
for(int i=y;i;i>>=1,t=t*t%P) if(i&1) r=r*t%P;
return r;
}
int phi(int x)
{
int r=x;
for(int i=1;i<=dCnt1;i++)
if(x%d1[i]==0) r=r/d1[i]*(d1[i]-1);
return r%P;
}
int main()
{
int task=read(); init(2e6);
while(task--)
{
n=read(),m=read(),k=read();
dCnt=0;
for(int i=1;i*i<=n;i++) if(n%i==0) d[++dCnt]=i;
int isSqr=(d[dCnt]*d[dCnt]==n);
for(int i=1;i<=dCnt;i++) d[dCnt*2+(isSqr^1)-i]=n/d[i];
dCnt+=dCnt-isSqr;
dCnt1=0;
for(int i=1,t=n;t>1;i++)
{
if(pr[i]*pr[i]>t) {d1[++dCnt1]=t; break;}
if(t%pr[i]==0) d1[++dCnt1]=pr[i];
while(t%pr[i]==0) t/=pr[i];
}
mtx M=mtx(m,m),Mk=mtx(m,m);
for(int i=1;i<=m;i++)
for(int j=1;j<=m;j++) M.x[i][j]=1;
for(int i=1;i<=k;i++)
{
int c1=read(),c2=read();
M.x[c1][c2]=M.x[c2][c1]=0;
}
int ans=0;
for(int i=1;i<=dCnt;i++)
{
int f=0; Mk=pow(M,d[i]);
for(int j=1;j<=m;j++) f=(f+Mk.x[j][j])%P;
ans+=f*phi(n/d[i])%P; while(ans>=P) ans-=P;
}
printf("%d
",ans*pow(n,P-2)%P);
}
return 0;
}
```]