LOJ6519. 魔力环（莫比乌斯反演＋生成函数）

题目链接

https://loj.ac/problem/6519

题解

这里给出的解法基于莫比乌斯反演。可以用群论计数的相关方法代替莫比乌斯反演，但两种方法的核心部分是一样的。

环计数的常见套路就是将环视为序列。我们统计所有不同的序列，那么对于最小循环节长度为 (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))，等价于求下面方程的整数解的数量：

[x_0 + x_1 + ... + x_{a - b} = b(0 leq x_i leq k, 0 leq x_0 + x_{a - b} leq k) ]

即被 (a - b) 颗白色珠子划分开的 (a - b + 1) 段黑色珠子的和为 (b)，且满足每连续一段长度不超过 (k) 的限制条件。运用生成函数的知识，求上面方程的解的数量等同于求如下多项式 (h(x)) 中 (x^b) 的系数：

[h(x) = left(sum_{i = 0}^{k} x^i ight) ^ {a - b - 1} left( left(sum_{i = 0}^{k} x^i ight)^2{ m mod} x^{k + 1} ight) ]

进一步地，有：

[h(x) = left(sum_{i = 0}^{k} x^i ight) ^ {a - b - 1} left(sum_{i = 0}^{k} (i +1)x^i ight) ]

由于 (sum_limits{i = 0}^k x^i = frac{1 - x^{k + 1}}{1 - x})，故：

[h(x) = left(frac{1 - x^{k + 1}}{1 - x} ight) ^ {a - b - 1} left(sum_{i = 0}^{k} (i +1)x^i ight) ]

再展开右侧的式子 (sum_limits{i = 0}^k(i + 1)x^i)：

[egin{aligned}sum_{i = 0}^k (i +1)x^i &= x^0 + 2x^1 + 3x^2 + cdots + (k + 1)x^k\ &= (x^0 + x^1 + cdots + x^k) + (x^1 +x^2 + cdots + x^k)+ cdots + x^k \ &= frac{x^0 - x^{k + 1}}{1 - x} + frac{x^1 - x^{k + 1}}{1 - x} + cdots + frac{x^k - x^{k - 1}}{1 - x} \ &= frac{frac{x^0 - x^{k + 1}}{1 - x} - (k + 1)x^{k + 1}}{1 - x} \ &= frac{1 - (k + 2)x^{k + 1} + (k + 1)x^{k + 2}}{(1 - x)^2}end{aligned} ]

因此，我们得到了：

[egin{aligned}h(x) &= left(frac{1 - x^{k + 1}}{1 - x} ight) ^ {a - b - 1} frac{1 - (k + 2)x^{k + 1} + (k + 1)x^{k + 2}}{(1 - x)^2} \ &= frac{(1 - x^{k + 1})^{a - b - 1}}{(1 - x)^{a - b + 1}}(1 - (k + 2)x^{k + 1} + (k + 1)x^{k + 2})end{aligned} ]

其中，((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) = left(sum_{i = 0}^{infty}inom{a - b - 1}{i}(-1)^ix^{(k + 1)i} ight)left(sum_{i = 0}^{infty}inom{a - b + i}{i}x^i ight)(1 - (k + 2)x^{k + 1} + (k + 1)x^{k + 2}) ]

当把 (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;
}