luogu2257 YY的GCD（莫比乌斯反演）

luogu2257

题目描述：神犇YY虐完数论后给傻×kAc出了一题
     给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对
     kAc这种傻×必然不会了，于是向你来请教……
     多组输入

输入格式：第一行一个整数T 表述数据组数
     接下来T行，每行两个正整数，表示N, M

输出格式：T行，每行一个整数表示第i组数据的结果

输入样例：
2
10 10
100 100

输出样例：
30
2791

解析：莫比乌斯反演中比较基础的题，也是我第一次做的莫比乌斯反演题。
   设f(n)表示gcd为n的数对的对数，F(n)表示gcd为n或n的倍数的数对的对数。
   那么显然 (F(n) = sum_{n|d}f(d))
   反演之后得 (f(n) = sum_{n|d}μ(⌊frac{d}{n}⌋)F(d))
   而 F(n) = (⌊frac{N}{n}⌋⌊frac{M}{n}⌋)
   所以可得 (f(n) = sum_{n|d}μ(⌊frac{d}{n}⌋)⌊frac{N}{n}⌋⌊frac{M}{n}⌋)
   所求的 ans = (sum_{n∈prim}f(n))
   ans = (sum_{n∈prim}sum_{n|d}μ(⌊frac{d}{n}⌋)⌊frac{N}{n}⌋⌊frac{M}{n}⌋)
   换一个枚举项，枚举 (⌊frac{d}{n}⌋) 那么 ans = (sum_{n∈prim}sum_{d=1}^{min(⌊frac{N}{n}⌋,⌊frac{M}{n}⌋)}μ(d)⌊frac{N}{dn}⌋⌊frac{M}{dn}⌋)
   将dn换成t，则 ans = (sum_{T=1}^{min(N,M)}sum_{t|T,t∈prim}μ(⌊frac{T}{t}⌋)⌊frac{N}{T}⌋⌊frac{M}{T}⌋)
   即 ans = (sum_{T=1}^{min(N,M)}⌊frac{N}{T}⌋⌊frac{M}{T}⌋sum_{t|T,t∈prim}μ(⌊frac{T}{t}⌋))
   到了这里，只要套上一个整除分块就做好了
   μ可以直接筛出来，维护一个前缀和即可

代码如下：

#include<cstdio>
#include<algorithm>
#define ll long long
using namespace std;

const int maxn = 1e7 + 5;
int q, mu[maxn], mark[maxn], prim[maxn], tot, g[maxn], n, m;
ll sum[maxn], ans;

int read(void) {
	char c; while (c = getchar(), c < '0' || c >'9'); int x = c - '0';
	while (c = getchar(), c >= '0' && c <= '9') x = x * 10 + c - '0'; return x;
}

void get_mu(void) { //筛mu函数 
	mu[1] = 1; mark[1] = 1;
	  for (int i = 2; i <= maxn - 5; ++ i) {
	  	if (!mark[i]) mu[i] = -1, prim[++ tot] = i;
	  	for (int j = 1; j <= tot; ++ j) {
	  		if (prim[j] * i > maxn - 5) break;
	  		mark[prim[j] * i] = 1;
	  		if (i % prim[j] == 0) break;
	  		  else mu[i * prim[j]] = -mu[i];
		  }
	  }
	  for (int i = 1; i <= tot; ++ i) //统计mu函数的前缀和 
	    for (int j = 1; j * prim[i] <= maxn - 5; ++ j) g[j * prim[i]] += mu[j];
	  for (int i = 1; i <= maxn - 5; ++ i) sum[i] = sum[i - 1] + g[i];
}

int main() {
	q = read();
	get_mu();
	  while (q --) {
	  	ans = 0; 
	  	n = read(); m = read();
	  	  for (int l = 1, r; l <= min(n, m); l = r + 1) { //整数分块统计答案 
	  	  	  r = min(n / (n / l), m / (m / l));
	  	  	  ans += 1ll * (n / l) * (m / l) * (sum[r] - sum[l - 1]);
			}
		printf("%lld
", ans);
	  }
	return 0;
}