POJ 2888Magic Bracelet

Posted awhitewall

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了POJ 2888Magic Bracelet相关的知识,希望对你有一定的参考价值。

博主的 bibi 时间

震惊!在加了两种欧拉函数求法优化 + 矩阵惊人优化后的代码竟然比不过朴素欧拉函数 + 矩阵无优化!!!博主吐血调题慌张分析仍然 (TLE)!!!这一切背后,究竟是人性的扭曲还是道德的沦丧?是 * 的爆发还是|饥|^渴的无奈??(我没有打错符号 /xyx

好了我就不 (bibi) 了,只是很好奇 (a+=b) 为什么会比 (a=a+b) 慢这么多,((a=a+b)\\%=mod) 也会比 (a=(a+b)\\%mod) 慢这么多,果然简便的写法会慢些吗 (QwQ)

Description

一个长度为N(N<=1000000000)的环形珍珠项链,有M(M<=10)种颜色,每个珍珠可以任选一种颜色染色,但是规定了K对禁止关系(a,b),即不能存在颜色分别为a,b的珍珠相邻,询问可以有多少种染色方法可以得到本质不同的珍珠项链(通过旋转得到的项链为本质相同)。

Solution

有了 (Polya) 定理后,我们现在需要算的其实就是置换 i(置换 i 表示会将项链旋转 i 下)的上色种类数。

这个玩意可以利用矩阵快速幂来完成。

首先,我们设循环节的个数为 t。我们发现,整条项链可以被分成 (frac{n}{t}) 块,而且每个块内的元素都包含在不同的循环里。显然每个循环内的所有元素颜色必须相同,那么我们只需要找出一个块涂色的所有情况就好了。

下面就是具体的方法:将一个数组 (a[m][m])(m 表示大小而不是下标)都赋为 1,再把不能相邻的颜色赋为 0。

接下来为了方便理解,我画一个示意图:

用这个来理解:

   for(int i = 1; i <= a.n; i ++)
      for(int j = 1; j <= a.n; j ++)
         for(int k = 1; k <= a.n; k ++)
             r.a[i][j] += a.a[i][k] * b.a[k][j];

这是图:(边缘的数字就是行和列)

    1 2 3     1 2 3

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

我们把左边的矩阵看作初始矩阵,右边的矩阵看成转移矩阵。

我们获得 (a[1][1]) 就是初始矩阵的第一行与转移矩阵的第一列相乘。
同理,获得 (a[1][2]) 就是初始矩阵的第一行与转移矩阵的第二列相乘。
再同理,获得 (a[2][1]) 就是就是初始矩阵的第二行与转移矩阵的第一列相乘。

理解了矩阵乘法,之后的定义就好懂了。

我们注意到之前的思路有一个问题:因为本质上是一个环,而这个环是好几个之前分的块组成的。那么两个块之间的元素也是相邻的,也就是第一个块的末元素和第二个块的首元素是相邻的。这点非常重要。

现在来说说 (a[i][j]) 的定义:以 i 为首元素,最后能到达 j 的总方案数。(到达是末元素不一定为 j,但是可以回到 j 的意思,这样就巧妙地处理了环)

这下用手模一下就可以了。(如果还是不懂可以私信或评论哟)

显然我们要求的就是 (sum_{i=0}^{m-1}a[i][i])

最后我们还有一个套路优化,我在这篇博客里证了,不会的可以去康康:【POJ 2154】Color

Code

我觉得真的????。

#include <cstdio>
#include <cstring>

const int mod = 9973, N = 1e5;

int n, m, p[N + 5], phi[N + 5], cnt;
struct Matrix {
    int n, m, a[10][10];
    Matrix() {n = m = 0; memset(a, 0, sizeof a);}
    Matrix operator * (const Matrix t) {
        Matrix r; r.n = n; r.m = t.m;
        for(int i = 0; i < n; ++ i)
            for(int j = 0; j < m; ++ j)
                if(a[i][j])
                    for(int k = 0; k < m; ++ k)
                        r.a[i][k] = (r.a[i][k] + a[i][j] * t.a[j][k]) % mod;
        return r;
    }
}A, B;

Matrix qkpow(Matrix x, int y) {
    Matrix r; r.n = r.m = x.n;
    for(int i = 0; i < r.n; ++ i) r.a[i][i] = 1;
    while(y) {
        if(y & 1) r = r * x;
        x = x * x; y >>= 1;
    }
    return r;
}

int read() {
    int x = 0, f = 1; char s;
    while((s = getchar()) < ‘0‘ || s > ‘9‘) if(s == ‘-‘) f = -1;
    while(s >= ‘0‘ && s <= ‘9‘) {x = (x << 1) + (x << 3) + (s ^ 48); s = getchar();}
    return x * f;
}

int Pow(int x, int y) {
    int ans = 1;
    while(y) {
        if(y & 1) ans = ans * x % mod;
        x = x * x % mod; y >>= 1;
    }
    return ans;
}

int Phi(int n) {
    if(n <= N) return phi[n] % mod;
    int ans = n;
    for(int i = 1; i <= cnt && p[i] * p[i] <= n; ++ i) {
        if(n % p[i] == 0) {
            ans = ans - ans / p[i];
            while(n % p[i] == 0) n /= p[i];
        }
    }
    if(n > 1) ans = ans - ans / n;
    return ans % mod;
}

int get(const int n) {
    int ans = 0;
    B = qkpow(A, n);
    for(int i = 0; i < m; ++ i) ans = ans + B.a[i][i];
    return ans % mod;
}

int Polya() {
    int ans = 0;
    for(int i = 1; i * i <= n; ++ i) {
        if(n % i) continue;
        ans = (ans + Phi(i) * get(n / i)) % mod;
        if(i * i == n) break;
        ans = (ans + Phi(n / i) * get(i)) % mod;
    }
    return ans;
}

void Euler() {
    phi[1] = 1;
    for(int i = 2; i <= N; ++ i) {
        if(! phi[i]) p[++ cnt] = i, phi[i] = i - 1;
        for(int j = 1; j <= cnt && p[j] * i <= N; ++ j) {
            if(i % p[j] == 0) {phi[p[j] * i] = phi[i] * p[j]; break;}
            else phi[p[j] * i] = phi[i] * (p[j] - 1);
        }
    }
}

int main() {
    int u, v;
    Euler();
    for(int T = read(); T; -- T) {
        n = read(), m = read(); A.n = A.m = m;
        for(int i = 0; i < m; ++ i)
            for(int j = 0; j < m; ++ j)
                A.a[i][j] = 1;
        for(int k = read(); k; -- k) {
            u = read() - 1, v = read() - 1;
            A.a[u][v] = A.a[v][u] = 0;
        }
        printf("%d
", Polya() * Pow(n % mod, mod - 2) % mod);
    }
    return 0;
}

以上是关于POJ 2888Magic Bracelet的主要内容,如果未能解决你的问题,请参考以下文章

POJ2888 Magic Bracelet

poj 2888 Magic Bracelet

POJ 2888 Magic Bracelet(burnside引理+矩阵)

POJ2888Magic Bracelet Burnside引理+欧拉函数+矩阵乘法

poj3624 Charm Bracelet(0-1背包 滚动数组)

POJ-3624-Charm Bracelet