[BZOJ2820]YY的GCD

题目大意：
　　对于给定的$n,m(n,m\leq10^7)$，求$为质数\displaystyle\sum_{x=1}^n\sum_{y=1}^m[\gcd(x,y)为质数]$。

思路：
　　设$p=\gcd(x,y),x=ap,y=bp$，则：
$$
\begin{align*}
原式&=\sum_{1\leq p\leq\min(n,m)且p为质数}\sum_{a=1}^{\lfloor\frac{n}{p}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{p}\rfloor}[\gcd(a,b)=1]\\
&=\sum_{1\leq p\leq\min(n,m)且p为质数}\sum_{a=1}^{\lfloor\frac{n}{p}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{p}\rfloor}\sum_{d|\gcd(a,b)}\mu(d)\\
&=\sum_{1\leq p\leq\min(n,m)且p为质数}\sum_{p=1}^{\lfloor\frac{\min(n,m)}{p}\rfloor}\mu(p)\sum_{a=1}^{\lfloor\frac{n}{dp}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{dp}\rfloor}1\\
&=\sum_{1\leq p\leq\min(n,m)且p为质数}\sum_{p=1}^{\lfloor\frac{\min(n,m)}{p}\rfloor}\mu(p)\lfloor\frac{n}{dp}\rfloor\lfloor\frac{m}{dp}\rfloor
\end{align*}
$$
　　令$k=dp$，则：
$$
\begin{align*}
原式&=\sum_{k=1}^{\min(n,m)}\sum_{p为质数且p|k}\mu(p)\lfloor\frac{n}{k}\rfloor\lfloor\frac{m}{k}\rfloor
\end{align*}
$$
　　其中$为质数且\displaystyle\sum_{p为质数且p|k}\mu(p)$可以预处理，$k$可以用数论分块枚举。

#include<cstdio>
#include<cctype>
#include<algorithm>
typedef long long int64;
inline int getint() {
    register char ch;
    while(!isdigit(ch=getchar()));
    register int x=ch^'0';
    while(isdigit(ch=getchar())) x=(((x<<2)+x)<<1)+(ch^'0');
    return x;
}
const int N=10000001,M=664580;
bool vis[N];
int mu[N],p[M];
int64 sum[N];
inline void sieve() {
    mu[1]=1;
    for(register int i=2;i<N;i++) {
        if(!vis[i]) {
            p[++p[0]]=i;
            mu[i]=-1;
        }
        for(register int j=1;j<=p[0]&&i*p[j]<N;j++) {
            vis[i*p[j]]=true;
            if(i%p[j]==0) {
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
    }
    for(register int i=1;i<=p[0];i++) {
        for(register int j=1;j*p[i]<N;j++) {
            sum[j*p[i]]+=mu[j];
        }
    }
    for(register int i=1;i<N;i++) {
        sum[i]+=sum[i-1];
    }
}
int main() {
    sieve();
    for(register int i=getint();i;i--) {
        const int n=getint(),m=getint(),lim=std::min(n,m);
        int64 ans=0;
        for(register int i=1,j;i<=lim;i=j+1) {
            j=std::min(n/(n/i),m/(m/i));
            ans+=(sum[j]-sum[i-1])*(n/i)*(m/i);
        }
        printf("%lld\n",ans);
    }
    return 0;
}