【HDU2888】BurningSide定理+因子分解+求逆元+矩阵快速幂

Magic Bracelet

Time Limit: 2000MS		Memory Limit: 131072K
Total Submissions: 4466		Accepted: 1458

Description

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 ≤ 10⁹, 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

Source

POJ Monthly--2006.07.30, cuiaoxiang

题意：用m种不同颜色的珠子连成一条长为n的项链，其中，有k对珠子不能相邻，问总共有多少种（mod 9973）n<10^9,m<=10

题解：组合计数也就burning和polya了，这题用的是Burning Side。

考虑在一种置换f下的稳定核方法，由于只有旋转对称，如果是旋转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条路径又回到自己的种数，其实，就是前面说的满足题意的排列数，矩阵乘法可以分治加个速，对于n的约数随便求一下，这道题就出来了。

题意：m种珠子串成由n个珠子构成的环，并且有一些限制，比如第i种与第j种不能相邻，旋转后相同算是同一种方案，求方案数。(m<=10,n<=10^9)

如果没有限制条件，可以很容易用Burnside定理+容斥得到答案。

如果只考虑限制，观察发现，m很小，n很大。那么可以通过快速幂得到从第i种到第j种，共有k个珠子的方案数。

再回到Burnside定理，显然有n种置换，现在问题是如何求每种置换能保持不变的着色方案数：

在【POJ】2154 Color已经推出，第i种置换（即旋转360/n*i度）会有GCD(n,i)种不同的颜色，循环节长度为GCD(n,i)。

就可以用快速幂得到从第i种回到第i种，共有GCD(n,i)个珠子的方案数。带回Burnside定理+容斥，问题就解决了。

【POJ】2409 Let it Bead->【POJ】2154 Color->【POJ】2888 Magic Bracelet。

#define DeBUG
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <string>
#include <set>
#include <sstream>
#include <map>
#include <list>
#include <bitset>
using namespace std ;
#define zero {0}
#define INF 0x3f3f3f3f
#define EPS 1e-6
#define TRUE true
#define FALSE false
typedef long long LL;
const double PI = acos(-1.0);
//#pragma comment(linker, "/STACK:102400000,102400000")
inline int sgn(double x)
{
    return fabs(x) < EPS ? 0 : (x < 0 ? -1 : 1);
}
#define N 32000
#define mod 9973
int n, m;
bool p[N];
std::vector<int> prime;
std::vector<int> factor;
std::vector<int> primefactor;
const int MAXN = 12;
struct Matrix
{
    int mat[MAXN][MAXN];
    void Zero()
    {
        memset(mat, 0, sizeof(mat));
    }
    void Unit()
    {
        memset(mat, 0, sizeof(mat));
        for (int i = 0; i < MAXN; i++)
            mat[i][i] = 1;
    }
} g;
Matrix operator*(Matrix &a, Matrix &b)
{
    Matrix tmp;
    tmp.Zero();
    for (int k = 0; k < MAXN; k++)
    {
        for (int i = 0; i < MAXN; i++)
        {
            if (!a.mat[i][k])
                continue;
            for (int j = 0; j < MAXN; j++)
            {
                tmp.mat[i][j] += a.mat[i][k] * b.mat[k][j] % mod;
            }
        }
    }
    return tmp;
}
Matrix operator ^(Matrix a, int k)
{
    Matrix tmp;
    tmp.Unit();
    for (; k; k >>= 1)
    {
        if (k & 1)
            tmp = tmp * a;
        a = a * a;
    }
    return tmp;
}
void Init()
{
    prime.clear();
    memset(p, 1, sizeof(p));
    for (int i = 2; i < 180; i++)
    {
        if (p[i])
            for (int j = i * i; j < N; j += i)
                p[j] = false;
    }
    for (int i = 2; i < N; i++)
    {
        if (p[i])
            prime.push_back(i);
    }
}
LL Ext_gcd(LL a, LL b, LL &x, LL &y)
{
    if (b == 0)
    {
        x = 1, y = 0;
        return a;
    }
    LL ret = Ext_gcd(b, a % b, y, x);
    y -= a / b * x;
    return ret;
}
LL Inv(LL a, int m)   ///求逆元a相对于m
{
    LL d, x, y, t = (LL)m;
    d = Ext_gcd(a, t, x, y);
    if (d == 1) return (x % t + t) % t;
    return -1;
}
void Factor(int x)
{
    int i, tmp;
    factor.clear();
    tmp = (int)(sqrt((double)x) + EPS);
    for (i = 1; i <= tmp; i++)
    {
        if (x % i == 0)
        {
            factor.push_back(i);
            if (i == tmp && i * i == x)
                continue;
            factor.push_back(n / i);
        }
    }
}
int NoChange(int x)
{
    int i, ans;
    Matrix tmp;
    tmp = g ^ x;
    for (i = ans = 0; i < m; i++)
    {
        ans += tmp.mat[i][i] % mod;
    }
    return ans;
}
void Prime(int x)
{
    int i, tmp;
    primefactor.clear();
    tmp = (int)(sqrt((double)x) + EPS);
    for (i = 0; prime[i] <= tmp&&i<prime.size(); i++)
    {
        if (x % prime[i] == 0)
        {
            primefactor.push_back(prime[i]);
            while (x % prime[i] == 0)
            {
                x /= prime[i];
            }
        }
    }
    if (x > 1)
        primefactor.push_back(x);
}
int Mul(int x, int &k)
{
    int i, ans;
    ans = 1;
    for (i = k = 0; x; x >>= 1, i++)
    {
        if (x & 1)
        {
            k++;
            ans *= primefactor[i];
        }
    }
    return ans;
}
int Count(int x)
{
    int i, j, t, ans, tmp;
    Prime(x);
    ans = 0;
    t = (int)primefactor.size();
    for (i = 1; i < (1 << t); i++)
    {
        tmp = Mul(i, j);
        if (j & 1)
            ans += x / tmp;
        else
            ans -= x / tmp;
    }
    return (x - ans) % mod;
}
int Burnside()
{
    int i, ans;
    Factor(n);
    for (i = ans = 0; i < (int)factor.size(); i++)
    {
        ans += Count(n / factor[i]) * NoChange(factor[i]) % mod;
    }
    return ans * Inv(n, mod) % mod;
}
int main()
{
#ifdef DeBUGs
    freopen("C:\\Users\\Sky\\Desktop\\1.in", "r", stdin);
#endif
    int T;
    int x, y, k;
    Init();
    scanf("%d", &T);
    while (T--)
    {
        scanf("%d%d%d", &n, &m, &k);
        g.Zero();
        for (int i = 0; i < m; i++)
        {
            for (int j = 0; j < m; j++)
                g.mat[i][j] = 1;
        }
        while (k--)
        {
            scanf("%d%d", &x, &y);
            x--;
            y--;
            g.mat[x][y] = g.mat[y][x] = 0;
        }
        printf("%d\n", Burnside());
    }

    return 0;
}