P2257 YY的GCD（莫比乌斯反演，多组询问）

题目传送门

tips：本题是题P2568的强化版，如果不会本题可以先试试P2568。

题意：

给出$T(1\le T\le 10^4)$ 组询问，每组询问一对$n,m$，让你求出：
a n s = ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) ∈ { p r i m e } ] ans=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\left[gcd(i,j)∈\{prime\}\right] ans=i=1∑n​j=1∑m​[gcd(i,j)∈{prime}]

思路：

假定你已经会了弱化版，并且不知道如何处理多组数据。
我们从这个式子开始： ∑ p ∈ { p r i m e } ∑ d = 1 ⌊ n p ⌋ μ ( d ) ⌊ n p d ⌋ ∗ ⌊ m p d ⌋ \sum\limits_{p∈\{prime\}}\sum\limits_{d=1}^{\lfloor \frac{n}{p} \rfloor}\mu(d)\lfloor\cfrac{n}{pd}\rfloor*\lfloor\cfrac{m}{pd}\rfloor p∈{prime}∑​d=1∑⌊pn​⌋​μ(d)⌊pdn​⌋∗⌊pdm​⌋。
当我们推到这个式子的时候，已经可以通过弱化版的题目了，直接枚举质数然后整除分块求答案就行了。
但是本题有$10^4$组数据，我们不可能再这样枚举，因为算一下时间复杂度明显太大了，于是我们想办法优化。我们设$Q=pd$，则
原 式 = ∑ p ∈ { p r i m e } ∑ d = 1 ⌊ n p ⌋ μ ( d ) ⌊ n Q ⌋ ∗ ⌊ m Q ⌋ 原式=\sum\limits_{p∈\{prime\}}\sum\limits_{d=1}^{\lfloor \frac{n}{p} \rfloor}\mu(d)\lfloor\cfrac{n}{Q}\rfloor*\lfloor\cfrac{m}{Q}\rfloor 原式=p∈{prime}∑​d=1∑⌊pn​⌋​μ(d)⌊Qn​⌋∗⌊Qm​⌋
然后我们枚举$Q$，可以得到式子 ∑ Q = 1 n ⌊ n Q ⌋ ∗ ⌊ m Q ⌋ ∗ ∑ p ∈ { p r i m e } , p ∣ Q μ ( Q p ) \sum\limits_{Q=1}^{n}\lfloor\cfrac{n}{Q}\rfloor*\lfloor\cfrac{m}{Q}\rfloor*\sum\limits_{p∈\{prime\},p|Q}\mu\left(\cfrac{Q}{p}\right) Q=1∑n​⌊Qn​⌋∗⌊Qm​⌋∗p∈{prime},p∣Q∑​μ(pQ​)。
我们设 f ( Q ) = ∑ p ∈ { p r i m e } , p ∣ Q μ ( Q p ) f(Q)=\sum\limits_{p∈\{prime\},p|Q}\mu\left(\cfrac{Q}{p}\right) f(Q)=p∈{prime},p∣Q∑​μ(pQ​)。
然后我们可以利用欧拉筛$O(n)$先把$\mu (i)$求出来，然后利用埃氏筛的思想，在$O(nloglogn)$的复杂度内把$f(Q)$求出来，并做一下前缀和，为整除分块做准备。

前面那一坨当然可以用整除分块求了，预处理的复杂度为$O(nloglogn)$，计算的复杂度为$O(T\ast \sqrt{n})$。

$Code$

// Author : ACfunhsl
// Time : 2021/5/17 14:13:11
#define int long long
const int N = 1e7+50;
const int inf = 0x3f3f3f3f;
const int mod = 1e9+7;
bool ok[N];
int p[N],cnt=0,mu[N],f[N];
void euler()
{
	mu[1] = 1;
	for(int i=2; i<N; i++)
	{
		if(!ok[i])
		{
			p[++cnt] = i;
			mu[i] = -1;
		}
		for(int j=1; j<=cnt&&i*p[j]<N; j++)
		{
			ok[i*p[j]] = 1;
			if(i%p[j]==0) break;
			mu[i*p[j]] = -mu[i];
		}
	}
	for(int i=1;i<=cnt;i++)
		for(int j=1;j*p[i]<N;j++)
			f[j*p[i]] += mu[j];
	for(int i=1; i<N; i++)
		f[i] += f[i-1];
}
int cal(int n,int m)
{
	int res = 0;
	for(int l=1,r; l<=n; l=r+1)
	{
		r = min(n/(n/l),m/(m/l));
		res += (n/l)*(m/l)*(f[r] - f[l-1]);
	}
	return res;
}
signed main()
{
	euler();
	int t;
	cin>>t;
	while(t--)
	{
		int n,m;
		cin>>n>>m;
		if(n>m)
			swap(n,m);
		cout<<cal(n,m)<<endl;
	}

	return 0;
}