[HAOI2011]Problem b&&[POI2007]Zap

题目大意：
　　$q(q\leq50000)$组询问，对于给定的$a,b,c,d(a,b,c,d\leq50000)$，求$\displaystyle\sum_{i=a}^b\sum_{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)$，有：
$$
\begin{align*}
&\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=k]\\
&=\sum_{k=1}^{\min(n,m)}\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[\gcd(i,j)=1]\\
&=\sum_{k=1}^{\min(n,m)}\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum_{d|\gcd(i,j)}\mu(d)\\
&=\sum_{k=1}^{\min(n,m)}\sum_{d=1}^{\min(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor)}\mu(d)\sum_{1\leq i\leq\lfloor\frac{n}{k}\rfloor且d|i}\sum_{1\leq i\leq\lfloor\frac{m}{k}\rfloor且d|j}1\\
&=\sum_{k=1}^{\min(n,m)}\sum_{d=1}^{\min(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor)}\mu(d)\lfloor\frac{\lfloor\frac{n}{k}\rfloor}{d}\rfloor\lfloor\frac{\lfloor\frac{m}{k}\rfloor}{d}\rfloor
\end{align*}
$$
　　$\mu$的前缀和可以预处理。因此到最后一步时，复杂度已经是$O(n)$的了，然而$q\leq50000$，还是会TLE。考虑某些时候对于某一范围内的$k$，$\lfloor\frac{\lfloor\frac{n}{k}\rfloor}{d}\rfloor$和$\lfloor\frac{\lfloor\frac{m}{k}\rfloor}{d}\rfloor$分别各自相等，因此我们可以将相等的情况一起考虑，这样总共有$O(\sqrt n)$种情况，时间复杂度$O(q\sqrt 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\n",calc(b,d,k)-calc(a-1,d,k)-calc(b,c-1,k)+calc(a-1,c-1,k));
    }
    return 0;
}