[BZOJ2301][HAOI2011]Problem b（莫比乌斯反演）

=== ===

这里放传送门

=== ===

题解

感觉自己在懵（fei）懂（chang）无（sha）知（bi）的时候用某种诡异的瞎推式子方法做了一些实际上是反演的题目？！！比如这个题。。。

这道题要求满足$a \leq x \leq b，c \leq y \leq d$的$gcd(x,y)=k$的数对$(x,y)$的数目，那么很容易想到可以转化成满足${a \over k} \leq x \leq {b \over k}，{c \over k} \leq y \leq {d \over k}$的$gcd(x,y)=1$的数对$(x,y)$的数目。再把每个询问用二维前缀和的方式拆成四个，就变成了要求解$1 \leq x \leq n，1 \leq y \leq m$范围内的互质数对个数。接下来我们要开始画柿子啦~

要求的东西是：

$$\sum_{i=1}^n\sum_{j=1}^m[(x,y)==1]$$
利用公式 $[n=1]=\sum\limits_{d|n}\mu(d)$（诶话说把这公式套上来是不是就是反演的一种），我们可以把式子化成： $$\sum_{i=1}^n\sum_{j=1}^m\sum_{d|(i,j)}\mu(d)$$
可以发现 $d|(i,j)$的充要条件实际上就是d是i和j的公约数，那么可以继续把式子化成： $$\sum_{i=1}^n\sum_{j=1}^m\sum_{d=1}^{min(n,m)}[d|i][d|j]\mu(d)$$
把那些 $\sum$移动一下排列成一个好看的样子： $$\sum_{d=1}^{min(n,m)}\sum_{i=1}^n[d|i]\sum_{j=1}^m[d|j]\mu(d)$$
可以发现 $\sum\limits_{i=1}^n[d|i]$就是在求 $1..n$的范围内有多少d的倍数，这个东西显然就是 $\lfloor{n \over d}\rfloor$。那么式子可以变成一个愉悦的形态： $$\sum_{d=1}^{min(n,m)}\lfloor{n \over d}\rfloor\lfloor{m \over d}\rfloor\mu(d)$$
用分块枚举除法取值的方法，我们就可以在 $O(\sqrt n + \sqrt m)$的时间复杂内解决单组询问辣！

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 60000
using namespace std;
int T,a,b,c,d,k,mu[60010],prm[60010];
bool ext[60010];
long long solve(long long n,long long m){
    long long ans,tail;
    if (n>m) swap(n,m);
    ans=0;
    for (int i=1;i<=n;i=tail+1){
        tail=min(n/(n/i),m/(m/i));
        ans+=(n/i)*(m/i)*(mu[tail]-mu[i-1]);
    }
    return ans;
}
int main()
{
    mu[1]=1;
    for (int i=2;i<=N;i++){
        if (ext[i]==false){
            prm[++prm[0]]=i;
            mu[i]=-1;
        }
        for (int j=1;j<=prm[0];j++){
            if (i*prm[j]>N) break;
            ext[i*prm[j]]=true;
            if (i%prm[j]==0){
                mu[i*prm[j]]=0;break;
            }else mu[i*prm[j]]=-mu[i];
        }
        mu[i]+=mu[i-1];
    }
    scanf("%d",&T);
    for (int wer=1;wer<=T;wer++){
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        a--;c--;b=b/k;d=d/k;a=a/k;c=c/k;
        if (b==0||d==0) {printf("0\n");continue;}
        printf("%I64d\n",solve(b,d)-solve(a,d)-solve(c,b)+solve(a,c));
    }
    return 0;
}