题目链接
题解
这里给出的解法基于莫比乌斯反演。可以用群论计数的相关方法代替莫比乌斯反演，但两种方法的核心部分是一样的。
环计数的常见套路就是将环视为序列。我们统计所有不同的序列，那么对于最小循环节长度为 (d) 的序列对应的环就会被统计 (d) 次。因此假设最小循环节长度为 (x) 的合法序列数为 (f(x))，那么答案即为 (sum_limits{d | { m gcd}(n, m)} frac{1}{d}f(d))。
直接求 (f(x)) 不好求，使用莫比乌斯反演转化为求 (g(x))， (g(x)) 表示最小循环节长度为 (x) 的因子的合法序列个数。那么有 (g(x) = sum_limits{d | x} f(d))，即 (f(x) = sum_limits{d | x} mu(d) g(frac{x}{d}))。
我们假设整个长度为 (n) 的序列被分成了 (l) 段，每一段都是原序列的一个循环节。令 (a = frac{n}{l}, b = frac{m}{l})（(a) 即为每一段的珠子数，(b) 即为每一段的黑色珠子数量），那么我们求 (g(a))，等价于求下面方程的整数解的数量：
即被 (a - b) 颗白色珠子划分开的 (a - b + 1) 段黑色珠子的和为 (b)，且满足每连续一段长度不超过 (k) 的限制条件。运用生成函数的知识，求上面方程的解的数量等同于求如下多项式 (h(x)) 中 (x^b) 的系数：
进一步地，有：
由于 (sum_limits{i = 0}^k x^i = frac{1 - x^{k + 1}}{1 - x})，故：
再展开右侧的式子 (sum_limits{i = 0}^k(i + 1)x^i)：
因此，我们得到了：
其中，((1 - x^{k + 1})^{a - b - 1}) 可化为 (sum_limits{i = 0}^{infty}inom{a - b - 1}{i}(-1)^ix^{(k + 1)i})，而 (frac{1}{(1 - x)^{a - b + 1}}) 即 ((1 - x)^{-(a - b + 1)})，可通过负整数次幂的二项式定理化为 (sum_limits{i = 0}^{infty}inom{a - b + i}{i}x^i)，那么可得：
当把 (h(x)) 化成该形式后，要求 (h(x)) 中 (x^b) 的系数就变得非常简单了。记 (s_1 = sum_limits{(k + 1)i + j = b}(-1)^iinom{a - b - 1}{i}inom{a - b+ j}{j})，(s_2 = (k + 2)sum_limits{(k + 1)i + j = b - k - 1}(-1)^iinom{a - b - 1}{i}inom{a - b+ j}{j})，(s_3 = (k + 1)sum_limits{(k + 1)i + j = b - k - 2}(-1)^iinom{a - b - 1}{i}inom{a - b+ j}{j})，(x^b) 的系数为 (w)，那么有 (w = s_1 - s_2 + s_3)。
求 (s_1, s_2, s_3) 只需按照 (s_1, s_2, s_3) 的式子枚举 (i) 即可，因为 (i) 确定 (j) 也就确定了。因此，我们可以在 (frac{b}{k + 1}) 的时间内求出 (h(x)) 中 (x^b) 的系数。
除去反演部分，我们就能够在 (frac{sigma(n)}{k + 1}) 的时间内解决此题，其中，(sigma(n)) 表示 (n) 的约数和。由于 (sigma(n)) 可近似看作 (n { m log} { m log} n)，接近线性，因此时间复杂度是非常优秀的。
代码
#include<bits/stdc++.h>
using namespace std;
#define X first
#define Y second
#define mp make_pair
#define pb push_back
#define debug(...) fprintf(stderr, __VA_ARGS__)
typedef long long ll;
typedef long double ld;
typedef unsigned int uint;
typedef pair<int, int> pii;
typedef unsigned long long ull;
template<typename T> inline void read(T& x) {
char c = getchar();
bool f = false;
for (x = 0; !isdigit(c); c = getchar()) {
if (c == '-') {
f = true;
}
}
for (; isdigit(c); c = getchar()) {
x = x * 10 + c - '0';
}
if (f) {
x = -x;
}
}
template<typename T, typename... U> inline void read(T& x, U&... y) {
read(x), read(y...);
}
template<typename T> inline bool checkMax(T& a, const T& b) {
return a < b ? a = b, true : false;
}
template<typename T> inline bool checkMin(T& a, const T& b) {
return a > b ? a = b, true : false;
}
const int N = 1e5 + 10, mod = 998244353, G = 3;
int n, m, k;
inline void add(int& x, int y) {
if (y < 0) {
y += mod;
}
x = (x + y) % mod;
}
inline void mul(int& x, int y) {
x = 1ll * x * y % mod;
}
inline int qpow(int v, int p) {
int res = 1;
for (; p; p >>= 1, mul(v, v)) {
if (p & 1) {
mul(res, v);
}
}
return res;
}
inline int gcd(int x, int y) {
return !y ? x : gcd(y, x % y);
}
int p[N], pri[N], mu[N], fac[N], invfac[N], c;
inline int binom(int n, int m) {
return n < 0 || m < 0 || n < m ? 0 : 1ll * fac[n] * invfac[m] % mod * invfac[n - m] % mod;
}
void sieve(int n) {
mu[1] = 1;
for (register int i = 2; i <= n; ++i) {
p[i] = 1;
}
for (register int i = 2; i <= n; ++i) {
if (p[i]) {
mu[i] = -1, pri[++c] = i;
}
for (register int j = 1, d; j <= c && (d = pri[j] * i) <= n; ++j) {
p[d] = 0;
if (i % pri[j] == 0) {
mu[d] = 0; break;
} else {
mu[d] = -mu[i];
}
}
}
fac[0] = invfac[0] = 1;
for (register int i = 1; i <= n; ++i) {
mul(fac[i] = fac[i - 1], i);
}
invfac[n] = qpow(fac[n], mod - 2);
for (register int i = n - 1; i; --i) {
mul(invfac[i] = invfac[i + 1], i + 1);
}
}
inline int calc(int n, int m) {
int res = 0;
for (register int i = 0; i * (k + 1) <= m; ++i) {
int j = m - i * (k + 1);
add(res, 1ll * binom(n - m - 1, i) * binom(n - m + j, j) % mod * (i & 1 ? mod - 1 : 1) % mod);
j = m - i * (k + 1) - k - 1;
if (j >= 0) {
add(res, 1ll * (k + 2) * binom(n - m - 1, i) % mod * binom(n - m + j, j) % mod * (i & 1 ? 1 : mod - 1) % mod);
}
j = m - i * (k + 1) - k - 2;
if (j >= 0) {
add(res, 1ll * (k + 1) * binom(n - m - 1, i) % mod * binom(n - m + j, j) % mod * (i & 1 ? mod - 1 : 1) % mod);
}
}
return res;
}
int f[N], g[N];
int main() {
read(n, m, k);
sieve(n);
int d = gcd(n, m);
for (register int i = 1; i <= d; ++i) {
if (d % i) {
continue;
}
g[n / i] = calc(n / i, m / i);
}
for (register int i = 1; i <= n; ++i) {
for (register int j = i; j <= n; j += i) {
add(f[j], mu[i] * g[j / i]);
}
}
int ans = 0;
for (register int i = 1; i <= n; ++i) {
if (n % i == 0) {
add(ans, 1ll * f[i] * qpow(i, mod - 2) % mod);
}
}
printf("%d
", ans);
return 0;
}