题解 OpenJ_Bailian-4052【Necklace】

好题分享

GDSYZX cjh

题目描述

https://vjudge.csgrandeur.cn/problem/OpenJ_Bailian-4052

你将制作一条项链。项链由 \(m\) 颗宝石组成，有 \(n\) 种宝石可供选用。对于第 \(i\) 种宝石，它在项链上的出现次数是如下四种限制中的一种：

• A 要求第 \(i\) 种宝石的出现次数为奇数。
• B 要求第 \(i\) 种宝石的出现次数为偶数。
• C 要求第 \(i\) 种宝石至少出现一次。
• D 不对出现次数作要求。

求可以制作多少种不同的项链。两条项链不同，当且仅当组成项链的宝石排成一排后不能旋转同构。

模 \(10007\)，多测 \(T\leq 50, n\leq 25, m\leq 10^5\)。

Polya 原理

置换群

考虑所有宝石都是 D 限制的情况。用置换群的概念刻画旋转同构，我们记一个群 \(G\)，\(G\) 里面是所有的旋转置换，例如当 \(m=5\) 时就有

\[G=\left\{\binom{1, 2, 3, 4, 5}{1, 2, 3, 4, 5}, \binom{1, 2, 3, 4, 5}{2, 3, 4, 5, 1}, \binom{1, 2, 3, 4, 5}{3, 4, 5, 1, 2}, \binom{1, 2, 3, 4, 5}{4, 5, 1, 2, 3}, \binom{1, 2, 3, 4, 5}{5, 1, 2, 3, 4}\right\} \]

这五个置换。同时要注意群的其他性质。将运算定义为置换的复合，则这个群有单位元 \(\binom{1, 2, 3, 4, 5}{1, 2, 3, 4, 5}\)，运算满足结合律，\(G\) 中元素互相复合后还在 \(G\) 中，\(G\) 确实是一个群。

Polya 原理

Polya 原理告诉我们，若想知道项链在 \(G\) 中置换的作用下有多少个本质不同的项链，则它是：

\[ans = \frac{1}{|G|}\sum_{g\in G}cyc(g) \]

其中 \(cyc(g)\) 表示有多少条项链使得被 \(g\) 作用之后仍不变。对于旋转置换 \(i\to (i+d)\bmod n\)，\(cyc(g)\) 就是将项链（这时候看作链）拆成 \(n/\gcd(n, d)\) 段，要求每一段都相等，这样旋转之后不变。

指数生成函数

不考虑旋转同构的限制，如何求出方案数？使用指数生成函数 EGF。

引入

数列 \(a = \{a_0, a_1, a_2, \cdots\}\) 的 EGF 为 \(A(x) = \sum_{i=0}^\infty a_i\dfrac{x^i}{i!}\)。两个 EGF 的二项卷积 \(A*B\) 定义为：

\[[x^n/n!](A*B)(x)=\sum_{i=0}^n\binom n i[x^i/i!]A(x)[x^{n-i}/(n-i)!]B(x). \]

也就是 \(c_n=\sum_{i}\binom n ia_ib_{n-i}\)。这个卷积的意义在于，如果 \(a_i, b_i\) 表示某种大小为 \(i\) 的组合对象的方案数，那么 \(c_i\) 就表示大小为 \(i\)，由 \(a, b\) “归并”而成的组合对象的方案数。例如 AA 和 CD 归并有 \(\binom{2+2}{2}=6\) 种结果：AACD, ACAD, ACDA, CADA, CAAD, CDAA。就是原来的每一种方案，保持原来顺序的前提下插入到另一种方案，而 EGF 刻画的就是这个东西。

另外有 \(e^x=\sum_{i=0}^\infty x_i/i!\)。

应用

考虑题目中的宝石怎么刻画限制。

• 出现次数为奇数，可以认为是所有奇数个数都有一个方案，偶数个数没有方案，所以尝试找到 \(\{0, 1, 0, 1, 0, 1, \cdots\}\) 的 EGF。它为 \(\dfrac{e^x-e^{-x}}{2}\)，可以感受一下正确性。
• 出现次数为偶数：\(\dfrac{e^x+e^{-x}}{2}\)。
• 至少出现一次：\(e^x-1\)。
• 没有限制：\(e^x\)。

记限制 ABCD 的出现次数为 \(a, b, c, d\)，则答案为

\[[x^n/n!]\left(\frac{e^x-e^{-x}}{2}\right)^a\left(\frac{e^x+e^{-x}}{2}\right)^b\left(e^x-1\right)^c e^{dx} \]

使用二项式定理展开这些式子并暴力相乘，最终得到类似于 \(\sum_kc_ke^{d_ix}\) 的东西。提取 \([x^n/n!]\) 系数，实际上就是分开每一项，因为 \([x^n/n!]ce^{dx}=cd^n\)，所以答案就是 \(\sum_k c_kd_k^n\)。

