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

测试地址：Magic Bracelet
题目大意：要用$N(\le 10^9)$颗珠子连成环形的手镯，共有$M(\le 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)=\sum_{k=1}^{M} m(j,k)\times 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(N^2M)$，显然不能通过此题，需要考虑优化。此题需要从两个方面入手，一是计算Burnside公式的复杂度，而是动态规划的复杂度。
首先优化计算Burnside公式的时间复杂度。我们知道$\gcd(N,i)|N$，而且当$\gcd(N,i)$相同时，$|C(f)|$相同，那么我们完全可以枚举$N$的因子$d$，计算$\gcd(N,i)=d$时的$|C(f)|$，再乘上使得$\gcd(N,i)=d$的$i$的个数，累加起来得到答案。后面那个个数显然就是$\varphi(N/d)$了（不懂的可以去看看我写的POJ2154的题解），那么枚举的时间复杂度就大大降低了。
接下来优化DP的时间复杂度，注意到$M$很小，那么二元函数$m$显然可以构造成一个矩阵$M'$，如果我们把$f(i,1),f(i,2),...,f(i,M)$从上到下排成一个列向量$F_i$，那么可以看出：$F_i=M'F_{i-1}$，所以$F_i=M'^{i-1}F_1$，于是我们可以先用矩阵快速幂算出$M'^d$，然后就可以算$F_{d+1}$了，但其实我们没有必要算出全部的$F_{d+1}$，因为我们只要求$f(d+1,k)$，这是第$k$行第$1$列的元素，那么就用$M'^{d}$的第$k$行与$F_1$的第$1$列做运算，根据定义，$F_1$中只有第$k$行为$1$，其余都为$0$，所以$f(d+1,k)=M'^{d}_{kk}$，那么$|C(f)|$就应该等于$M'^d$对角线上的元素和。
经过了以上两个优化，外层枚举的复杂度降到差不多$O(\sqrt N)$，DP的时间复杂度从$O(NM)$降到$O(M^3\log N)$，因此总的时间复杂度为$O(M^3\sqrt N\log N)$，可以通过此题。
以下是本人代码：

#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;
}