[莫比乌斯反演][三元环计数]「SDOI2018」旧试题
-
一个结论：\(d(ijk)=\sum_{a|i,b|j,c|k}[(a,b)=(a,c)=(b,c)=1]\)
-
证明：分每个质因子独立考虑，因为 \(a,b,c\) 两两互质，故对于任意一个质因子 \(p\) 这三个数中最多存在一个数出现过 \(p\) 这个质因子，这样一个质因子 \(p\) 造成的乘积贡献为 \(i,j,k\) 中 \(p\) 的个数之和再加 \(1\)，也就是 \(ijk\) 中 \(p\) 的个数之和再加 \(1\)
-
于是开始反演
-
\[\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{a|i,b|j,c|k}\sum_{x|a,x|b}\sum_{y|a,y|c}\sum_{z|b,z|c}\mu(x)\mu(y)\mu(z) \]
-
\[=\sum_{x,y,z}\mu(x)\mu(y)\mu(z)\sum_{lcm(x,y)|a}\sum_{lcm(x,z)|b}\sum_{lcm(y,z)|c}\lfloor\frac Aa\rfloor\lfloor\frac Bb\rfloor\lfloor\frac Cc\rfloor \]
-
设 \(f(n,x)=\sum_{x|i}\lfloor\frac ni\rfloor\)（对于一个 \(n\) 可以 \(O(n\log n)\) 预处理出所有 \(f(n,x)\)）：
-
\[\sum_{x,y,z}\mu(x)\mu(y)\mu(z)f(A,lcm(x,y))f(B,lcm(x,z))f(C,lcm(y,z)) \]
-
这看上去还是三方的，但我们注意到 \(f(n,>n)=0\)，而这个式子中 \(f\) 的第二个参数又是 \(lcm\)，我们可以直观感觉到贡献不为 \(0\) 的 \((x,y,z)\) 是很少的
-
由于 \((x,y)(x,z)(y,z)\) 同时出现，考虑建一个图，点集为 \([1,\max(A,B,C)]\) 中所有 \(\mu\) 不为 \(0\) 的数，如果 \(lcm(x,y)\le\max(A,B,C)\) 就连边 \((x,y)\)，我们要考虑的就是这个图所有的三元环，注意自环的特殊处理
-
先考虑如何建图：设 \(n=\max(A,B,C)\)，先枚举 \(d=(x,y)\) 再枚举 \(i=\frac xd,j=\frac yd\)，满足 \(\frac{xy}d=ijd\le n\) 也就是 \(ij\le\frac nd\)，如果 \((i,j)=1\) 就连边 \((id,jd)\)
-
这样的复杂度为 \(\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac nd\rfloor}\lfloor\frac n{id}\rfloor=\sum_{d=1}^nO(\frac nd\log\frac nd)<O(n\log^2n)\)。
-
建完这张图之后发现 \(n=10^5\) 时边数只有 \(760741\)，于是可以枚举这张图的三元环
-
三元环计数黑科技：把所有的点按照度数从大到小排序，度数相同情况下任意。然后把所有边定向，排序后在前面（度数大）的点（连向在后面（度数小）的点，形成一个 DAG
-
枚举三元环上排在最前面的点 \(u\)，把 \(u\) 的所有出点标记，然后枚举 \(u\) 的出边 \(u\rightarrow v\) 和 \(v\) 的出边 \(v\rightarrow w\)，如果 \(w\) 被标记则找到了一个三元环
-
复杂度为 \(O(m\sqrt m)\)，\(m\) 为边数。分析（\(d_u\) 为点 \(u\) 的度数，\(cnt_u\) 为点 \(u\) 的出度数）：
-
复杂度相当于对每条边的出点 \(u\) 求 \(cnt_u\) 的和
-
若一条边的出点 \(u\) 满足 \(d_u\le\sqrt m\)，则 \(cnt_u\) 必然也不超过 \(\sqrt m\)，这一部分的复杂度 \(O(m\sqrt m)\)
-
否则设入点为 \(v\)，则必然 \(d_v\ge d_u>sqrt_m\)，由于图中 \(d>\sqrt m\) 的点不超过 \(\sqrt m\) 个，故每个 \(cnt_u\) 贡献不超过 \(\sqrt m\) 次，又因为所有点的 \(cnt\) 之和为 \(m\)，故这一部分复杂度也是 \(O(m\sqrt m)\)
Code
#include <bits/stdc++.h>
template <class T>
inline void read(T &res)
{
res = 0; bool bo = 0; char c;
while (((c = getchar()) < '0' || c > '9') && c != '-');
if (c == '-') bo = 1; else res = c - 48;
while ((c = getchar()) >= '0' && c <= '9')
res = (res << 3) + (res << 1) + (c - 48);
if (bo) res = ~res + 1;
}
typedef long long ll;
const int N = 1e5 + 5, M = 2e6 + 5, djq = 1e9 + 7;
int n, A, B, C, tot, pri[N], miu[N], ecnt, nxt[M], adj[N], go[M], val[M], d[N],
fA[N], fB[N], fC[N], X[M], Y[M], Z[M], ans, m, a[N], id[N], vis[N];
bool mark[N];
void sieve()
{
mark[0] = mark[miu[1] = 1] = 1;
for (int i = 2; i <= 100000; i++)
{
if (!mark[i]) miu[pri[++tot] = i] = -1;
for (int j = 1; j <= tot; j++)
{
if (1ll * i * pri[j] > 100000) break;
mark[i * pri[j]] = 1;
if (i % pri[j] == 0) break;
else miu[i * pri[j]] = -miu[i];
}
}
}
void add_edge(int u, int v, int w)
{
nxt[++ecnt] = adj[u]; adj[u] = ecnt; go[ecnt] = v; val[ecnt] = w;
}
inline bool comp(int a, int b) {return d[a] > d[b];}
void work()
{
read(A); read(B); read(C);
n = std::max(std::max(A, B), C);
ecnt = m = ans = 0; int ec = 0;
for (int i = 1; i <= n; i++) adj[i] = 0;
for (int d = 1; d <= n; d++)
for (int i = 1; i <= n / d; i++) if (miu[i * d])
for (int j = 1; j <= n / d / i; j++)
if (miu[j * d] && std::__gcd(i, j) == 1)
X[++ec] = i * d, Y[ec] = j * d, Z[ec] = i * j * d;
for (int i = 1; i <= n; i++)
{
fA[i] = fB[i] = fC[i] = d[i] = 0;
for (int j = i; j <= n; j += i)
fA[i] += A / j, fB[i] += B / j, fC[i] += C / j;
if (miu[i]) a[++m] = i;
}
std::sort(a + 1, a + m + 1);
for (int i = 1; i <= m; i++) id[a[i]] = i;
for (int i = 1; i <= ec; i++) if (id[X[i]] < id[Y[i]])
add_edge(X[i], Y[i], Z[i]);
for (int i = 1; i <= m; i++)
{
int u = a[i];
for (int e = adj[u], v = go[e]; e; e = nxt[e], v = go[e])
vis[v] = val[e];
for (int e = adj[u], v = go[e]; e; e = nxt[e], v = go[e])
{
ans = (ans + 1ll * fA[u] * fB[val[e]] % djq * fC[val[e]] % djq
* miu[v] + djq) % djq;
ans = (ans + 1ll * fB[u] * fA[val[e]] % djq * fC[val[e]] % djq
* miu[v] + djq) % djq;
ans = (ans + 1ll * fC[u] * fA[val[e]] % djq * fB[val[e]] % djq
* miu[v] + djq) % djq;
ans = (ans + 1ll * fA[v] * fB[val[e]] % djq * fC[val[e]] % djq
* miu[u] + djq) % djq;
ans = (ans + 1ll * fB[v] * fA[val[e]] % djq * fC[val[e]] % djq
* miu[u] + djq) % djq;
ans = (ans + 1ll * fC[v] * fA[val[e]] % djq * fB[val[e]] % djq
* miu[u] + djq) % djq;
for (int x = adj[v], w = go[x]; x; x = nxt[x], w = go[x])
if (vis[w])
{
ll sum = 0;
sum += 1ll * fA[val[e]] * fB[val[x]] % djq * fC[vis[w]];
sum += 1ll * fA[val[e]] * fC[val[x]] % djq * fB[vis[w]];
sum += 1ll * fB[val[e]] * fA[val[x]] % djq * fC[vis[w]];
sum += 1ll * fB[val[e]] * fC[val[x]] % djq * fA[vis[w]];
sum += 1ll * fC[val[e]] * fA[val[x]] % djq * fB[vis[w]];
sum += 1ll * fC[val[e]] * fB[val[x]] % djq * fA[vis[w]];
if (sum %= djq, miu[u] * miu[v] * miu[w] == 1)
ans = (ans + sum) % djq;
else ans = (ans - sum + djq) % djq;
}
}
for (int e = adj[u], v = go[e]; e; e = nxt[e], v = go[e])
vis[v] = 0;
if (miu[u] == 1) ans = (1ll * fA[u] * fB[u] % djq * fC[u] + ans) % djq;
else ans = (ans - 1ll * fA[u] * fB[u] % djq * fC[u] % djq + djq) % djq;
}
printf("%d\n", ans);
}
int main()
{
sieve();
int T; read(T);
while (T--) work();
return 0;
}
