题解 P1587 【[NOI2016]循环之美】
知识点：莫比乌斯反演 积性函数 杜教筛
废话前言：
我是古明地恋，写这篇题解的人已经被我
请各位读者自行无视搞事的恋恋带有删除线的内容，谢谢茄子。
这道题目本身并不难，但是公式推导/代码过程中具有迷惑性的内容比较多，因此，像我这种什么都不会还眼瞎的马上就该退役回归文化课的辣鸡废柴蒟蒻/kel/kel/kk就会在debug上浪费无谓地消耗大量时间。
因此，在平时写代码时，养成良好的代码习惯（比如：确定常用的变量名，不轻易修改；即使是十分熟悉的函数，也认真写完每一行等）是非常重要的。考试时省下的10分钟debug时间，或许就决定了你没有大学上个人的竞赛生涯是否能继续下去。
此外，我写这篇题解的目的并不是简单地告诉大家如何AC这道水题，而是将我的思考过程分享给大家，让大家下一次遇到类似的数论题目的时候能够真正做到独立思考。因此大家不要过多依赖题解中的代码现成的公式，而是利用纸笔自行推演。
PS:第一次用\(L_AT^EX\)写公式，真好用qwq
那么，开始正题：
注意：此处默认大家已经对莫比乌斯函数等常见的数论函数有了一定的掌握，不再赘述。这你都不会还学数论？
一句话题意：求\(k\)进制下，满足\(1 \leq x \leq n, 1 \leq y \leq m,\)且 \(\frac{x}{y}\)是纯循环小数(或整数)的既约分数的二元组\((x,y)\)的个数。
那么很显然打表或者自行探索规律发现，当分母\(y\)与\(k\)互质时，\(\frac{x}{y}\)是纯循环小数。而既约分数需满足分子\(x\)与分母\(y\)互质。
由此我们得到下面的柿子：
然而直接枚举的复杂度无法接受。
观察到柿子中出现了\([\gcd=1]\)的形式，这启发我们使用莫比乌斯反演的一个常用的变形:
引入到上式当中，得到：
考虑把枚举\(t\)提取到枚举\(j\)的前面，将较复杂的枚举约数转化为简便的枚举倍数。得到（此处枚举的\(j\)实际上是\("t\)的\(j\)倍\("\)）:
即
注意到此处枚举\(t\)的时候，仍然有\(t|i\)的限制。为了将这个限制消除，我们试图将枚举\(t\)转移到更前面的位置，同样把枚举\(i\)转化为\("t\)的\(i\)倍\("\)：
(由于枚举\(i\)转化为了\("t\)的\(i\)倍\("\)，后面的\(d|i\)变成了\(d|it\))
解决完了\(t\)的问题，枚举\(d\)的问题又该怎么办呢？
管那么多干什么 往前提不就行了ww
我们发现，如果确定了\(t\)和\(d\)的值，那么合法的\(i\)的数目是可以\(O(1)\)计算的，并且\(d|k\)的条件十分容易满足。那么我们就考虑把枚举\(d\)提到前面来，把\(d|it\)的限制条件交给\(i\)来完成：
设\(w=\gcd(t,d)\)，则:
此时我们已经得到了一个比较优秀的暴力做法：由于\(k\leq2000\) ，使得\(\mu(d)\not=0\)的\(d\)最多只有\(16\)个，线性筛\(\mu\)，暴力枚举\(d\)和\(t\)，时间复杂度约为\(O(16\min(n,m))\)，期望得分\(70\)左右。
思考一下，暴力为什么会被限制在线性时间呢？
仍然是\(d|it\)的限制，让我们无法对\(t\)进行整除分块。
仔细观察暴力柿子，再撕烤一下，我们发现：由于\(w=\gcd(t,d)\)，如果确定了\(w\)的值，\(i\)的限制就与\(\frac{t}{w}\)无关了。
那么枚举\(\frac{t}{w}\)整除分块？
诶等等，那个奇怪的\([\gcd(t,d)=1]\)是哪来的？
老套路了，此处枚举的\(t\)(我们叫它\(t'\)吧qwq)，实际上是\("w\)的\(t'\)倍\("\)。又因为原式中\(w=\gcd(t,d)\)，所以只有\(t'\perp d\)，才不会重复枚举。
那怎么办？再拆一个\(\mu\)出来？
[本能]本我的解放(物理).png 我觉得不行。
继续观察上式，我们发现，本来\(\sum\mu(t)\)就需要杜教筛，那能不能顺便把\(\gcd\)的限制一起解决呢？
想得美
可以。
于是某个蒟蒻发明了这个函数：
把这个函数代入上式：
说了这么多废话，这不和上面的柿子一模一样？
然而并不是。如果\(d\)是一个给定的常数，那么\(\omega(d,x)\)是积性的。积性，意味着可以杜教筛。
那就来吧：
设\(g(x)\)是积性函数，要求\(\sum_{i=1}^xg(x)\)和\(\sum_{i=1}^x(g*\omega)(x)\)可以快速求解，自然想到\(g=I\)。(此处由于把\(d\)当做了常数，所以\(\omega\)中省略了这个参数)
然后——我们遇到了前所未有的巨大挑战 （其实就是我菜） ：
\(\sum_{i=1}^x(I*\omega)(x)\)的值究竟是多少？
凭着以往的经验，我自信的写下了
int ret = 1;
然后获得了\(64\)分的好成绩。所有需要杜教筛的点无一幸免。
这...
冷静思考一下，我们发现\(I*\omega\not=\epsilon\)。
对于\(I*\mu\)来说，任何\(x\)只要拥有至少一个质因子，就可以用二项式定理证明\((I*\mu)(x)=0\)。
但是\(I*\omega\)呢？由于引入了\(\gcd(d,x)\)的限制，只有当\(x\)拥有至少一个与\(d\)互质的质因子时，才可以得到\((I*\omega)(x)=0\)。因此我们每次杜教筛求\(\omega(d,x)\)的时候，需要用一个\(dfs\)得到满足\(1\leq i \leq x, \gcd(i,d)=1\)的\(i\)的个数，才能正常杜教筛。
那么枚举\(d\)，\(w\)，对\(t\)整除分块，杜教筛\(\omega(d, t)\)，就可以在\(O(3^{h(k)}\sqrt[\frac{2}{3}]{\min(n,m)})\)的时间复杂度上界(\(h(k)\)指\(k\)的质因子数量，并且实际运行远远达不到这个上界)中完成了。
那么，如果不考虑写错变量名调试的一个小时这道题就顺利的解决啦~☆
附上代码：
#include <cstdio>
#include <cctype>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <unordered_map>
#define R register
#define ll long long
using namespace std;
const int N = 1e6;
int n, m, k, omg[16][N + 10], ispr[N + 10], prime[N], nopr, num, p[5], somg[16][N + 10], miu[20], d[20];
unordered_map<int, int>sum[16];
inline void pre() {
int tp = k;
for (R int i = 2; i * i <= k; ++i) {
if (tp % i == 0) {
p[++num] = i;
while (tp % i == 0)
tp /= i;
}
}
if (tp != 1)
p[++num] = tp;
for (R int i = 0; i < (1 << num); ++i) {
d[i] = 1, miu[i] = 1;
for (R int j = 0; j < num; ++j)
if (i & (1 << j))
d[i] *= p[j + 1], miu[i] *= -1;
}
for (R int i = 0; i < (1 << num); ++i)
somg[i][1] = 1;
for (R int i = 2; i <= N; ++i) {
if (!ispr[i]) {
for (R int j = 0; j < (1 << num); ++j)
if (d[j] % i)
omg[j][i] = -1;
prime[++nopr] = i;
}
for (R int j = 0; j < (1 << num); ++j)
somg[j][i] = somg[j][i - 1] + omg[j][i];
for (R int j = 1, t; j <= nopr && (t = prime[j] * i) <= N; ++j) {
ispr[t] = 1;
if (i % prime[j] == 0) {
for (R int u = 0; u < (1 << num); ++u)
omg[u][t] = 0;
break;
}
for (R int u = 0; u < (1 << num); ++u)
omg[u][t] = omg[u][i] * (d[u] % prime[j] ? -1 : 0);
}
}
return;
}
int dfs(int u, int step, int mul, int lim) {
if (step == num)
return 1;
if (~u & (1 << step))
return dfs(u, step + 1, mul, lim);
int ret = 0;
for (R ll i = 1, j; (j = i * mul) <= lim; i *= p[step + 1])
ret += dfs(u, step + 1, j, lim);
return ret;
}
int calc(int u, int x) {
if (x <= N)
return somg[u][x];
if (sum[u].find(x) != sum[u].end())
return sum[u][x];
int ret = dfs(u, 0, 1, x), i;
for (i = 2; i * i <= x; ++i)
ret -= calc(u, x / i);
for (R int k = x / i, j; i <= x; i = j + 1, --k)
j = x / k, ret -= (j - i + 1) * calc(u, k);
return sum[u][x] = ret;
}
int main() {
cin >> n >> m >> k;
pre();
ll ans = 0;
for (R int i = 0; i < (1 << num); ++i) {
for (R int j = 0; j < (1 << num); ++j) {
if (d[i] % d[j])
continue;
ll nx = n / d[j], mx = m / d[i];
if (nx > mx)
swap(nx, mx);
for (R int u = 1, v; u <= nx; u = v + 1) {
v = min(nx / (nx / u), mx / (mx / u));
ans += miu[i] * miu[j] * (calc(i, v) - calc(i, u - 1)) * (nx / u) * (mx / u);
}
}
}
cout << ans;
return 0;//呐，想着复制代码的大哥哥，恋恋就在你身后哟
}
浙公网安备 33010602011771号