Magic Bracelet（POJ 2888，Burnside 引理 + 矩阵快速幂优化 dp）

最新推荐文章于 2022-07-27 22:29:04 发布

The___Flash

最新推荐文章于 2022-07-27 22:29:04 发布

阅读量340

点赞数

分类专栏： # 数据结构优化 DP

本文链接：https://blog.csdn.net/The___Flash/article/details/112765971

版权

数据结构优化 DP 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

一、题目链接

Magic Bracelet

二、题目大意

给定 $m$ 种不同颜色的魔法珠子，每种颜色的珠子的个数都足够多。

现在要从中挑选 $n$ 个珠子，串成一个环形魔法手链。

魔法珠子之间存在 $k$ 对排斥关系，互相排斥的两种颜色的珠子不能相邻，否则会发生爆炸。（同一种颜色的珠子之间也可能存在排斥）

请问一共可以制作出多少种不同的手链。

注意，如果两个手链经旋转后能够完全重合在一起，对应位置的珠子颜色完全相同，则视为同一种手链。

答案对 $99739973$ 取模。

$\leq T \leq 10 \quad 1 \leq n \leq 10^9 \quad gcd(n,9973)=1 \quad 1 \leq m \leq 10 \quad 0 \leq k \leq \frac{m(m+1)}{2} \quad 1 \leq a,b \leq m$ .

三、分析

由于循环之间存在限制关系，故无法使用 $P\acute olya$ 定理，而只能使用 $B u r n s i d e$ 引理.

旋转操作共有 $n$ 种置换，考虑一个置换 $P_i$ ：顺时针旋转 $i$ 个珠子（ $\leq i \leq n$ ）.

通过二元一次不定方程解得循环长度为 $\frac{n}{gcd(n,i)}$ ，于是循环个数为 $g c d (n, i)$ .

因此置换 $P_i$ 的不动点个数等于长度为 $g c d (n, i)$ 的环的合法方案数，记为 $F (g c d (n, i))$ .

根据 $B u r n s i d e$ 引理可得答案为 $\begin{aligned} ANS = \frac{1}{n}\sum_{i=1}^n{F(gcd(n,i))} \end{aligned}$ .

由于 $n$ 最大 $10^9$ ，直接枚举 $i$ 是无法接受的. 而 $g c d (n, i)$ 最多只有 $1600$ 种，所以我们去枚举 $g c d (n, i)$ .

于是 $\begin{aligned} ANS = \frac{1}{n}\sum_{d|n}\varphi(\frac{n}{d})F(d) \end{aligned}$ ，其中 $\varphi(\frac{n}{d})$ 可以先预处理出 $10^5$ 以内的质数再 $O(\frac{\sqrt{n/d}}{ln(n/d)})$ 求出.

考虑如何计算 $F (d)$ ，记 $c [i] [j]$ 表示颜色 $i$ 和颜色 $j$ 的珠子能否相邻， $f_x[i][j]$ 为第一颗珠子的颜色为 $x$ 且确定了前 $i$ 颗珠子的颜色且第 $i$ 颗珠子颜色为 $j$ 的方案数.

则有递推式 $\begin{aligned} f_x[i][j] = \sum_{c[j][k]} f_x[i - 1][k] \quad (2 \leq i \leq n \quad 1 \leq j,k \leq m) \end{aligned}$ .

于是 $\begin{aligned} F(d) = \sum_{x}\sum_{c[x][j]} f_x[d][j] \end{aligned}$

直接 $d p$ 的时间复杂度是 $O(nm^3)$ 的，不难发现可以使用矩阵快速幂加速 $d p$ 递推.

即

$\left[ \begin{matrix} f_1[i][1] & f_1[i][2] & \cdots & f_1[i][m] \\ f_2[i][1] & f_2[i][2] & \cdots & f_2[i][m] \\ \vdots & \vdots & & \vdots \\ f_m[i][1] & f_m[i][2] & \cdots & f_m[i][m] \end{matrix} \right] = \left[ \begin{matrix} f_1[i-1][1] & f_1[i-1][2] & \cdots & f_1[i-1][m] \\ f_2[i-1][1] & f_2[i-1][2] & \cdots & f_2[i-1][m] \\ \vdots & \vdots & & \vdots \\ f_m[i-1][1] & f_m[i-1][2] & \cdots & f_m[i-1][m] \end{matrix} \right] \left[ \begin{matrix} c[1][1] & c[1][2] & \cdots & c[1][m] \\ c[2][1] & c[2][2] & \cdots & c[2][m] \\ \vdots & \vdots & & \vdots \\ c[m][1] & c[m][2] & \cdots & c[m][m] \end{matrix} \right]$

