POJ 2888Magic Bracelet

Posted 2021-02-28 awhitewall

博主的 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;
}