[BZOJ2820]YY的GCD

2820: YY的GCD

Time Limit: 10 Sec  Memory Limit: 512 MB Submit: 2248  Solved: 1207 [Submit][Status][Discuss]

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

 

由对称性，不妨设$n\le m$

那么$ans=\sum_{\left(p\right)IsPrime}\sum_{i=1}^n\sum_{j=1}^m\left[gcd\left(i,j\right)=p\right]$

由其它一些经典的题目（比如[POI2007]ZAP）很轻易地可以得到

$ans=\sum_{\left(p\right)IsPrime}\sum_{d=1}^{\lfloor\frac{n}{p}\rfloor}\mu\left(d\right)\lfloor\frac{n}{pd}\rfloor\lfloor\frac{m}{pd}\rfloor$

由于在前$n$个数中质数的个数约为$\frac{n}{lnn}$个

那么即使是加上了分块优化时间复杂度也是$O\left(\frac{n\sqrt{n}}{lnn}\right)$肯定会TLE

考虑设$pd=T$

那么$ans=\sum_{T=1}^{n}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{p\mid T}\mu\left(\frac{T}{p}\right)$

那么加上分块优化和前缀和处理后单次询问的时间复杂度就是$O\left(\sqrt{n}\right)$

总时间复杂度$O\left(T\sqrt{n}\right)$

#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long ll;
const int maxn = 10000000 + 10;
bool mark[maxn] = {0};
int pri[664580], prn = 0;
int mu[maxn], sum[maxn];
void shai(){
    mu[1] = 1;
    for(int i = 2; i < maxn; i++){
        if(!mark[i]){
            mu[i] = -1;
            pri[++prn] = i;
        }
        for(int j = 1; j <= prn && i * pri[j] < maxn; j++){
            mark[i * pri[j]] = true;
            if(i % pri[j] == 0){
                mu[i * pri[j]] = 0;
                break;
            }
            else mu[i * pri[j]] = -mu[i]; // = mu[pri[j]] * mu[i]
        }
    }
    for(int i = 1; i <= prn; i++)
        for(int j = 1; j * pri[i] < maxn; j++) sum[pri[i] * j] += mu[j];
    sum[0] = 0;
    for(int i = 1; i < maxn; i++)
        sum[i] += sum[i - 1];
}
int main(){
    shai();
    int T, n, m;
    ll ans;
    scanf("%d", &T);
    while(T--){
        scanf("%d %d", &n, &m);    
        if(n > m) swap(n, m);
        ans = 0;
        for(int p, i = 1; i <= n; i = p + 1){
            p = min(n / (n / i), m / (m / i));
            ans += (ll) (sum[p] - sum[i - 1]) * (n / p) * (m / p);
        }
        printf("%lld\n", ans);
    }
    return 0;
}