[BZOJ2820]YY的GCD

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

Solution

假定n<m
$$\sum_{isprime(p)}\sum_{a=1}^n\sum_{b=1}^mgcd(a,b)==p$$
$$\sum_{isprime(p)}\sum_{a=1}^{\left \lfloor \frac{n}{p} \right \rfloor}\sum_{b=1}^{\left \lfloor \frac{m}{p} \right \rfloor}gcd(a,b)==1$$

$$\sum_{isprime(p)}\sum_{a=1}^{\left \lfloor \frac{n}{p} \right \rfloor}\sum_{b=1}^{\left \lfloor \frac{m}{p} \right \rfloor}\sum_{d|gcd(a,b)}\mu(d)$$

$$\sum_{isprime(p)}\sum_{a=1}^{\left \lfloor \frac{n}{p} \right \rfloor}\sum_{b=1}^{\left \lfloor \frac{m}{p} \right \rfloor}\sum_{d|a \land d|b}\mu(d)$$
$$\sum_{isprime(p)}\sum_{d=1}^{\left \lfloor \frac{n}{p} \right \rfloor}\mu(d){\left \lfloor \frac{n}{pd} \right \rfloor}{\left \lfloor \frac{m}{pd} \right \rfloor}$$
将pd设为k
$$\sum_{k=1}^{n}\sum_{isprime(p) \land p|k}\mu(\frac{k}{p}){\left \lfloor \frac{n}{k} \right \rfloor}{\left \lfloor \frac{m}{k} \right \rfloor}$$
$$\sum_{k=1}^{n}F(k){\left \lfloor \frac{n}{k} \right \rfloor}{\left \lfloor \frac{m}{k} \right \rfloor}$$
可以线筛预处理F,或者暴力枚举质数（这个复杂度我不能确定）,按照素数粗略个数n/logn以及调和级数求和nlogn来看暴力的复杂度接近On
处理完F以后就是喜闻乐见的下底函数分块

Code

#include <bits/stdc++.h>
using namespace std;
#define rep(i, l, r) for (int i = (l); i <= (r); i++)
#define per(i, r, l) for (int i = (r); i >= (l); i--)
#define X first
#define Y second
#define MS(_) memset(_, 0, sizeof(_))
#define PB push_back
#define MP make_pair
#define debug(...) fprintf(stderr, __VA_ARGS__)
template<typename T> inline void read(T &x){
    x = 0; T f = 1; char ch = getchar();
    while (!isdigit(ch)) {if (ch == '-') f = -1; ch = getchar();}
    while (isdigit(ch))  {x = x * 10 + ch - '0'; ch = getchar();}
    x *= f;
}
typedef long long ll;
const int N = 11111111;

bool check[N];
int prime[N];
ll mu[N], f[N], sum[N];
inline void prework(){ int tot = 0;
    MS(check); mu[1] = f[1] = 1ll;
    rep(i, 2, N){
        if (!check[i]){prime[++tot] = i; mu[i] = -1; f[i] = 1ll;}
        rep(j, 1, tot){
            if (1ll * i * prime[j] > 1ll * N) break;
            check[i * prime[j]] = true;
            if (i % prime[j] == 0){
                mu[i * prime[j]] = 0; f[i * prime[j]] = mu[i]; break;
            }else{
                mu[i * prime[j]] = -mu[i]; f[i * prime[j]] = mu[i] - f[i];
            }
        }
    }
    rep(i, 1, N) sum[i] = sum[i - 1] + f[i] * 1ll;
}

inline ll cal(ll n, ll m){
    ll ans = 0; int pos = 0;
    for (int i = 2; i <= n; i = pos + 1){
        pos = min(n/(n/i), m/(m/i));
        ans += (sum[pos] - sum[i-1]) * (n/i) * (m/i) * 1ll;
    }
    return ans;
}
int main(){
    prework();
    int case_cnt; read(case_cnt);
    while (case_cnt--){ ll n, m;
        read(n); read(m); if (n > m) swap(n, m);
        printf("%lld\n", cal(n, m));
    }
    return 0;
}