YY的GCD [Bzoj 2820]

题目地址请点击——

---

YY的GCD

---

Description

神犇 YY 虐完数论后给傻× kAc 出了一题。
给定 $N$, $M$,求 $1\leq x\leq N$ , $1\leq y\leq 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 \leq 10000$
$N, M \leq 10000000$

---

Solution

令 $M<N$，则

ans=∑p∑p|iμ(ip)(⌊Mi⌋⋅⌊Ni⌋)=∑p∑i=1⌊Mi⌋μ(i)⋅(⌊Mip⌋⋅⌊Nip⌋)=∑i=1M⌊Mi⌋⋅⌊Ni⌋⋅∑p|iμ(ip)

$$\begin{align*} ans&=\sum_{p}\sum_{p|i}\mu(\frac{i}{p})(\lfloor\frac{M}{i}\rfloor·\lfloor\frac{N}{i}\rfloor)\\ &=\sum_{p}\sum_{i=1}^{\lfloor\frac{M}{i}\rfloor}\mu(i)·(\lfloor\frac{M}{ip}\rfloor·\lfloor\frac{N}{ip}\rfloor)\\ &=\sum_{i=1}^M\lfloor\frac{M}{i}\rfloor·\lfloor\frac{N}{i}\rfloor·\sum_{p|i}\mu(\frac{i}{p}) \end{align*}$$
设 $f_i=\sum_{p|i}\mu(\frac{i}{p})$，
暴力枚举每个质数，更新其倍数的 $f$ 值，由 $\sum^{n}_{i=1} \frac{n}{i} ≈ n·lnn$ 这个结论易知每个质数更新时是均摊 $O(lnn)$ 的，而质数个数接近 $\frac{n}{lnn}$，故暴力枚举的时间复杂度为 $O(n)$。

令 $g_i=\sum_{j=1}^if_j$，即 $f_i$ 的前缀和，则可以进行分块优化。

---

Code

#include <iostream>
#include <cstdio>

#define LL long long
#define Min(x,y) ((x)<(y)?(x):(y))

using namespace std;

LL T;
LL n,m;
bool no_prime[(10000000)+10];
LL prime[(10000000)+10];
short mu[(10000000)+10];
LL sum[(10000000)+10];
LL f[(10000000)+10];

void start(){
    mu[1]=1;
    for(LL i=2;i<=(10000000);i++){
        if(!no_prime[i]){
            prime[++prime[0]]=i;
            mu[i]=-1;
        }
        for(LL j=1;prime[j]*i<=(10000000);j++){
            no_prime[prime[j]*i]=true;
            if(i%prime[j]==0){mu[prime[j]*i]=0;break;}
            mu[prime[j]*i]=-mu[i];
        }
    }
    for(LL i=1;i<=prime[0];i++)
        for(LL j=1;prime[i]*j<=(10000000);j++)
            f[prime[i]*j]+=mu[j];
    for(LL i=1;i<=(10000000);i++)sum[i]=sum[i-1]+f[i];
}

int main(){
    start();
    scanf("%lld",&T);
    for(LL i=1;i<=T;i++){
        scanf("%lld%lld",&m,&n);
        LL ans=0;
        if(m>n)swap(m,n);
        for(LL j=1,it;j<=m;j=it+1){
            it=Min(m/(m/j),n/(n/j));
            ans+=(m/j)*(n/j)*(sum[it]-sum[j-1]);
        }
        printf("%lld\n",ans);
    }
}