[HAOI2011]Problem b&&[POI2007]Zap

题目大意：
　　$q(qleq50000)$组询问，对于给定的$a,b,c,d(a,b,c,dleq50000)$，求$displaystylesum_{i=a}^bsum_{j=c}^d[gcd(i,j)=k]$。

思路：
　　首先可以利用容斥原理，将$(a,b,c,d)$的询问拆成$(b,d)$、$(a-1,d)$、$(b,c-1)$和$(a-1,c-1)$四个询问，对于询问$(n,m)$，有：
$$
egin{align*}
&sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=k]\
&=sum_{k=1}^{min(n,m)}sum_{i=1}^{lfloorfrac{n}{k} floor}sum_{j=1}^{lfloorfrac{m}{k} floor}[gcd(i,j)=1]\
&=sum_{k=1}^{min(n,m)}sum_{i=1}^{lfloorfrac{n}{k} floor}sum_{j=1}^{lfloorfrac{m}{k} floor}sum_{d|gcd(i,j)}mu(d)\
&=sum_{k=1}^{min(n,m)}sum_{d=1}^{min(lfloorfrac{n}{k} floor,lfloorfrac{m}{k} floor)}mu(d)sum_{1leq ileqlfloorfrac{n}{k} floor且d|i}sum_{1leq ileqlfloorfrac{m}{k} floor且d|j}1\
&=sum_{k=1}^{min(n,m)}sum_{d=1}^{min(lfloorfrac{n}{k} floor,lfloorfrac{m}{k} floor)}mu(d)lfloorfrac{lfloorfrac{n}{k} floor}{d} floorlfloorfrac{lfloorfrac{m}{k} floor}{d} floor
end{align*}
$$
　　$mu$的前缀和可以预处理。因此到最后一步时，复杂度已经是$O(n)$的了，然而$qleq50000$，还是会TLE。考虑某些时候对于某一范围内的$k$，$lfloorfrac{lfloorfrac{n}{k} floor}{d} floor$和$lfloorfrac{lfloorfrac{m}{k} floor}{d} floor$分别各自相等，因此我们可以将相等的情况一起考虑，这样总共有$O(sqrt n)$种情况，时间复杂度$O(qsqrt n)$。

#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=50001,M=5134;
bool vis[N];
int mu[N],sum[N],p[M];
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<N;i++) {
        sum[i]=sum[i-1]+mu[i];
    }
}
inline int64 calc(int n,int m,const int &k) {
    int64 ret=0;
    n/=k,m/=k;
    const int lim=std::min(n,m);
    for(register int i=1,j;i<=lim;i=j+1) {
        j=std::min(n/(n/i),m/(m/i));
        ret+=(int64)(sum[j]-sum[i-1])*(n/i)*(m/i);
    }
    return ret;
}
int main() {
    sieve();
    for(register int i=getint();i;i--) {
        const int a=getint(),b=getint(),c=getint(),d=getint(),k=getint();
        printf("%lld
",calc(b,d,k)-calc(a-1,d,k)-calc(b,c-1,k)+calc(a-1,c-1,k));
    }
    return 0;
}