容斥

题意不在说明。

解法

求$x<=a,y<=b$且$gcd(x,y)=d$，那么对式子除以个$d$，得到 $x’<=\frac{a}{d},y’<=\frac{b}{d},gcd(x’,y’)=1$
所以就可以利用容斥求。
那么方案数是多少呢？我们令$a/d=a’,b/d=b’$，那么答案应该就是$a’ \cdot b’-\frac{a’}{2} \cdot \frac{b’}{2}-\frac{a’}{3} \cdot \frac{b’}{3}-……+\frac{a’}{6} \cdot \frac{b’}{6}+……-……$。
就是减去有一个公因子的个数，加上有两个的，也就是容斥原理的应用。
但是可以发现，如果我们这样做的话，每次询问的复杂度是$O(N)$，所以会爆掉。
思考如何去优化，根据答案的形式，可以发现这其实是一个莫比乌斯函数。
所以答案可以转化成$a’ \cdot b’-\sum\limits_{i=1}^{min(a’,b’)}\frac{a’}{i} \cdot \frac{b’}{i} \cdot mobius(i)$。可见复杂度是$O(N)$。
考虑优化。我们将$\frac{a’}{i}$拿出来看，就是$\frac{a’}{1} \cdot \frac{a’}{2} \cdot …… \cdot \frac{a’}{n}$。虽然有n项相乘，但其中最多只有$2\sqrt{a’}$项，所以就可以将原来分成$\sqrt{n}$段，利用前缀和，理论上可以将求和一部分简化到$O(\sqrt{N})$的复杂度。

证明

为什么要分成$\sqrt{n}$段呢？很好证明，将其分成两段$\frac{a’}{1} \cdot \frac{a’}{2} \cdot \frac{a’}{\sqrt{a’}}$与$\frac{a’}{\sqrt{a’}+1} \cdot …… \cdot \frac{a’}{n}$，前面只有只有$sqrt{a’}$个不同答案，后面只有$sqrt{a’}$中答案，所以总共就有$2\sqrt{a’}$

如何去求

这里有一个经典的函数，我们设$g(x)$，他的值就是能够使$\frac{a}{x}=\frac{a}{g(x)}$成立的最大的数，也就是$a$除以一个数，值与$\frac{a}{x}$相等，且最大的一个数，也就有$\frac{a}{x}>\frac{a}{g(x)+1}$。
这个$g(x)=\frac{a}{\frac{a}{x}}$。

请注意，本文所有除法都带有向下取整。

证明$\frac{a}{x}=\frac{a}{g(x)}$。

因为有$g(x)=\frac{a}{\frac{a}{x}}$，将分子中的向下取整去掉，那么分子变大整个式子变小，化简就得$g(x)=x$。那么$g(x)>=x$。所以$\frac{a}{g(x)}<=\frac{a}{x}$。那就只需证明$\frac{a}{\frac{a}{\frac{a}{x}}}>=\frac{a}{x}$。将第二个取整符号去掉，分子变大，原式变小，化简成$\frac{a}{x}=\frac{a}{x}$，所以$\frac{a}{g(x)}>=\frac{a}{x}$。
即证。

证明$\frac{a}{x}>\frac{a}{g(x)+1}$

这个证明很有必要，不证的话就只代表一个，就不能保证复杂度为$O(\sqrt{N})$。
先令$a=kx+r,0<=r< x$
那么原式就是$k>\frac{a}{\frac{a}{k}+1}<=>k \cdot (\frac{a}{k}+1)>a$。
再令$a=pk+q,0<=q< k$。
带回刚刚的式子中，就是$k \cdot (p+1)>pk+q->kp+k>pk+q->k>q$。
即证
那么所以，对于答案的求和，只需要跳$\sqrt{N}$次。

参考代码

#include<bits/stdc++.h>
using namespace std;
const int N=50010;
typedef long long ll;
int T;
int n;
int cnt;
int a,b,d;
int sum[N];
bool vis[N];
int pri[N],mobius[N];
void Mobius() {
    mobius[1]=1;
    for(int i=2; i<N; ++i) {
        if(!vis[i]) {
            pri[cnt++]=i;
            mobius[i]=-1;
        }
        for(int j=0; pri[j]*i<N; ++j) {
            int t=pri[j]*i;
            vis[t]=1;
            if(i%pri[j]==0) {
                mobius[t]=0;
                break;
            }
            mobius[t]=mobius[i]*-1;
        }
    }
    for(int i=1; i<N; ++i) sum[i]=sum[i-1]+mobius[i];
    return ;
}
int main() {
    Mobius();
    scanf("%d",&T);
    while(T--) {
        scanf("%d%d%d",&a,&b,&d);
        a/=d,b/=d;
        n=min(a,b);
        ll res=0;
        for(int l=1,r; l<=n; l=r+1) {
            r=min(n,min(a/(a/l),b/(b/l)));
            res+=(sum[r]-sum[l-1])*1ll*(a/l)*(b/l);
        }
        printf("%lld\n",res);
    }
    return 0;
}