P4619 [SDOI2018]旧试题

P4619 [SDOI2018]旧试题

题目描述

求下面表达式的值

\[\large\sum_{i=1}^{A}\sum_{j=1}^{B}\sum_{k=1}^{C}\sigma_{0}(ijk) \mod(10^9+7) \]

数据范围

\(1 \le A,B,C \le 100000\)

解题思路

这题非常困难是因为推着推着式子变成了图论问题？

让我们先来推一波式子把

要用到的结论

\[\large \sigma_{0}(ijk)=\sum_{x|i}\sum_{y|j}\sum_{z|k}[\gcd(x,y)=1][\gcd(y,z)=1][\gcd(x,z)=1]\\ [x=1] = \sum_{d|x}\mu(d) \]

第一个式子的理解

单独考虑一个质因子 P，i 中最大有 \(P^a\)，j 中最大有 \(P^b\)，k 中最大有 \(P^c\)，那么 (ijk) 中有 \(P^{a+b+c}\)， 现在我们尝试建立起映射，选因数时，如果 P 选了小于等于 a 个，那么全在 i 中选，显然满足全部互质的条件，大于 a 个但小于等于 a + b 个，那么优先在 i 中选，剩下的在 j 中选，然后将 i 清零，大于 a + b 也同理，也就是选质数时尽量靠前选，但我们只让最后一个被选到的有，前面的设为 0，这样可以保证互质，且是一种双射关系

第二个式子的理解

裸的莫比乌斯反演

正式开始推式子

\[\large\sum_{i=1}^{A}\sum_{j=1}^{B}\sum_{k=1}^{C}\sigma_{0}(ijk)\\ \large=\sum_{i=1}^{A}\sum_{j=1}^{B}\sum_{k=1}^{C}\sum_{x|i}\sum_{y|j}\sum_{z|k}[\gcd(x,y)=1][\gcd(y,z)=1][\gcd(x,z)=1]\\ \large=\sum_{x=1}^{A}\sum_{y=1}^{B}\sum_{z=1}^{C}[\gcd(x,y)=1][\gcd(y,z)=1][\gcd(x,z)=1]\sum_{i=1}^{\frac Ax}\sum_{j=1}^{\frac By}\sum_{k=1}^{\frac Cz}\\ \large=\sum_{x=1}^{A}\sum_{y=1}^{B}\sum_{z=1}^{C}[\gcd(x,y)=1][\gcd(y,z)=1][\gcd(x,z)=1]\lfloor\frac Ax \rfloor\lfloor\frac By \rfloor\lfloor\frac Cz \rfloor\\ \large=\sum_{x=1}^{A}\sum_{y=1}^{B}\sum_{z=1}^{C}\sum_{u|x,u|y}\mu(u)\sum_{v|y,v|z}\mu(v)\sum_{w|x,w|z}\mu(w)\lfloor\frac Ax \rfloor\lfloor\frac By \rfloor\lfloor\frac Cz \rfloor\\ \large=\sum_{u=1}^{\min(A,B)}\sum_{v=1}^{\min(B,C)}\sum_{w=1}^{\min(A,C)}\mu(u)\mu(v)\mu(w)\sum_{x=lcm(u,w)}\lfloor\frac Ax \rfloor\sum_{y=lcm(u,v)}\lfloor\frac By\rfloor\sum_{z=lcm(v,w)} \lfloor\frac Cz \rfloor\\ \large=\sum_{u=1}^{\min(A,B)}\sum_{v=1}^{\min(B,C)}\sum_{w=1}^{\min(A,C)}\mu(u)\mu(v)\mu(w)F(lcm(u,w),A)F(lcm(u,v),B)F(lcm(v,w),C) \]

然后就是超级神奇的操作

我们把枚举的三元组抽象为一个三元环 (x, y, z)，两个点间连边权为 lcm 的边，这个三元环对答案有贡献当且仅当

• \(\mu(x)\neq0 ~~\mu(y) \neq0~\mu(z)\neq0\)
• \(lcm(x, y, z) \le \max(A, B, C)\)

你觉得它可能会快，但你认为它还是暴力，事实上它跑的飞快（比纯暴力来说

枚举三元环肯定不是暴力枚举，容易发现 lcm 的 \(\mu\) 也不为零，我们可以直接枚举 lcm，将其质因数分解，根据质因数枚举其中一个端点，然后暴力枚举另一个端点，这样枚举的每一个都是合法的点对，时间复杂度 \(\Theta(m)\)

打表发现有大概 70 多万条边，暴力枚举肯定不行，用上黑科技 三元环计数 可以做到 \(\Theta(m\sqrt m)\)，卡卡常勉强跑过。

一发过了，但除法真的比乘法慢一大堆，请谨慎使用除法！！

#pragma GCC optimize(3, "inline")
#include <queue>
#include <vector>
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define MP make_pair
#define ll long long
#define fi first
#define se second
using namespace std;

template <typename T>
void read(T &x) {
    x = 0; bool f = 0;
    char c = getchar();
    for (;!isdigit(c);c=getchar()) if (c=='-') f=1;
    for (;isdigit(c);c=getchar()) x=x*10+(c^48);
    if (f) x=-x;
}

