YY的GCD

$枚举gcd(x,y)=p，研究有多少对(x,y)满足：$

$p|x且p|y且1\le x\le n且1\le y\le m且gcd(x/p,y/p)=1$

$设x=p\ast a,y=p\ast b$

$约束变为1\le a\le \lfloor n/p\rfloor 且1\le b\le \lfloor m/p\rfloor 且gcd(a,b)=1$

${\sum }_{i=1}^{\lfloor n/p\rfloor }{\sum }_{j=1}^{\lfloor m/p\rfloor }[gcd(i,j)=1]={\sum }_{i=1}^{\lfloor n/p\rfloor }{\sum }_{j=1}^{\lfloor m/p\rfloor }{\sum }_{d|gcd(i,j)}\mu (d)$

$[gcd(i,j)=1]={\sum }_{d|gcd(i,j)}\mu (d)$

$原式=\sum (\mu (d)\ast (\lfloor n/(p\ast d)\rfloor )\ast (\lfloor m/(p\ast d)\rfloor ))$

$令t=p\ast d则d=t/p$

$原式={\sum }_{p|t}(\mu (t/p)\ast (\lfloor n/t\rfloor )\ast (\lfloor m/t\rfloor ))$

$定义函数h(x)={\sum }_{p|x}\mu (x/p)$

${\sum }_{p是质数}原式=\sum (h(t)\ast (\lfloor n/t\rfloor )\ast (\lfloor m/t\rfloor ))$

$只要我们能求h(x)的区间和，就能用整除分块$

#include <cstdio>
#include <iostream>
using namespace std;
#define LL long long
#define ULL unsigned long long

template <typename T> void read (T &x) { x = 0; T f = 1;char tem = getchar ();while (tem < '0' || tem > '9') {if (tem == '-') f = -1;tem = getchar ();}while (tem >= '0' && tem <= '9') {x = (x << 1) + (x << 3) + tem - '0';tem = getchar ();}x *= f; return; }
template <typename T> void write (T x) { if (x < 0) {x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0'); }
template <typename T> void print (T x, char ch) { write (x); putchar (ch); }
template <typename T> T Max (T x, T y) { return x > y ? x : y; }
template <typename T> T Min (T x, T y) { return x < y ? x : y; }
template <typename T> T Abs (T x) { return x > 0 ? x : -x; }

const int Maxn = 1e7;
const int Maxm = 1e6;

int cnt, primes[Maxm + 5], mu[Maxn + 5], minp[Maxn + 5], h[Maxn + 5];
bool vis[Maxn + 5];
LL pre[Maxn + 5];
void Euler () {
	mu[1] = 1;
	for (int i = 2; i <= Maxn; i++) {
		if (vis[i] == 0) {
			primes[++cnt] = i;
			mu[i] = -1;
			minp[i] = i;
		}
		for (int j = 1; j <= cnt; j++) {
			if (primes[j] > Maxn / i) break;
			vis[i * primes[j]] = 1;
			minp[i * primes[j]] = primes[j];
			if (i % primes[j] == 0) {
				mu[i * primes[j]] = 0;
				break;
			}
			mu[i * primes[j]] = mu[i] * mu[primes[j]];
		}
	}
	
	for (int i = 1; i <= Maxn; i++) {
		int p = i;
		while (p > 1) {
			h[i] += mu[i / minp[p]];
			int tmp = minp[p];
			while (p % tmp == 0) 
				p /= tmp;
		}
		pre[i] = pre[i - 1] + h[i];
	}
}

int t, n, m;
LL Calc (int l, int r) {
	return pre[r] - pre[l - 1];
}

int main () {
	Euler ();
	
	read (t);
	while (t--) {
		read (n); read (m);
		int l = 1, r;
		LL res = 0;
		while (l <= Min (n, m)) {
			r = Min (n / (n / l), m / (m / l));
			res = (res + (LL)(n / l) * (m / l) * Calc (l, r));
			l = r + 1;
		}
		print (res, '\n');
	}
	return 0;
}