[POJ 2888]Magic Bracelet Burnside与矩阵加速

传送门·Teleport!

题目大意

给你长度为N的项链(可以旋转),问用M种颜色染色有多少种本质不同的方案。
其中还有k对限制关系,表示颜色i和颜色j不能放在一起,答案要取模。
N<=1e9, M<=10

解题报告

先不考虑限制,就变成了 这道题目
考虑限制,使用burnside定理,那么现在难点在于求D(Gi)的值。
由于置换是旋转,循环(转i次的置换下)共有gcd(N, i),那么每个循环选一个元素出来其实就可以代表整个环的状态。
脑补一下(其实是因为懒得去严谨证明。。),发现只要从项链的1号位置开始连续选取gcd(N, i)个元素就刚好能表示整个项链啦!
把每个颜色看成节点,结合限制做成一个状态机A,那么A^gcd(N, i)就是D(Gi),D(Gi)可以通过矩阵加速求得了。
也就是:[latex] D(G_i) = A^{gcd(N,i)} [/latex]
但是N太大了我们不能直接暴力累和,考虑借鉴没有限制的做法。
如果设Gi为转的长度为i的置换,那么有答案等于:
$$ \frac{1}{|G|} \sum_{i=0}^{N-1}D(G_i) = \frac{1}{N} \sum_{d|N} D(G_d) \times \sum_{i=0}^{N-1} (gcd(i, N) == d) $$
继续变换,把i乘上d, 把状态机A代入:
$$ = \frac{1}{N} \sum_{d|N} A^{gcd(N, d)} \times \sum_{i=1}^{N/d-1} (gcd(i, N) == 1) = \frac{1}{N} \sum_{d|N} A^d \times \varphi(\frac{N}{d})$$
矩阵算出来后,考虑这是一个环,所以首尾也要满足限制关系,再稍微搞一下。
然后答案统计完乘个逆元即可,很有价值的题目。

但 — 只是可惜出题人卡常

用c++交tle,g++就AC了什么鬼。。

代码

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <cassert>
const int MAXM = 12;
const int MOD  = 9973;
int N, M, K;
struct Matrix {
    int mat[MAXM][MAXM];
    Matrix () {memset(mat, 0, sizeof mat);}
    Matrix operator * (const Matrix & b) {
        Matrix ret;
        for(int i = 1; i<=M; i++) {
            int *p = ret.mat[i]+1;
            for(int j = 1; j<=M; j++, p++)
                for(int k = 1; k<=M; k++) {
                    *p += (mat[i][k]*b.mat[k][j])%MOD;
                    if(*p >= MOD) *p -= MOD;
                }
        }
        return ret;
    }
} g;
inline int advPow (int a, int b) {
    int ret = 1;
    a %= MOD;
    while(b) {
        if(b & 1) ret = ret*a%MOD;
        a = a*a%MOD;
        b >>= 1;
    }
    return ret;
}
inline Matrix advPow (Matrix a, int b) {
    Matrix ret;
    for(int i = 1; i<=M; i++) ret.mat[i][i] = 1;
    while(b) {
        if(b&1) ret = ret*a;
        a = a*a;
        b >>= 1;
    }
    return ret;
}

inline int varphi (int x) {
    int ret = x;
    for(int i = 2; i*i<=x; i++)
        if(x % i == 0) {
            ret -= ret/i;
            while(x%i == 0) x /= i;
        }
    if(x != 1) ret -= ret/x;
    return ret%MOD;
}

int process (int d) {
    Matrix tmp = advPow(g, d-1);
    int phi = varphi(N/d);
    int ret = 0;
    for(int i = 1; i<=M; i++)
        for(int j = 1; j<=M; j++) {
            if(!g.mat[i][j]) continue ;
            ret += tmp.mat[i][j];
            if(ret >= MOD) ret -= MOD;
        }
    ret = ret*phi%MOD;
    return ret;
}

inline int getInt () {
    char ch; int ret = 0;
    while((ch = getchar()) < '0' || ch > '9') ;
    ret = ch - '0';
    while((ch = getchar()) >= '0' && ch <= '9')
        ret = ret*10 + ch - '0';
    return ret;
}
int main () {
    int cas = getInt();
    for(int tc = 1; tc <= cas; tc++) {
        N = getInt(), M = getInt(), K = getInt();
        for(int i = 1; i<=M; i++)
            for(int j = 1; j<=M; j++)
                g.mat[i][j] = 1;
        for(int i = 1; i<=K; i++) {
            int a, b;
            scanf("%d%d", &a, &b);
            g.mat[a][b] = g.mat[b][a] = 0;
        }
        int sum = 0;
        Matrix tmp;
        for(int d = 1; d*d <= N; d++)
            if(N%d == 0) {
                sum += process(d);
                if(d*d != N) sum += process(N/d);
                sum %= MOD;
            }
        printf("%d\n", sum*advPow(N, MOD-2)%MOD);
    }
}