【hdu 6960】Necklace of Beads （Polya计数）

题目链接

题目大意

$n(3\le n\le 10^6)$ 个珠子的的环形项链，每个珠子可以涂红色、蓝色或绿色，要求相邻珠子颜色不同，且绿色珠子数目不超过 $k$ 。旋转后相同的方案视作同一种。求方案数。

解题思路

可将不允许旋转时的所有着色方案看作一个集合，我们所要求的则是这个集合在旋转作用群作用下的轨道数目。

根据 Burnside 引理，集合 $X$ 在群 $G$ 作用下的轨道数目为
$$|X/G|=\frac{1}{|G|}\sum _{g\in G}|X^g|$$

记长度为 $n$ 的项链在不允许旋转时，恰好染 $k$ 个绿珠的染色方案数为 $g(n,k)$。（规定$n=1$时为$0$）

则答案为
$$\frac{1}{n}\sum _{i=1}^n\sum _{j=0}^{[k\frac{\gcd (n,i)}{n}]}g(\gcd (n,i),j)=\frac{1}{n}\sum _{d\mid n}\phi (\frac{n}{d})\sum _{j=0}^{[\frac{kd}{n}]}g(d,j)$$

下面的问题就是求 $g(n,k)$。

设在长度为 $n$ 的环上选 $k$ 个不相邻位置的方案数为 $h(n,k)$ ，则 $g(n,k)=2^{k}h(n,k)$

而 $h(n,k)=\{\begin{array}{ll}0& n=1\\ 2& n>1,k=0\\ \left(\genfrac{}{}{0ex}{}{n-k}{k}\right)+\left(\genfrac{}{}{0ex}{}{n-k-1}{k-1}\right)& 其他情况\end{array}$

答案为

$$\frac{1}{n}\sum _{d\mid n}\phi (\frac{n}{d})\sum _{j=0}^{[\frac{kd}{n}]}{2}^{j}h(d,j)$$

代码

#include <bits/stdc++.h>
#define rep(i, l, r) for (int i = l; i <= r; ++i)
#define per(i, r, l) for (int i = r; i >= l; --i)
using namespace std;
const int N = 1000006;
const int mod = 998244353;

int T;
int n, k;
int prime[N], tot = 0;
bool iscmp[N];
int phi[N];
int pw2[N];
int fac[N], invfac[N];
int pw(int x, int y) {
    int ret = 1;
    while (y) {
        if (y & 1) ret = 1ll * ret * x % mod;
        x = 1ll * x * x % mod;
        y >>= 1;
    }
    return ret;
}
int inv(int x) { return pw(x, mod - 2); }

void init() {
    pw2[0] = 1;
    rep(i, 1, 1000000) pw2[i] = 2ll * pw2[i - 1] % mod;

    fac[0] = 1;
    rep(i, 1, 1000000) fac[i] = 1ll * fac[i - 1] * i % mod;

    invfac[1000000] = inv(fac[1000000]);
    per(i, 999999, 0) invfac[i] = 1ll * invfac[i + 1] * (i + 1) % mod;

    phi[1] = 1;
    rep(i, 2, 1000000) {
        if (!iscmp[i]) prime[++tot] = i, phi[i] = i - 1;
        for (int j = 1; j <= tot && i * prime[j] <= 1000000; ++j) {
            iscmp[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            } else {
                phi[i * prime[j]] = phi[i] * (prime[j] - 1);
            }
        }
    }
    // rep(i, 1, 100) printf("%d ", phi[i]);
}
int c(int n, int k) {
    if (k < 0 || k > n) return 0;
    return 1ll * fac[n] * invfac[k] % mod * invfac[n - k] % mod;
}
int h(int n, int k) {
    // printf("  %d %d\n", n, k);
    return (c(n - k, k) + c(n - k - 1, k - 1)) % mod;
}
int solve(int n, int k) {
    if (n == 1) return 1;
    int ret = 0;
    rep(d, 2, n) {  // cyc length
        if (n % d) continue;
        int tmp = 2 * (d % 2 == 0);
        rep(j, 1, 1ll * k * d / n) {
            tmp = (tmp + 1ll * pw2[j] * h(d, j)) % mod;
        }
        // printf("%d %d\n", d, tmp);
        ret = (ret + 1ll * tmp * phi[n / d]) % mod;
    }
    return 1ll * ret * inv(n) % mod;
}
int main() {
    init();
    scanf("%d", &T);
    while (T--) {
        scanf("%d%d", &n, &k);
        printf("%d\n", solve(n, k));
    }
    return 0;
}