[莫比乌斯反演][三元环计数]「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}^Asum_{j=1}^Bsum_{k=1}^Csum_{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}lfloorfrac Aa floorlfloorfrac Bb floorlfloorfrac Cc floor ]

• 设 (f(n,x)=sum_{x|i}lfloorfrac ni floor)（对于一个 (n) 可以 (O(nlog 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)lemax(A,B,C)) 就连边 ((x,y))，我们要考虑的就是这个图所有的三元环，注意自环的特殊处理

• 先考虑如何建图：设 (n=max(A,B,C))，先枚举 (d=(x,y)) 再枚举 (i=frac xd,j=frac yd)，满足 (frac{xy}d=ijdle n) 也就是 (ijlefrac nd)，如果 ((i,j)=1) 就连边 ((id,jd))

• 这样的复杂度为 (sum_{d=1}^nsum_{i=1}^{lfloorfrac nd floor}lfloorfrac n{id} floor=sum_{d=1}^nO(frac ndlogfrac nd)<O(nlog^2n))。

• 建完这张图之后发现 (n=10^5) 时边数只有 (760741)，于是可以枚举这张图的三元环

• 三元环计数黑科技：把所有的点按照度数从大到小排序，度数相同情况下任意。然后把所有边定向，排序后在前面（度数大）的点（连向在后面（度数小）的点，形成一个 DAG

• 枚举三元环上排在最前面的点 (u)，把 (u) 的所有出点标记，然后枚举 (u) 的出边 (u ightarrow v) 和 (v) 的出边 (v ightarrow w)，如果 (w) 被标记则找到了一个三元环

• 复杂度为 (O(msqrt m))，(m) 为边数。分析（(d_u) 为点 (u) 的度数，(cnt_u) 为点 (u) 的出度数）：

• 复杂度相当于对每条边的出点 (u) 求 (cnt_u) 的和

• 若一条边的出点 (u) 满足 (d_ulesqrt m)，则 (cnt_u) 必然也不超过 (sqrt m)，这一部分的复杂度 (O(msqrt m))

• 否则设入点为 (v)，则必然 (d_vge d_u>sqrt_m)，由于图中 (d>sqrt m) 的点不超过 (sqrt m) 个，故每个 (cnt_u) 贡献不超过 (sqrt m) 次，又因为所有点的 (cnt) 之和为 (m)，故这一部分复杂度也是 (O(msqrt 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
", ans);
}

int main()
{
	sieve();
	int T; read(T);
	while (T--) work();
	return 0;
}