于是

$\left[ \begin{matrix} f_1[d][1] & f_1[d][2] & \cdots & f_1[d][m] \\ f_2[d][1] & f_2[d][2] & \cdots & f_2[d][m] \\ \vdots & \vdots & & \vdots \\ f_m[d][1] & f_m[d][2] & \cdots & f_m[d][m] \end{matrix} \right] = \left[ \begin{matrix} c[1][1] & c[1][2] & \cdots & c[1][m] \\ c[2][1] & c[2][2] & \cdots & c[2][m] \\ \vdots & \vdots & & \vdots \\ c[m][1] & c[m][2] & \cdots & c[m][m] \end{matrix} \right]^{d-1}$

四、代码实现

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

typedef long long ll;
typedef unsigned long long ull;

const int M = (int)1e5;
const int N = (int)1e5;
const int inf = 0x3f3f3f3f;
const double eps = 1e-8;
const int mod = (int)9973;

int n, m, k;
bool c[10][10];

struct Matrix
{
    int D[10][10];

    Matrix operator*(const Matrix& B)const
    {
        Matrix C;
        for(int i = 0; i < m; ++i)
        {
            for(int j = 0; j < m; ++j)
            {
                C.D[i][j] = 0;
                for(int k = 0; k < m; ++k)
                {
                    C.D[i][j] = (C.D[i][j] + D[i][k] * B.D[k][j]) % mod;
                }
            }
        }
        return C;
    }
};

Matrix quick(Matrix A, int b)
{
    Matrix S;
    if(b == 0)
    {
        for(int i = 0; i < m; ++i)
            for(int j = 0; j < m; ++j) S.D[i][j] = (i == j);
    }
    else
    {
        S = A; --b;
        while(b)
        {
            if(b & 1) S = S * A;
            A = A * A;
            b >>= 1;
        }
    }
    return S;
}

bool is_prime[M + 5];
int prime[M + 5], cnt;

void init()
{
    memset(is_prime, 1, sizeof(is_prime));
    is_prime[0] = is_prime[1] = 0;
    for(int i = 2; i <= M; ++i)
    {
        if(is_prime[i]) prime[++cnt] = i;
        for(int j = 1; j <= cnt && i * prime[j] <= M; ++j)
        {
            is_prime[i * prime[j]] = 0;
            if(i % prime[j] == 0) break;
        }
    }
}

int F(int d)
{
    int ans = 0;
    Matrix A;

    for(int i = 0; i < m; ++i)
    {
        for(int j = 0; j < m; ++j)
        {
            A.D[i][j] = c[i][j];
        }
    }
    A = quick(A, d - 1);
    for(int i = 0; i < m; ++i)
    {
        for(int j = 0; j < m; ++j)
        {
            if(c[i][j]) ans = (ans + A.D[i][j]) % mod;
        }
    }
    return ans;
}

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

int quick(int a, int b, int p = mod)
{
    a %= p;
    int s = 1;
    while(b)
    {
        if(b & 1) s = s * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return s % p;
}

int inv(int n, int p = mod)
{
    return quick(n, p - 2, p);
}

void work()
{
    scanf("%d %d %d", &n, &m, &k);
    for(int i = 0; i < m; ++i) for(int j = 0; j < m; ++j) c[i][j] = 1;
    for(int i = 1, a, b; i <= k; ++i)
    {
        scanf("%d %d", &a, &b);
        --a, --b;
        c[a][b] = c[b][a] = 0;
    }

    int ans = 0;
    for(int d = 1; d * d <= n; ++d)
    {
        if(n % d) continue;
        ans = (ans + phi(n / d) % mod * F(d)) % mod;
        if(d != n / d) ans = (ans + phi(d) % mod * F(n / d)) % mod;
    }
    ans = ans * inv(n) % mod;
    printf("%d\n", ans);
}

int main()
{
//    freopen("input.txt", "r", stdin);
//    freopen("output.txt", "w", stdout);
    init();
    int T; scanf("%d", &T);
    while(T--) work();
//    work();
    return 0;
}