【bzoj1101】Zap【神奇的∑】

传送门：
http://www.lydsy.com/JudgeOnline/problem.php?id=1101
求$\sum_{i=1}^a\sum_{j=1}^b[gcd(i,j)==d]$.
设$a \le b$.
令$a'=\lfloor \frac ad \rfloor,b'=\lfloor \frac bd \rfloor$.
除以$d$得：
$\quad\sum_{i=1}^{a'}\sum_{i=1}^{b'}[gcd(i,j)==1]$
$=\sum_{i=1}^{a'}\sum_{i=1}^{b'}\sum_{d|gcd(i,j)}\mu(d)$
$=\sum_{i=1}^{a'}\sum_{i=1}^{b'}\sum_{d|i \;\text{and}\;d|j}\mu(d)$
$=\sum_{d=1}^{a'}\mu(d)\sum_{i=1}^{\lfloor a'/d\rfloor}\sum_{j=1}^{\lfloor b'/d\rfloor} 1$
$=\sum_{d=1}^{a'}\mu(d)\left\lfloor\frac{a'}d\right\rfloor\left\lfloor\frac{b'}d\right\rfloor$
容易知道$\lfloor a'/d\rfloor$单调不上升，且最多有种$2\sqrt {a'}$不同的取值。所以按取值分成$O(\sqrt n)$个段分别处理，一个连续段内的和可以用预处理出的莫比乌斯函数前缀和求出.

if(a>b)swap(a,b);
int nex;
ll ans=0;
for (int i=1;i<=a;i=nex+1){
    nex=min(a/(a/i),b/(b/i));
    ans+=1ll*(s[nex]-s[i-1])*(a/i)*(b/i);
}

本文一些地方借鉴了策爷的文章（orz jcvb）

//FGD正在破解一段密码，他需要回答很多类似的问题：
//对于给定的整数a,b和d，有多少正整数对x,y，
//满足x<=a，y<=b，并且gcd(x,y)=d。作为FGD的同学，
//FGD希望得到你的帮助。
#include<cstdio>
#include<cstring>
#include<iostream>
using namespace std;
const int maxn=50001;
int tot,prime[maxn],mu[maxn];
bool isprime[maxn];
void genPrime(int n){
    memset(isprime,true,sizeof isprime);
    mu[1]=1;
    for(int i=2;i<=n;++i){
        if(isprime[i]){
            prime[tot++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<tot&&i*prime[j]<=n;++j){
            isprime[i*prime[j]]=false;
            if(i%prime[j]==0){
                //phi[i*prime[j]]=phi[i]*prime[j];
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=-mu[i];//phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
}
int gcd(int a,int b){return b?gcd(b,a%b):a;}
int main(){
    int a,b,d,n;
    genPrime(50000);
    for(int i=1;i<=50000;++i) mu[i]+=mu[i-1];
    scanf("%d",&n);
    while(n--){
        scanf("%d%d%d",&a,&b,&d);
        a/=d;b/=d;
        if(a>b) swap(a,b);
        int ans=0;

        for(int i=1,nex;i<=a;i=nex+1){
            nex=min(a/(a/i),b/(b/i));
            ans+=(mu[nex]-mu[i-1])*(a/i)*(b/i);
        }
        printf("%d\n",ans);
    }
//  ans=0;
//  for(int i=1;i<=a;++i) for(int j=1;j<=b;++j) ans+=(gcd(i,j)==1);
//  cout<<ans<<"\n";
}