【香蕉OI】GCD 和 LCM （莫比乌斯反演）

题意

给出 $T$ 组询问，每组询问求 $$\sum _{i\le n}\sum _{j\le m}[gcd(i,j)\le a]lcm(i,j)$$

$T\le 10^4,n,m,a\le 10^5$

思路

比赛的送分题，看起来很套路，但是我套路都忘光了。

首先枚举 $gcd(i,j)$

$$ans(n,m,a)=\sum _{d\le a}\sum _{i\le \lfloor \frac{n}{d}\rfloor }\sum _{j\le \lfloor \frac{m}{d}\rfloor }[gcd(i,j)=1]\cdot i\cdot j\cdot d$$

套路转化 $[gcd(i,j)=1]={\sum }_{p|i,p|j}\mu (p)$

$$ans(n,m,a)=\sum _{d\le a}\sum _{i\le \lfloor \frac{n}{d}\rfloor }\sum _{j\le \lfloor \frac{m}{d}\rfloor }\sum _{p|i,p|j}\mu (p)\cdot i\cdot j\cdot d$$

换个顺序

$$ans(n,m,a)=\sum _{d\le a}d\sum _{p\le \lfloor \frac{n}{d}\rfloor }\mu (p)\cdot p^2\sum _{i\le \lfloor \frac{n}{dp}\rfloor }i\sum _{j\le \lfloor \frac{m}{dp}\rfloor }j$$

套路转化，令 $k=d\cdot p$，枚举 $k$

$$\begin{array}{rl}ans(n,m,a)& =\sum _{k\le n}\sum _{i\le \lfloor \frac{n}{k}\rfloor }i\sum _{j\le \lfloor \frac{m}{k}\rfloor }j\sum _{d\le a,d|k}d\cdot \mu (\frac{k}{d})\cdot (\frac{k}{d}{)}^2\\ & =\sum _{k\le n}{k}^2\sum _{i\le \lfloor \frac{n}{k}\rfloor }i\sum _{j\le \lfloor \frac{m}{k}\rfloor }j\sum _{d\le a,d|k}\frac{\mu (\frac{k}{d})}{d}\end{array}$$

假如忽略 $d\le a$ 的限制，后面那个部分是一个积性函数，可以预处理。

但是加上了限制之后就不行了。那么我们考虑把询问离线，按照 $a$ 排序，用树状数组维护后面的部分。每次将 $d=a$ 暴力插入树状数组。修改复杂度 $O(a\log a\log n)$。查询的时候整除分块一下。

总复杂度 $O(a\log a\log n+T\sqrt{n}\log n)$

代码

#include<bits/stdc++.h>
using namespace std;
const int TN = 1e4 + 10, N = 1e5 + 10, mod = 1e9 + 7;
int T, ans[TN], t[N];
int mu[N], b[N], p[N], pn;
struct que{
	int n, m, a, id;
}q[TN];

void init(){
	mu[1] = 1;
	for (int i = 2; i < N; ++ i){
		if (!b[i])
			p[++pn] = i, mu[i] = mod-1;
		for (int j = 1; j <= pn; ++ j){
			if (1LL*i*p[j] >= N) break;
			b[i*p[j]] = 1;
			if (i%p[j] == 0){mu[i*p[j]] = 0; break;}
			else mu[i*p[j]] = mod-mu[i];
		}
	}
}

bool cmp(que x, que y){
	return x.a < y.a;
}

void add(int &x, int y){x += y; if (x >= mod) x -= mod;}

int lowbit(int x){return x&(-x);}

void upd(int x, int y){
	for (int i = x; i < N; i += lowbit(i))
		add(t[i], y);
}

int qry(int x){
	int ret = 0;
	for (int i = x; i; i -= lowbit(i))
		add(ret, t[i]);
	return ret;
}

int sum(int x){return 1LL*(x+1)*x/2%mod;}

int solve(int n, int m, int a){
	if (n > m) swap(n, m);
	int ret = 0;
	for (int i = 1, j; i <= n; i = j+1){
		j = min(n, min(n / (n / i), m / (m / i)));
		add(ret, 1LL * sum(n/i) * sum(m/i) % mod * (qry(j)-qry(i-1)+mod) % mod);
	}
	return ret;
}

int main()
{
	init();
	scanf("%d", &T);
	for (int i = 1; i <= T; ++ i)
		scanf("%d%d%d", &q[i].n, &q[i].m, &q[i].a), q[i].id = i;
	sort(q + 1, q + T + 1, cmp);
	for (int i = 1, k = 1; i <= T; ++ i){
		for (; k <= q[i].a; ++ k)
			for (int j = k; j < N; j += k)
				upd(j, 1LL * k * mu[j/k] % mod * (j/k) % mod * (j/k) % mod);
		ans[q[i].id] = solve(q[i].n, q[i].m, q[i].a);
	}
	for (int i = 1; i <= T; ++ i)
		printf("%d\n", ans[i]);
	return 0;
}