template<typename F>
inline void write(F x)
{
	static short st[30];short tp=0;
	if(x<0) putchar('-'),x=-x;
	do st[++tp]=x%10,x/=10; while(x);
	while(tp) putchar('0'|st[tp--]);
	putchar('\n');
}

template <typename T>
inline void Mx(T &x, T y) { x < y && (x = y); }

template <typename T>
inline void Mn(T &x, T y) { x > y && (x = y); }

#define re register
const int Maxn = 300000, N = Maxn + 5;
const int P = 1e9+7;
ll f[N], mu[N], e[N], prime[N], p[N], tot;
void prework(void) {
	mu[1] = 1; f[1] = 1;
	for (re int i = 2;i <= Maxn; i++) {
		if (!e[i]) prime[++tot] = e[i] = i, mu[i] = -1, f[i] = p[i] = 2;
		for (re int j = 1; prime[j] * i <= Maxn && j <= tot; j++) {
			int t = prime[j] * i; e[t] = prime[j];
			if (prime[j] == e[i]) { p[t] = p[i] + 1, f[t] = f[i] / p[i] * p[t]; break; } 
			f[t] = f[i] * (p[t] = 2), mu[t] = mu[i] * -1;
		}
//		cout << f[i] << endl;
		f[i] = f[i] + f[i-1]; if (f[i] >= P) f[i] -= P;
	}
}


struct node {
	int x, y, lcm;
}ed[N*20];

ll ans;
ll vis[N], lm[N], a[N], b[N], c[N], A, B, C;
ll st[N], deg[N], lim, tp, T;
vector<pair<ll, ll> > v[N];
inline void clear(void) {
	lim = max(max(A, B), C), tot = ans = 0;
	for (int i = 1;i <= lim; i++) 
		a[i] = f[A/i], b[i] = f[B/i], c[i] = f[C/i], v[i].clear();
}

inline ll calc(int x, int y, int z) {
	return a[x] * b[y] * c[z];
}

inline void build(void) {
//	cout << calc(2, 1, 2) + calc(2, 2, 1) + calc(1, 2, 2) << endl;
	memset(deg, 0, lim * 4 + 200);
	for (int i = 1;i <= lim; i++) {
		if (!mu[i]) continue; int x = i; tp = 0;
//		if (i % 1000 == 0) cout << i << endl;
		while (x != 1) st[++tp] = e[x], x /= e[x];
		for (re int j = 0; j < (1 << tp); j++) {
			int u = 1, t, v;
			for (re int k = 1;k <= tp; k++)
				if (j & (1 << (k - 1))) u *= st[k];
			t = i / u;
			for (int s = j; ; s = (s - 1) & j) {
				v = t;
				for (int p = 0;p < tp; p++)
					if (s >> p & 1) v *= st[p + 1];
//				cout << u << ' ' << v << endl;
				if (u > v) {
					ed[++tot] = (node) {u, v, i}; deg[u]++, deg[v]++;
//					cout << u << ' ' << i << endl;
//					cout << calc(u, i, i) + calc(i, u, i) + calc(i, i, u) << endl;
					ans += mu[v] * (calc(u, i, i) + calc(i, u, i) + calc(i, i, u));
					ans += mu[u] * (calc(v, i, i) + calc(i, v, i) + calc(i, i, v));
				}
				if (!s) break;
			}
		}
		ans += mu[i] * calc(i, i, i);
	}
	for (int i = 1;i <= tot; i++) {
		int x = ed[i].x, y = ed[i].y;
		if (deg[x] > deg[y] || (deg[x] == deg[y] && x > y)) swap(x, y);
		v[x].emplace_back(y, ed[i].lcm);
	}
}

void work(void) {
	memset(vis, 0, lim * 4 + 200);
//	cout << calc(10, 5, 2) + calc(10, 2, 5) + calc(2, 5, 10) + calc(2, 10, 5) << endl;
	for (int i = 1;i <= lim; i++) {
		for (pair<int, int> to: v[i]) 
			vis[to.fi] = i, lm[to.fi] = to.se;
		for (pair<int, int> to: v[i]) {
			for (pair<int, int> To: v[to.fi]) {
				if (vis[To.fi] != i) continue;
//				int pre = ans;
//				cout << lm[To.fi] << ' ' << to.se << ' ' << To.se << endl;
				ans += mu[i] * mu[to.fi] * mu[To.fi] * 
					( calc(lm[To.fi], to.se, To.se) 
					+ calc(lm[To.fi], To.se, to.se)
					+ calc(to.se, lm[To.fi], To.se) 
					+ calc(to.se, To.se, lm[To.fi])
					+ calc(To.se, lm[To.fi], to.se)
					+ calc(To.se, to.se, lm[To.fi]) );
//				cout << ans - pre << endl;
			}
		}
		ans %= P;
	}
	ans = (ans % P + P) % P;
}

int main() {
	for (prework(), read(T); T; T--) {
		read(A), read(B), read(C);
		clear(), build(), work(), write(ans);
	}
	return 0;
}
/*

1
100000 100000 100000

*/