SDOI2018 旧试题（莫比乌斯反演）（三元环计数）

传送门

旧试题 ---- 约数个数和 ？
利用那道题的推论，发现
$$d(ijk)=\sum _{p|i}\sum _{q|j}\sum _{r|k}[(p,q)=1][(q,r)=1][(p,r)=1]$$
证明类似不再赘述
接下来用你高超的反演技巧可以得到
$$Ans=\sum _{i=1}^{min(A,C)}\sum _{j=1}^{min(A,B)}\sum _{k=1}^{min(B,C)}\mu (i)\mu (j)\mu (k)\sum _{lcm(i,k)|a}\lfloor \frac{A}{a}\rfloor \sum _{lcm(i,j)|b}\lfloor \frac{B}{b}\rfloor \sum _{lcm(j,k)|c}\lfloor \frac{C}{c}\rfloor$$
令 $f(x)={\sum }_{x|a}\lfloor \frac{A}{a}\rfloor$ 就可以在 $nlog$ 的时间内解决后面一坨
发现 $lcm$ 其实把个数限制的很死
我们考虑将 $\mu (i)\ne 0,\mu (j)\ne 0,lcm(i,j)\le max(A,B,C)$ 的那出来，在 $i,j$ 之间建一条边
惊讶的发现边数非常少
那么有贡献的三元组就是一个三圆环，按照三元环计数的方法对其排序统计贡献即可
需要特殊处理右两个相同或者三个相同的情况

#include<bits/stdc++.h>
#define cs const
using namespace std;
cs int N = 1e5 + 5;
cs int Mod = 1e9 + 7;
int read(){
	int cnt = 0, f = 1; char ch = 0;
	while(!isdigit(ch)){ ch = getchar(); if(ch == '-') f = -1; }
	while(isdigit(ch)) cnt = cnt*10 + (ch-'0'), ch = getchar();
	return cnt * f;
}
typedef long long ll;
int gcd(int a, int b){ return !b ? a : gcd(b, a % b); }
int mu[N], prim[N], pc; bool isp[N];
void prework(int n){
	mu[1] = 1; 
	for(int i = 2; i <= n; i++){
		if(!isp[i]) prim[++pc] = i, mu[i] = -1;
		for(int j = 1; j <= pc; j++){
			if(i * prim[j] > n) break; isp[i * prim[j]] = 1;
			if(i % prim[j] == 0) break; mu[i * prim[j]] = -mu[i];
		}
	}
}
int T, A, B, C, n, fa[N], fb[N], fc[N];
struct edge{ int to, lcm; edge(int x=0, int y=0){ to=x, lcm=y; }} ; 
int u[N*10], v[N*10], vl[N*10], tot, deg[N];
vector<edge> E[N];
void Clear(){
	for(int i = 0; i <= n; i++) 
	E[i].clear(), deg[i] = fa[i] = fb[i] = fc[i] = 0;
	tot = 0;
}
bool cmp(int a, int b){ return deg[a]>deg[b]||(deg[a]==deg[b]&&a>b); }
void Solve(){
	A = read(), B = read(), C = read(); 
	n = max(A, max(B, C));
	for(int i = 1; i <= n; i++)
	for(int j = i; j <= n; j += i) fa[i] += A/j, fb[i] += B/j, fc[i] += C/j;
	ll ans = 0;
	for(int i = 1; i <= n; i++) ans += (ll)mu[i]*fa[i]*fb[i]*fc[i];
	for(int i = 1; i <= n; i++) if(mu[i])
	for(int j = 1, lj = n/i; j <= lj; j++) if(mu[i*j])
	for(int k = 1, lk = n/i/j; k <= lk; k++) if(k!=j && mu[i*k] && gcd(j,k) == 1){
		int x = i*j, y = i*k, z = i*j*k, tmp = mu[x]*mu[x]*mu[y];
		ans += (ll)tmp*fa[x]*fb[z]*fc[z];
		ans += (ll)tmp*fa[z]*fb[x]*fc[z];
		ans += (ll)tmp*fa[z]*fb[z]*fc[x];
		if(x>y) ++tot, u[tot]=x, v[tot]=y, vl[tot]=z, ++deg[x], ++deg[y];
	}
	for(int i = 1; i <= tot; i++){
		if(cmp(u[i],v[i])) E[u[i]].push_back(edge(v[i], vl[i]));
		else E[v[i]].push_back(edge(u[i], vl[i]));
	}
	static int lc[N], vis[N]; int TIME = 0;
	for(int x = 1; x <= n; x++) if(mu[x]){
		++TIME;
		for(int i = 0; i < E[x].size(); i++) lc[E[x][i].to] = E[x][i].lcm, vis[E[x][i].to] = TIME;
		for(int e = 0, y, w1; e < E[x].size(); e++){
			y = E[x][e].to; w1 = E[x][e].lcm;
			for(int o = 0, z, w2; o < E[y].size(); o++){
				z = E[y][o].to, w2 = E[y][o].lcm; 
				if(vis[z] ^ TIME) continue;
				int tmp = mu[x]*mu[y]*mu[z];
				ans += (ll)tmp*fa[w1]*fb[w2]*fc[lc[z]];
				ans += (ll)tmp*fa[w1]*fb[lc[z]]*fc[w2];
				ans += (ll)tmp*fa[w2]*fb[w1]*fc[lc[z]];
				ans += (ll)tmp*fa[w2]*fb[lc[z]]*fc[w1];
				ans += (ll)tmp*fa[lc[z]]*fb[w1]*fc[w2];
				ans += (ll)tmp*fa[lc[z]]*fb[w2]*fc[w1];
			}
		}
	} cout << ans % Mod << '\n';
}
int main(){
	T = read(); prework(1e5);
	while(T--) Solve(), Clear();
	return 0;
}