POJ2888 Magic Bracelet（Burnside引理+矩阵快速幂优化DP）

POJ2888 Magic Bracelet

原题地址：http://poj.org/problem?id=2888

题意：
T组数据。
有n个珠子串成一个环，珠子共有m种，有k个要求说第x种和第y种珠子不能相邻。
我们认为两个可以通过旋转得到的环是同一种，问一共能串出多少本质不同的环。对9973取模。

数据范围
1 ≤ n ≤ 1e9， gcd(n, 9973) = 1，1 ≤ m ≤ 10， 1 ≤ k ≤ m(m − 1) ⁄ 2

题解：
这篇题解讲的很清楚。

重新学习burnside引理：
$L=\frac{1}{|G|}\sum_{k=1}^{n}|Z_k|=\frac{1}{|G|}\sum_{i=1}^{s}D(a_i)$
倘若没有限制的话，由此可以得到Polya定理：
$L = \frac{1}{|G|}(m^{c(g_1)}+m^{c(g_2)}+\dots+m^{c(g_s)})$

对于此题，对于n个珠子的一个环，他的置换有 0,1,2…n-1共n个，
要求每个置换的$D(a_i)$，置换$a_i$下不变元素的个数。

观察发现，一个长为n的环的置换i，其循环节个数为$gcd(n,i)$
例如：
对于一个长为6的环，其置换$\{ ^{1,2,3,4,5,6}_{5,6,1,2,3,4}\}$

因为1可以转到3，3可以转到5，2可以转到4，4可以转到6，
有两个循环{1,3,5}{2,4,6}，每个循环的长度为3。$gcd(6,2)=2$。每两个分成一组每组就是每个循环中各取了一个元素。
如果要是在这个置换下不变，那么一个循环必须染同一种颜色。

于是求在置换i下不变元素的个数，就是求这个小长为$gcd(n,i)$的小块的染色方案，记为$f(gcd(n,i))$，这个可以用矩阵快速幂优化DP求。
答案$=\frac{1}{n}\sum_{i=1}^{n}f( gcd(n,i))$
$=\frac{1}{n}\sum_{d|n}f(d)\phi(n/d)$

我觉着这个复杂度有点虚，但是因为d的个数最多不到2000，所以没没毛病。

代码：

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=11;
const int mod=9973;
struct Matrix
{
    int n,a[N][N];
    void init() {for(int i=1;i<=10;i++) a[i][i]=1; n=10;}
}e,a;
Matrix operator*(const Matrix &A,const Matrix &B)
{
    Matrix C; C.n=A.n;
    memset(C.a,0,sizeof(C.a));
    for(int i=1;i<=A.n;i++) 
    for(int j=1;j<=A.n;j++)
    for(int k=1;k<=A.n;k++)
    C.a[i][j]+=A.a[i][k]*B.a[k][j];                         
    for(int i=1;i<=A.n;i++) 
    for(int j=1;j<=A.n;j++)
    C.a[i][j]%=mod;
    return C;
}
int modpow(int A,int B)
{
    int ret=1; int base=A;
    for(;B;B>>=1)
    {
        if(B&1) ret=(1LL*ret*base)%mod;
        base=(1LL*base*base)%mod;
    }
    return ret;
}
int T,n,m,k;
int phi(int x)
{
    int ret=x; 
    for(int i=2;i*i<=x;i++)
    {
        if(x%i==0)
        {
            ret=ret/i; ret=ret*(i-1);
            while(x%i==0) x=x/i;
        }
    }
    if(x>1) {ret=ret/x; ret=ret*(x-1);}
    return ret;
}
int cal(Matrix base,int x)
{
    Matrix ret=e; ret.n=m;  int ans=0;
    for(;x;x>>=1)
    {
        if(x&1) ret=ret*base;
        base=base*base;
    }
    for(int i=1;i<=m;i++) {ans=(ans+ret.a[i][i]); if(ans>=mod) ans-=mod;}
    return ans; 
}
int getans(int x)
{
    int ans=cal(a,x);
    ans=(1LL*ans*phi(n/x))%mod;
    return ans;
}
int main()
{   
    scanf("%d",&T);
    e.init(); int ret;
    while(T--)
    {
        scanf("%d%d%d",&n,&m,&k);
        a.n=m; ret=0;
        for(int i=1;i<=m;i++) for(int j=1;j<=m;j++) a.a[i][j]=1;
        for(int i=1;i<=k;i++)
        {
            int x,y; scanf("%d%d",&x,&y);
            a.a[x][y]=a.a[y][x]=0;
        }
        for(int i=1;i*i<=n;i++)
        if(n%i==0) {ret+=getans(i); if(n/i!=i) ret+=getans(n/i); if(ret>=mod) ret%=mod; }           
        ret=(1LL*ret*modpow(n,mod-2))%mod;
        printf("%d\n",ret);
    }
    return 0;
}