[BZOJ]1101: [POI2007]Zap

Description
　　FGD正在破解一段密码，他需要回答很多类似的问题：对于给定的整数a,b和d，有多少正整数对x,y，满足x<=a，y<=b，并且gcd(x,y)=d。作为FGD的同学，FGD希望得到你的帮助。
Input
　 第一行包含一个正整数n，表示一共有n组询问。（1<=n<= 50000）接下来n行，每行表示一个询问，每行三个正整数，分别为a,b,d。（1<=d<=a,b<=50000）
Output
　　对于每组询问，输出到输出文件zap.out一个正整数，表示满足条件的整数对数。
Sample Input
2
4 5 2
6 4 3
Sample Output
3
2
//对于第一组询问，满足条件的整数对有(2,2)，（2,4），（4，2）。对于第二组询问，满足条件的整数对有（6,3），（3,3）。

参照别人提供的莫比乌斯函数的一个性质

$$\sum_{d|n}μ(d)=[n==1]$$
在这一题中，我们需要利用这一个性质。
如果不经变换，这题要求的东西就是
$$\sum_{i=1}^a\sum_{j=1}^b[gcd(i,j)==1]$$
我们先令
$$a'= \lfloor a/d\rfloor,b'=\lfloor b/d\rfloor$$
然后将原式转为
$$\sum_{i=1}^a'\sum_{j=1}^b'[gcd(i,j)==1]$$
再用前面提到的性质将式子再变换一下
$$\sum_{i=1}^a'\sum_{j=1}^b'\sum_{d|gcd(i,j)}\mu(d)$$
因为
$$d|gcd(i,j)\Leftrightarrow(d|i~\And~d|j)$$
然后结果就可以用一个简单的式子表示出来了
$$\sum_{i=1}^{min(a',b')}\mu(i)\lfloor a'/i\rfloor \lfloor b'/i\rfloor$$
其中会有一大段μ(i)要加的数量一致，用前缀和直接一起搞定
上述一切思想都是别人的智慧结晶，我只是一个搬运工,打完可以去水2301,双倍经验

#include <cstdio>
using namespace std;
inline char tc(void)
{
    static char fl[10000],*A=fl,*B=fl;
    return A==B&&(B=(A=fl)+fread(fl,1,10000,stdin),A==B)?EOF:*A++;
}
inline int read(void)
{
    int a=0;static char c;
    while((c=tc())<'0'||c>'9');
    while(c>='0'&&c<='9')
        a=a*10+c-'0',c=tc();
    return a;
}
inline void swap(int &a,int &b)
{
    int t=a;a=b,b=t;return ;
}
inline int min(int a,int b)
{
    return a<b?a:b;
}
inline int max(int a,int b)
{
    return a>b?a:b;
}
int n,cnt,mu[50001],sum[50001],p[50001],ans,pos;
char b[50001];
inline void prepare(void)
{
    mu[1]=1;
    for (int i=2;i<=50000;++i)
    {
        if(!b[i])
            p[++cnt]=i,mu[i]=-1;
        for (int j=1;j<=cnt&&p[j]*i<=50000;++j)
        {
            b[i*p[j]]=1;
            if(i%p[j]==0) 
            {
                mu[i*p[j]]=0;
                break;
            }
            else
                mu[i*p[j]]=-mu[i];
        }
    }
    for(int i=1;i<=50000;i++)
        sum[i]=sum[i-1]+mu[i];
    return ;
}
int cal(int n,int m)
{
    if(n>m)swap(n,m);
    ans=0;
    for(int i=1;i<=n;i=pos+1)
    {
        pos=min(n/(n/i),m/(m/i));
        ans+=(sum[pos]-sum[i-1])*(n/i)*(m/i);
    }
    return ans;
}
int main(void)
{
    register int i,a,b,c;
    prepare(),n=read();
    while(n--)
        a=read(),b=read(),c=read(),printf("%d\n",cal(a/c,b/c));
    return 0;
}