Magic Bracelet POJ - 2888 Burnside引理+数论+DP矩阵优化

题目链接

---

Ginny’s birthday is coming soon. Harry Potter is preparing a birthday present for his new girlfriend. The present is a magic bracelet which consists of n magic beads. The are m kinds of different magic beads. Each kind of beads has its unique characteristic. Stringing many beads together a beautiful circular magic bracelet will be made. As Harry Potter’s friend Hermione has pointed out, beads of certain pairs of kinds will interact with each other and explode, Harry Potter must be very careful to make sure that beads of these pairs are not stringed next to each other.

There infinite beads of each kind. How many different bracelets can Harry make if repetitions produced by rotation around the center of the bracelet are neglected? Find the answer taken modulo 9973.

Input

The first line of the input contains the number of test cases.

Each test cases starts with a line containing three integers n (1 ≤ n ≤ 109, gcd(n, 9973) = 1), m (1 ≤ m ≤ 10), k (1 ≤ k ≤ m(m − 1) ⁄ 2). The next k lines each contain two integers a and b (1 ≤ a, b ≤ m), indicating beads of kind a cannot be stringed to beads of kind b.

Output

Output the answer of each test case on a separate line.

Sample Input

4
3 2 0
3 2 1
1 2
3 2 2
1 1
1 2
3 2 3
1 1
1 2
2 2

Sample Output

4
2
1
0

---

题意:

给出n个珠子的项链和m种珠子；

珠子之间有k对关系，这些珠子不能相邻；

无法通过旋转变成相同的项链视为本质不同；

求本质不同的项链个数，答案对9973取模；

n<=10^9，gcd(n,9973)=1，m<=10；

思路:

Burnside引理

由于只有旋转对称，如果是旋转k个珠子，那么循环节也就是gcd(n,k)=r，枚举k的话是不现实的，那么只有枚举r，即n的所有约数。gcd(n,k)=r，即gcd(n/r,k/r)=1，也就是与n/r互质的数的个数（欧拉函数）就是循环节为r的置换个数。

对循环节为r的情况，需要考虑循环节内部珠子的排列，使它们满足题目要求，还要考虑第一个珠子与最后一个珠子是否满足要求（最后一个珠子的下一个珠子也是下一轮的第一个珠子），由于珠子种类只有10，可以用邻接矩阵map[i][j]表示i,j两种珠子是否能相邻，如果能，map[i][j]=1,反之，map[i][j]=0，这样的话，离散数学老师应该说过，如果用0,1矩阵A来表示无向图的连通情况的话，A^k代表的就是一个点经过k条路后能到达的地方的方法数。

因此，对于循环节为r的情况，A^r就是任意点经过r条路能到达的地方，与之对应的map[i][i]就是一个珠子经过可行路径转了r条路径又回到自己的种数，其实，就是前面说的满足题意的排列数，那么矩阵快速幂加速一下就可以了..

代码:

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
using namespace std;
#define maxm 55555
#define maxn 11
#define mod 9973
int prime[maxm],is_prime[maxm],res;
void get_prime(int n)
{
    res=0;
    memset(is_prime,0,sizeof(is_prime));
    for(int i=2;i<n;i++)
        if(!is_prime[i])
        {
            prime[res++]=i;
            for(int j=i;j<n;j+=i)
                is_prime[j]=1;
        }
}
int mod_pow(int a,int b,int p)
{
    int ans=1;
    a%=p;
    while(b)
    {
        if(b&1) ans=(ans*a)%p;
        a=(a*a)%p;
        b>>=1;
    } 
    return ans;
}
int get_euler(int n)
{
    int ans=n;
    for(int i=0;i<res&&prime[i]*prime[i]<=n;i++)
        if(n%prime[i]==0)
        {
            ans=ans/prime[i]*(prime[i]-1);
            while(n%prime[i]==0)
                n/=prime[i];
        }
    if(n>1)ans=ans/n*(n-1);
    return ans;
} 
struct Mat
{
    int mat[maxn][maxn];
    int row,col;
};
Mat mod_mul(Mat a,Mat b,int p)
{
    Mat ans;
    ans.row=a.row;
    ans.col=b.col;
    memset(ans.mat,0,sizeof(ans.mat));
    for(int i=0;i<ans.row;i++)      
        for(int k=0;k<a.col;k++)
            if(a.mat[i][k])
                for(int j=0;j<ans.col;j++)
                {
                    ans.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                    ans.mat[i][j]%=p;
                }
    return ans;
}
Mat mod_pow(Mat a,int k,int p) 
{
    Mat ans;
    ans.row=a.row;
    ans.col=a.col;
    for(int i=0;i<a.row;i++)
        for(int j=0;j<a.col;j++)
            ans.mat[i][j]=(i==j);
    while(k)
    {
        if(k&1)ans=mod_mul(ans,a,p);
        a=mod_mul(a,a,p);
        k>>=1;
    }
    return ans;
}
int get_track(Mat M,int n,int p)
{
    M=mod_pow(M,n,p);
    int ans=0;
    for(int i=0;i<M.row;i++)
        ans=(ans+M.mat[i][i])%p;
    return ans;
}
int main()  
{  
    get_prime(maxm);
    int T,n,m,k;   
    scanf("%d",&T);  
    while(T--)  
    {  
        scanf("%d%d%d",&n,&m,&k);
        Mat M;
        M.row=M.col=m;
        for(int i=0;i<m;i++)
            for(int j=0;j<m;j++)
                M.mat[i][j]=1;
        while(k--)
        {
            int a,b;
            scanf("%d%d",&a,&b);
            M.mat[a-1][b-1]=M.mat[b-1][a-1]=0;
        }
        int ans=0;
        for(int i=1;i*i<=n;i++)
            if(n%i==0)
            {
                if(i*i!=n)
                    ans=(ans+get_euler(i)%mod*get_track(M,n/i,mod)%mod+get_euler(n/i)%mod*get_track(M,i,mod)%mod)%mod;
                else 
                   ans=(ans+get_euler(i)%mod*get_track(M,n/i,mod)%mod)%mod;
            }
        ans=ans*mod_pow(n,mod-2,mod)%mod;
        printf("%d\n",ans);
    }  
    return 0;  
}