BZOJ 4314 倍数？倍数！

Description

要求你从$Z_n$（模$n$的剩余系）里选出$k$个不重复元素，使得他们的和模$n$为0。求方案数。

$1\leq n\leq 10^9, 1\leq k\leq 10^3$。

Solution

看上去就很像单位根反演的题？

先放一下众所周知的单位根反演的式子： $$\frac{1}{n}\sum_{i = 0}^{n - 1}(\xi_{n}^i)^k=[n | k]$$ 证明很容易：假如有$n | k$，那么和式中每一项都是1，加起来除$n$后就事1了（这种情况下等比数列公比为1，所以不可以采用等比数列求和公式）；反之，我们采用等比数列求和公式可知原式中的和式为$\frac{1 - \xi_n^{kn}}{1 - \xi_n^k}$，这个东西的分子为0，那么自然和式左边为0。

那么考虑构造二元答案多项式$f(x, y)$： $$f(x, y) = \prod_{i = 0}^{n - 1}(1 + x^iy)$$ 这个多项式的每一项中$x$的次数表示和，$y$表示选数的数量。我们显然要求所有含$y^k$且$x$的次数为$n$的倍数的项的系数之和，那么考虑将$y$视为常数，对$x$进行单位根反演： $$f(\xi_n^t, y)=\prod_{i = 0}^{n - 1}(1 + (\xi_n^{t})^iy)$$ 观察到如果$t$和$n$不互质，那么两者可以同时约去一约数（反之，若有$t\perp n$，那么称$\xi_n^t$为一个$n$次本原单位根）。假设约完之后的单位根为$\xi_a^b$，那么我们发现当积式中枚举的$i$大于$a$时，会出现循环。所以我们只需要取循环节（显然长度为$a$，因为我们知道所有$a$阶单位根构成一个循环群，而$a$阶本原单位根一定是该群的生成元。原因事首先显然用$\xi_a^1$可以生成所有$a$阶单位根，所以所有$a$阶单位根构成一个循环群；然后你根据数论里的欧拉定理可以知道$x^{\phi(p)}\equiv 1\pmod{p}(x\perp p)$，那么对于任意本原单位根$\xi_a^b$有$\xi_a^{b^{\phi(a)}} = \xi_a^1$，式子左边的东西是$\xi_a^b$的若干次方，所以用$\xi_a^b$可以生成$\xi_a^1$，自然就可以生成整个循环群）的若干次方就行了，那么有： $$f(\xi_a^b, y) = (\prod_{i = 0}^{a - 1}(1 + (\xi_a^b)^iy))^\frac{n}{a}$$

根据上面的说法，考虑里面的积式一定可以枚举到所有$a$阶单位根，所以这个积式也可以写成： $$f(\xi_a^b, y) = (\prod_{i = 0}^{a - 1} (1 + \xi_a^iy))^\frac{n}{a}$$ 这么一来，我们发现对于任意$a$阶本原单位根$\xi_a^b$，对答案的贡献都是一样的，而我们知道$a$阶本原单位根有$\phi(a)$个。只要我们对所有$n$的约数求出其所有本原单位根的贡献的和，那么这题就做完了。

那么我们考虑，如果我们知道了$a$，那么怎么快速求一个本原单位根的贡献呢？

考虑那个多项式（里面的积式），首先很容易发现其0次项为1。然后我们从零点入手……

然后我们发现这个多项式几乎没法求零点……那么我们做一步小转换吧： $$\prod_{i = 0}^{a - 1}\xi_a^i(\xi_a^{a - i} + y)$$

然后我们发现$\xi_a^{a - i}$就枚举了所有单位根，所以零点集合就是所有单位根相反数的集合。

如果$a$为偶数，那么把所有单位根取相反数之后得到的集合其实和原集合事一样的……这个大致可以理解为所有单位根构成的图形关于原点中心对称。而我们知道$x^a - 1 = 0$的解就是全体$a$阶单位根，因此答案多项式和$x^a - 1$只有常数项的不同，而取反之后得到了$1 - x^a$，其零点、常数项都符合我们的要求。

如果$a$为奇数，那么考虑$x^a$，往里面带一个单位根的相反数的话，因为$a$为奇数所以会得到$-1$，所以该多项式的零点就是$x^a + 1 = 0$的全体解。而左边多项式的常数项也符合我们的要求，因此此时多项式为$1 + x^a$。

综上可得： $$\prod_{i = 0}^{a - 1}(1 + \xi_a^iy) = 1 - (-y)^a$$ 而最后的答案多项式为： $$\frac{1}{n}\sum_{d | n}\phi(d)(1 - (-y)^d)^\frac{n}{d}$$ 这个东西的$y^k$项的系数就是答案了……

Code

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <functional>
#include <utility>
typedef long long ll;
const ll ha = 1000000007LL;
ll pow_mod(ll a, ll b) {
  ll ans = 1LL, res = a;
  while(b) {
    if(1LL & b) ans = (ans * res) % ha;
    res = (res * res) % ha; b >>= 1;
  }
  return ans;
}
ll inv(ll x) {
  return pow_mod(x, ha - 2LL);
}
 
ll sqrt(ll n) {
  ll ans = 1;
  while((ans + 1LL) * (ans + 1LL) <= n) ans ++;
  return ans;
}
int phi(int n) {
  int ans = n, m = sqrt(n);
  for(int i = 2; i <= m; i ++) {
    if(n % i == 0) {
      ans = (ans / i) * (i - 1);
      while(n % i == 0) n /= i;
    }
  }
  if(n > 1) ans = (ans / n) * (n - 1);
  return ans;
}
ll C(int n, int m) {
  if(n < m) return 0;
  ll ans = 1;
  for(int i = 1; i <= m; i ++) {
    ans = (ans * (ll(n - i + 1))) % ha;
    ans = (ans * inv(i)) % ha;
  }
  return ans;
}
ll query(int d, int nd, int k) {
  if(k % d != 0) return 0;
  ll ret = phi(d);
  if(d % 2 == 0 && (k / d) % 2 != 0) {
    ret = (ha - ret) % ha;
  }
  ret = (ret * C(nd, k / d)) % ha;
  return ret;
}
ll calc(int n, int k) {
  int m = sqrt(n); ll ans = 0;
  for(int i = 1; i <= m; i ++) {
    if(n % i != 0) continue;
    ans = (ans + query(i, n / i, k)) % ha;
    if(i * i != n) ans = (ans + query(n / i, i, k)) % ha;
  }
  ans = (ans * inv(n)) % ha;
  return ans;
}
 
int main() {
  int n, k; scanf("%d%d", &n, &k);
  printf("%lld\n", calc(n, k));
  return 0;
}