SP4491 PGCD - Primes in GCD Table

题目大意

求 $\sum\limits_{p \in primes}\sum\limits_{x=1}^{n}\sum\limits_{y=1}^{m}[\gcd(x,y)=p]$。

$T$ 组数据。

其中 $1 \leq T \leq 10,1 \leq n,m \leq 10^7$。

解题思路

求：

$$\sum\limits_{p \in primes}\sum\limits_{x=1}^{n}\sum\limits_{y=1}^{m}[\gcd(x,y)=p]$$

化简得：

$$\sum\limits_{p \in primes}\sum\limits_{x=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\sum\limits_{y=1}^{\left\lfloor\frac{m}{p}\right\rfloor}[\gcd(x,y)=1]$$

根据：

$$\sum\limits_{d|n} \mu(d) = \begin{cases} 1 && n=1 \\ 0 && n \ne 1\end{cases}$$

有：

$$\sum\limits_{p \in primes}\sum\limits_{x=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\sum\limits_{y=1}^{\left\lfloor\frac{m}{p}\right\rfloor}\sum\limits_{k|i,k|j}\mu(k)$$

先枚举 $k$，有：

$$\sum\limits_{p \in primes}\sum\limits_{k=1}^{\min(\left\lfloor\frac{n}{p}\right\rfloor,\left\lfloor\frac{m}{p}\right\rfloor)}\sum\limits_{x=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\sum\limits_{y=1}^{\left\lfloor\frac{m}{p}\right\rfloor}\mu(k)[k|i][k|j]$$

那么就有：

$$\sum\limits_{p \in primes}\sum\limits_{k=1}^{\min(\left\lfloor\frac{n}{p}\right\rfloor,\left\lfloor\frac{m}{p}\right\rfloor)}\left\lfloor\frac{n}{pk}\right\rfloor\left\lfloor\frac{m}{pk}\right\rfloor\mu(k)$$

设 $T=pk$，有：

$$\sum\limits_{T=1}^{\min(n,m)}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum\limits_{p\in primes,p|T} \mu(\frac{T}{p})$$

求前缀和再数论分块即可。

CODE

#include <bits/stdc++.h>

using namespace std;

#define int long long

const int _ = 1e7;

int T, n, m;

bool vis[_ + 7];

int tot, pri[_ + 7], g[_ + 7], sum[_ + 7], mu[_ + 7];

void init()
{
	mu[1] = 1;
	for(int i = 2; i < _; ++i)
	{
		if(!vis[i])
		{
			pri[++tot] = i;
			mu[i] = -1;
			g[i] = 1;
		}
		for(int j = 1; j <= tot, i * pri[j] < _; ++j)
		{
			vis[i * pri[j]] = 1;
			if(i % pri[j])
			{
				mu[i * pri[j]] = -mu[i];
				g[i * pri[j]] = mu[i] - g[i];
			}
			else
			{
				mu[i * pri[j]] = 0;
				g[i * pri[j]] = mu[i];
				break;
			}
		}
		sum[i] = sum[i - 1] + g[i];
	}
}

signed main()
{
	init();
	scanf("%lld", &T);
	while(T--)
	{
		scanf("%lld%lld", &n, &m);
		if(n > m) swap(n, m);
		int ans = 0ll;
		for(int i = 1, last; i <= n; i = last + 1)
		{
			last = min(n / (n / i), m / (m / i));
			ans += (n / i) * (m / i) * (sum[last] - sum[i - 1]);
		}
		printf("%lld\n", ans);
	}
	return 0;
}