【POJ2888】Magic Bracelet-Burnside引理+数论+DP矩阵优化

测试地址：Magic Bracelet
题目大意：要用N(≤109)颗珠子连成环形的手镯，共有M(≤10)种不同的珠子，规定K个条件，每个条件规定某两种珠子不能相邻，旋转后相同的方案视作相同，问有多少种本质不同的方案，对9973取模（gcd(N,9973)=1）。
做法：这题非常经典，思路很有代表性，是进阶Burnside引理和Polya定理题的一块敲门砖，需要用到：Burnside引理，欧拉函数，DP+矩阵快速幂，乘法逆元。
首先一看题目是计数，就自然而然地联想到Burnside引理和Polya定理，然而这题有不能相邻的条件，所以不能直接用Polya定理，那么我们考虑使用Burnside引理来解决。
我们知道环形的旋转置换有N种，第i种置换（这里定义为旋转i颗珠子的置换）的循环数为gcd(N,i)，观察这样的置换，我们发现它很有规律：第1个元素属于第1个循环，第2个元素属于第2个循环，…，第gcd(N,i)个元素属于第gcd(N,i)个循环，第gcd(N,i)+1个元素属于第1个循环，以此类推。由于我们计算|C(f)|时每个循环内元素着色应该相同，那么我们实际上只需确定前gcd(N,i)个元素的填色方案数即可。设对前i个元素进行着色，且最后一个元素颜色为j的方案数为f(i,j)，并设二元函数m(i,j)，表示如果i,j两种颜色的珠子可以相邻，那么m(i,j)=1，否则m(i,j)=0，那么我们可以得到状态转移方程：
f(i,j)=∑Mk=1m(j,k)×f(i−1,k)
但是最后一个元素的填色除了和它前一个元素有关，还和第一个元素有关，那么这个状态转移方程是不是错的呢？实际上，我们只需一开始枚举第一个元素的颜色k，然后这一种情况的答案就是第gcd(N,i)+1个元素与第一个元素同色的方案数，即：f(gcd(N,i)+1,k)，累加每种情况即可得出|C(f)|。带进Burnside公式计算之后，最外面还要乘一个1/N，由于gcd(N,9973)=1，所以我们直接求N关于模9973的乘法逆元，然后再乘即可。
上面的方法的时间复杂度是O(N2M)，显然不能通过此题，需要考虑优化。此题需要从两个方面入手，一是计算Burnside公式的复杂度，而是动态规划的复杂度。
首先优化计算Burnside公式的时间复杂度。我们知道gcd(N,i)|N，而且当gcd(N,i)相同时，|C(f)|相同，那么我们完全可以枚举N的因子d，计算gcd(N,i)=d时的|C(f)|，再乘上使得gcd(N,i)=d的i的个数，累加起来得到答案。后面那个个数显然就是φ(N/d)了（不懂的可以去看看我写的POJ2154的题解），那么枚举的时间复杂度就大大降低了。
接下来优化DP的时间复杂度，注意到M很小，那么二元函数m显然可以构造成一个矩阵M′，如果我们把f(i,1),f(i,2),...,f(i,M)从上到下排成一个列向量Fi，那么可以看出：Fi=M′Fi−1，所以Fi=M′i−1F1，于是我们可以先用矩阵快速幂算出M′d，然后就可以算Fd+1了，但其实我们没有必要算出全部的Fd+1，因为我们只要求f(d+1,k)，这是第k行第1列的元素，那么就用M′d的第k行与F1的第1列做运算，根据定义，F1中只有第k行为1，其余都为0，所以f(d+1,k)=M′dkk，那么|C(f)|就应该等于M′d对角线上的元素和。
经过了以上两个优化，外层枚举的复杂度降到差不多O(N−−√)，DP的时间复杂度从O(NM)降到O(M3logN)，因此总的时间复杂度为O(M3N−−√logN)，可以通过此题。
以下是本人代码：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define ll long long
using namespace std;
ll n,ans;
int T,m,k;
struct matrix {ll s[11][11];} M[40];

void exgcd(ll a,ll b,ll &x,ll &y)
{
  ll x0=1,y0=0,x1=0,y1=1;
  while(b)
  {
    ll tmp,q;
    q=a/b;
    tmp=x0,x0=x1,x1=tmp-q*x1;
    tmp=y0,y0=y1,y1=tmp-q*y1;
    tmp=a,a=b,b=tmp%b;
  }
  x=x0,y=y0;
}

matrix mult(matrix A,matrix B)
{
  matrix S;
  memset(S.s,0,sizeof(S.s));
  for(int i=1;i<=m;i++)
    for(int j=1;j<=m;j++)
      for(int k=1;k<=m;k++)
        S.s[i][j]=(S.s[i][j]+A.s[i][k]*B.s[k][j])%9973;
  return S;
}

matrix power(ll x)
{
  int i=0;
  matrix S;
  memset(S.s,0,sizeof(S.s));
  for(int j=1;j<=m;j++) S.s[j][j]=1;
  while(x)
  {
    if (x&1) S=mult(S,M[i]);
    i++;x>>=1;
  }
  return S;
}

ll phi(ll x)
{
  ll s=x;
  for(ll i=2;i*i<=x;i++)
    if (!(x%i))
    {
      s=s/i*(i-1);
      while(!(x%i)) x/=i;
    }
  if (x>1) s=s/x*(x-1);
  return s;
}

void solve(ll x)
{
  matrix S=power(x);
  ll sum=0;
  for(int i=1;i<=m;i++)
    sum=(sum+S.s[i][i])%9973;
  sum=(sum*phi(n/x))%9973;
  ans=(ans+sum)%9973;
}

int main()
{
  scanf("%d",&T);
  while(T--)
  {
    scanf("%lld%d%d",&n,&m,&k);
    ans=0;
    memset(M[0].s,0,sizeof(M[0].s));
    for(int i=1;i<=m;i++)
      for(int j=1;j<=m;j++)
        M[0].s[i][j]=1;
    for(int i=1,a,b;i<=k;i++)
    {
      scanf("%d%d",&a,&b);
      M[0].s[a][b]=M[0].s[b][a]=0;
    }

    for(int i=1;i<=35;i++) M[i]=mult(M[i-1],M[i-1]);
    for(ll i=1;i*i<=n;i++)
      if (n%i==0)
      {
        solve(i);
        if (i!=n/i) solve(n/i);
      }

    ll x0,y0;
    exgcd(n,9973,x0,y0);
    x0=(x0%9973+9973)%9973;
    printf("%lld\n",(x0*ans)%9973);
  }

  return 0;
}