【洛谷2257/BZOJ2820】YY的GCD（数论/莫比乌斯函数）

题目：

洛谷2257

预备知识：莫比乌斯定理（懵逼乌斯定理）

(mu*1=epsilon)（证bu明hui略zheng）
其中（我校学长把(epsilon(x))叫单位函数但是为什么我没百度到qwq）

[epsilon(x)=egin{cases}1 & x=1\ 0 & x eq1\ end{cases}]

[mu(x)=egin{cases}1 & x=1\ 0 & 存在质数p使p^2|x\ (-1)^k & k是x质因数的个数 end{cases}]

那个(*)是迪利克雷卷积，换成人话就是

[epsilon(n)=sum_{d|n}mu(d) ]

我觉得用这种方式理解莫比乌斯定理比设两个函数容易

分析：

这题莫比乌斯定理的经典用例。
本文中默认(N>M)
默认(p)是质数

显然如果(gcd(i,j)=p)，那么(gcd(frac{i}{p}, frac{j}{p})=1)
那么题目所求可以转换成下面的式子

[sum_{p}^Nsum_i^{N/p}sum_j^{M/p}epsilon(gcd(i,j)) ]

其中（我校学长把这个叫单位函数但是我没百度到qwq）

[epsilon(x)=egin{cases}1 & x=1\ 0 & x eq1\ end{cases}]

根据莫比乌斯反演定理，上面的式子就可以变成

[sum_{p=2}^Nsum_i^{N/p}sum_j^{M/p}sum_{d|gcd(i,j)}mu(d) ]

改变一下枚举顺序，用(d·i)表示原来的(i)，(d·j)表示原来的(j)，得到

[sum_{p=2}^Nsum_d^{N/p}sum_i^{N/pd}sum_j^{M/pd}mu(d) ]

可以发现(mu(d))和(i)、(j)没半毛钱关系，仅仅是乘上(i)和(j)可以取的值的数量
也就是

[sum_{p=2}^Nsum_d^{N/p}mu(d)*lfloorfrac{N}{pd} floor*lfloorfrac{M}{pd} floor ]

令(T=pd)，枚举T，上式可变成

[sum_T^Nsum_{p|T}mu(frac{T}{p})*lfloorfrac{N}{T} floor*lfloorfrac{M}{T} floor ]

设$$g(x)=sum_{p|x}mu(frac{x}{p})$$
则上式就是

[sum_T^Nlfloorfrac{N}{T} floor*lfloorfrac{M}{T} floor*g(T) ]

现在考虑如何求(g(x))这个函数。
首先，对于任意质数(p)，显然(g(p)=mu(1)=1)
然后，对于任意合数(n=kp_0)（(p_0)是质数）(g(n))中显然存在(mu(frac{n}{p_0}))也就是(mu(k))这一项
当({p_0}|k)，也就是(p_0^2|n)，对于任意(p|k)且(p eq p_0)，(mu(frac{n}{p}))中一定有(p_0^2)这个质数平方因子。根据(mu(x))的定义，(mu(frac{n}{p})=0)
所以此时(g(n)=mu(k))
当(p_0)不能整除(k)，对于任意(p|k)，(mu(frac{n}{p}))比(mu(frac{k}{p}))多了(p_0)这个质因子。根据(mu(x))的定义(mu(frac{n}{p})=-mu(frac{k}{p}))
所以此时(g(n)=-g(k)+mu(k))
总结一下

[g(x)=egin{cases}1 & x是质数\ mu(k) & x=kp且p能整除k\ -g(k)+mu(k) & x=kp且p不能整除k end{cases}]

显然这个函数可以用线性筛求

[sum_T^Nlfloorfrac{N}{T} floor*lfloorfrac{M}{T} floor*g(T) ]

再来看这个式子，既然(g(T))可以直接预处理并(O(1))查询，那么计算这个式子的时间复杂度就是枚举(T)的复杂度(O(N))
我会做啦！

别急，这题还有(T)组询问，所以复杂度是O(不可过)(O(NT))，这个过不了。
但是我们可以发现(lfloorfrac{N}{T} floor*lfloorfrac{M}{T} floor)在(T)的一段区间内是不变的，所以可以给(g(T))算个前缀和然后分段计算，据说复杂度是(O(sqrt N T))的（我不会证），这样就可以过了

代码：

#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std;

namespace zyt
{
	typedef long long ll;
	const int N = 1e7 + 10, M = 7e5;
	bool mark[N];
	int cnt, prime[M], phi[N], mu[N];
   	ll g[N];
	void init()
	{
		mu[1] = 1;
		for (int i = 2; i < N; i++)
		{
			if (!mark[i])
				prime[cnt++] = i, mu[i] = -1, g[i] = 1;
			for (int j = 0; j < cnt && (ll)i * prime[j] < N; j++)
			{
				int k = i * prime[j];
				mark[k] = true;
				if (i % prime[j] == 0)
				{
					mu[k] = 0;
					g[k] = mu[i];
					break;
				}
				else
				{
					mu[k] = -mu[i];
					g[k] = -g[i] + mu[i];
				}
			}
		}
		for (int i = 2; i < N; i++)
			g[i] += g[i - 1];
	}
	void work()
	{
		int T;
		init();
		scanf("%d", &T);
		while (T--)
		{
			int n, m, pos = cnt;
			ll ans = 0;
			scanf("%d%d", &n, &m);
			if (n > m)
				swap(n, m);
			for (int t = 1; t <= n;)
			{
				int tmp = min(n / (n / t), m / (m / t));
				ans += (g[tmp] - g[t - 1]) * (n / t) * (m / t);
				t = tmp + 1;
			}
			printf("%lld
", ans);
		}
	}
}
int main()
{
	zyt::work();
	return 0;
}