POJ 2888 Magic Bracelet 有限制的polya

题意：给定n（n <= 10^9）颗珠子和m种颜色，这n颗珠子可以组成一串项链，每颗珠子可以用m中颜色来填充，如果两种涂色方法通过旋转项链可以得到，则两种涂色方法等价。现在再给定K组限制，每组限制a、b代表颜色a和颜色b不能涂在相邻的珠子上面。问一共有多少种涂色方法。

解题思路：这题和POJ2154 题目大致相同，不同的是这里多了一个限制条件，即可能两种颜色不能涂在相邻的珠子上面。对于这种情况我们可以通过矩阵连乘得到，先初始化矩阵array[i][j]为1.如果颜色a和颜色b不能涂在相邻的珠子，那么array[a][b] = array[b][a] = 0; 对于具有n/L个循环节的置换(L为每个循环节的长度)。先求出array[][]的n/L次幂，然后将这个矩阵的array[i][i] (1<=i<=m)全部加起来即为这种置换下涂色不变的方法数。

以上转自：http://hi.baidu.com/ofeitian/blog/item/694770caf3da47f052664fca.html

补充：

A的k次幂表示图的长度为k的通路，回到自身则表示长度为k的回路。学过离散数学都知道了。

但是为什么可以用这样一个通路来求”置换保持着色不变的个数“？

假如一个10阶置换f可以分解成5个轮换：[1,2] [3,4] [5,6] [7,8] [9,10]。同时我们又求出一个长度为5的颜色回路（a,b,c,d,e)。

那么我们就可以用a给[1,2]上色，b给[3,4]上色，c给[5,6]上色····。

所以每一个回路可以提供一种上色方式，并且这样的上色方式在通过f置换后保持不变。

所以把总的回路个数加起来就得到了，保持不变的总着色数。

另外千万别忘记除以总的置换个数，由于gcd(n,9973) = 1，求一下逆元就OK···

#include<cmath>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cstdio>
using namespace std;

const int MAXN = 1000000;
const int modulo = 9973;
int a[MAXN], p[MAXN], pn;

class CMatrix
{
public:
    int data[11][11], n;
    CMatrix (int, int);
    CMatrix operator* ( const CMatrix& ) const;
    CMatrix power( int ) const;
};

CMatrix::CMatrix ( int f, int size )
{
    n = size;
    memset(data,0,sizeof(data));
    if ( f == 0 ) return;
    for ( int i = 0; i < 11; i++ ) data[i][i] = 1;
}

//注意稀疏矩阵的优化
CMatrix CMatrix::operator* ( const CMatrix ¶m ) const
{
    CMatrix ret(0, n);
    for ( int i = 0; i < n; i++ )
        for ( int j = 0; j < n; j++ )
        {
            for ( int k = 0; k < n; k++ )
                ret.data[i][j] += data[i][k] * param.data[k][j];
            ret.data[i][j] %= modulo;
        }
    return ret;
}

//快速求幂
CMatrix CMatrix::power ( int exp ) const
{
    CMatrix ret(1, n);
    CMatrix tmp = (*this);
    while ( exp > 0 )
    {
        if ( exp & 1 )
            ret = ret * tmp;
        tmp = tmp * tmp;
        exp >>= 1;
    }
    return ret;
}

void Prime()
{
    int i, j; pn = 0;
    memset(a,0,sizeof(a));
    for ( i = 2; i < MAXN; i++ )
    {
        if ( !a[i] ) p[pn++] = i;
        for(j = 0; j < pn && i*p[j] < MAXN && (p[j]<=a[i] || a[i]==0); j++)
            a[i*p[j]] = p[j];
    }
}

int Euler ( int n )
{
    int i, ret = n;
    for ( i = 0; i < pn && p[i] * p[i] <= n; i++ )
    {
        if ( n % p[i] == 0 )
        {
            ret = ret - ret / p[i];
            while ( n % p[i] == 0 ) n /= p[i];
        }
    }
    if ( n > 1 )
        ret = ret - ret / n;
    return ret % modulo;
}

int loop ( const CMatrix &A, int l )
{
    int ret = 0;
    CMatrix tmp = A.power(l);
    for ( int i = 0; i < tmp.n; i++ )
        ret += tmp.data[i][i];
    return ret % modulo;
}

int ext_gcd ( int a, int b, int& x, int& y )
{
    int ret, tmp;
    if ( b == 0 )
    {
        x = 1, y = 0;
        return a;
    }
    ret = ext_gcd(b, a % b, x, y);
    tmp = x, x = y, y = tmp - a / b * y;
    return ret;
}

int inverse ( int n )
{
    int x, y;
    ext_gcd (n,modulo,x,y);
    x = x % modulo;
    return x >= 0 ? x : x + modulo;
}

int polya ( const CMatrix &A, int n, int m )
{
    int ret = 0;
    for( int l = 1; l * l <= n; l++ )
    {
       if ( n % l ) continue;
       ret = (ret + Euler(l) * loop(A,n/l)) % modulo;
       if ( l * l == n ) break;
       ret = (ret + Euler(n/l) * loop(A,l)) % modulo;
    }
    return ret * inverse(n) % modulo;
}

int main()
{
    Prime();
    int n, m, k, cs, x, y;
    scanf("%d",&cs);
    while ( cs-- )
    {
        scanf("%d%d%d",&n,&m,&k);
        CMatrix A(0,m);
        for ( int i = 0; i < m; i++ )
            for ( int j = 0; j < m; j++ )
                A.data[i][j] = 1;
        while ( k-- )
        {
            scanf("%d%d",&x,&y);
            A.data[x-1][y-1] = A.data[y-1][x-1] = 0;
        }
        printf("%d\n",polya(A,n,m));
    }
    return 0;
}