细节

• 最终统计答案的时候大概是这样：\(ans = \sum_{i=1}^mcalc(\gcd(i, m))/m\)。其中 \(calc(n)\) 计算项链长度为 \(n\) 时的答案，不考虑旋转同构。
• 如果 \(10007|m\)，\(m\) 没有逆元，考虑先计算答案模 \(10007^2\) 的结果 \(ans'\)，然后计算 \(ans=\frac{ans'/10007}{m/10007}\)，这时候 \(m/10007\) 在 \(\pmod {10007^2}\) 意义下有逆元，可以计算。
• 时间复杂度视具体实现大约为 \(O(Tn^2\log n\sqrt m)\)。

code

#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define debug(...) void(0)
#endif
typedef long long LL;
template <unsigned umod>
struct modint {
  static constexpr int mod = umod;
  unsigned v;
  modint() : v(0) {}
  template <class T>
  modint(T x) {
    x %= mod;
    if (x < 0) x += mod;
    v = x;
  }
  modint operator+() const { return *this; }
  modint operator-() const { return modint() - *this; }
  friend int raw(const modint &self) { return self.v; }
  modint &operator+=(const modint &rhs) {
    v += rhs.v;
    if (v >= umod) v -= umod;
    return *this;
  }
  modint &operator-=(const modint &rhs) {
    v -= rhs.v;
    if (v >= umod) v += umod;
    return *this;
  }
  modint &operator*=(const modint &rhs) {
    v = static_cast<unsigned>(1ull * v * rhs.v % umod);
    return *this;
  }
  modint &operator/=(const modint &rhs) {
    assert(rhs.v);
    return *this *= qpow(rhs, mod - 2);
  }
  template <class T>
  friend modint qpow(modint a, T b) {
    modint r = 1;
    for (; b; b >>= 1, a *= a)
      if (b & 1) r *= a;
    return r;
  }
  friend modint operator+(modint lhs, const modint &rhs) { return lhs += rhs; }
  friend modint operator-(modint lhs, const modint &rhs) { return lhs -= rhs; }
  friend modint operator*(modint lhs, const modint &rhs) { return lhs *= rhs; }
  friend modint operator/(modint lhs, const modint &rhs) { return lhs /= rhs; }
  friend bool operator==(const modint &lhs, const modint &rhs) {
    return lhs.v == rhs.v;
  }
  friend bool operator!=(const modint &lhs, const modint &rhs) {
    return lhs.v != rhs.v;
  }
};
constexpr int mod = 10007;
typedef modint<mod * mod> mint;
template <int N>
struct C_prime {
  mint fac[N + 10], ifac[N + 10];
  C_prime() {
    for (int i = raw(fac[0] = 1); i <= N; i++) fac[i] = fac[i - 1] * i;
    ifac[N] = 1 / fac[N];
    for (int i = N; i >= 1; i--) ifac[i - 1] = ifac[i] * i;
  }
  mint operator()(int n, int m) {
    return n >= m ? fac[n] * ifac[m] * ifac[n - m] : 0;
  }
};
int gcd(int x, int y) { return !y ? x : gcd(y, x % y); }
int n, m, cnt[4];
char cols[30];
C_prime<30> binom;
mint expand(int d, int n) {
  //[x^n] e^{dx}
  return qpow(mint(d), n);
}
mint query(int n) {
  mint ans = 0;
  for (int i = 0; i <= cnt[0]; i++) {
    for (int j = 0; j <= cnt[1]; j++) {
      for (int k = 0; k <= cnt[2]; k++) {
        ans += ((cnt[0] - i + cnt[2] - k) & 1 ? -1 : 1) *
                binom(cnt[0], i) * binom(cnt[1], j) * binom(cnt[2], k) *
                expand(i * 2 - cnt[0] + j * 2 - cnt[1] + k + cnt[3], n);
      }
    }
  }
  return ans;
}
mint Ans[100010];
int mian() {
  cin >> cols;
  memset(cnt, 0, sizeof cnt);
  for (int i = 0; i < n; i++) cnt[cols[i] - 'A']++;
  for (int i = 1; i * i <= m; i++)
    if (m % i == 0) {
      Ans[i] = query(i);
      Ans[m / i] = query(m / i);
    }
  mint ans = 0;
  for (int i = 1; i <= m; i++) {
    int k = gcd(m, i);
    if (cnt[0] && (m / k) % 2 == 0) continue;
    ans += Ans[k];
  }
  ans /= qpow(mint(2), cnt[0] + cnt[1]);
  if (m % mod == 0)
    cout << raw(mint(raw(ans) / mod) / mint(m / mod)) % mod << endl;
  else
    cout << raw(ans / m) % mod << endl;
  return 0;
}
int main() {
#ifndef NF
  cin.tie(nullptr)->sync_with_stdio(0);
#endif
  while (cin >> n >> m, n && m) mian();
  return 0;
}