bzoj 2820 YY的GCD - 莫比乌斯反演 - 线性筛

Description

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

Input

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

Output

T行，每行一个整数表示第i组数据的结果

Sample Input

2

10 10

100 100

Sample Output

30

2791

HINT

T = 10000

N, M <= 10000000

　　题目大意 (题目太简洁不需要大意)。

　　根据常用套路，我们有

　　


　　显然TLE，为了更好的分段计算，所以考虑把后面两个向下取整取整的分数挪出来。于是得到了下面这个优美的式子:



　　现在就来考虑后面那个求和。

　　考虑莫比乌斯函数的性质和p为质数。

　　1)当T存在一个质因子的指数超过2，显然对答案贡献为0

　　2)当T存在超过1个质因子的指数为2，显然对答案贡献为0

　　3)当T只存在1个质因子的指数为2，那么当且仅当p的值为这个质因子时，对答案才有贡献，这个贡献为-1的不同质因子次幂。

　　4)当T的质因子的指数均为1，那么对答案的贡献为T的质因子个数乘-1的不同质因子个数减一 次幂。

　　显然可以用线性筛预处理出这个的值的前缀和，这样就可以做到单次查询根号级别。

　　至于如何用线性筛预处理，我记录了一个enable表示是否存在1个质因子的指数大于等于2和defp表示不同的质因子的个数。详细的过程请看代码(假设读者的数学水平可以轻松把这个秒掉)。

Code

/**
* bzoj
* Problem#2820
* Accepted
* Time:4708ms
* Memory:141912k
*/
#include <bits/stdc++.h>
#ifndef WIN32
#define Auto "%lld"
#else
#define Auto "%I64d"
#endif
using namespace std;
typedef bool boolean;

#define LL long long

const int limit = 1e7;

int num = 0;
int prime[1000000];
boolean vis[limit + 1];
boolean enable[limit + 1];
LL g[limit + 1];
int defp[limit + 1];
inline void Euler() {
memset(vis, false, sizeof(vis));
memset(enable, true, sizeof(enable));
g[0] = g[1] = 0;
for(int i = 2; i <= limit; i++) {
if(!vis[i])    prime[num++] = i, g[i] = 1, defp[i] = 1;
for(int j = 0; j < num && i * 1LL * prime[j] <= limit; j++) {
int c = i * prime[j];
vis[c] = true;
if(i % prime[j]) {
defp[c] = defp[i] + 1;
enable[c] = enable[i];
g[c] = (enable[i]) ? ((defp[c]) * ((defp[c] & 1) ? (1) : (-1))) : (-g[i]);
} else {
defp[c] = defp[i];
enable[c] = false;
g[c] = (enable[i]) ? ((defp[c] & 1) ? (-1) : (1)) : (0);
break;
}
}
}
for(int i = 2; i <= limit; i++)
g[i] += g[i - 1];
//    cout << num << endl;
//    for(int i = 1; i <= 20; i++)
//        cout << i << " " << g[i] << " " << defp[i] << endl;
}

int n, m;
inline void init() {
scanf("%d%d", &n, &m);
}

inline void solve() {
if(n > m)    swap(n, m);
LL res = 0;
for(int i = 1, j; i <= n; i = j + 1) {
j = min(n / (n / i), m / (m / i));
res += (n / j) * 1LL * (m / j) * (g[j] - g[i - 1]);
}
printf(Auto"\n", res);
}

int T;
int main() {
Euler();
scanf("%d", &T);
while(T--) {
init();
solve();
}
return 0;
}