[POI2007]ZAP-Queries （莫比乌斯反演+整除分块）

[POI2007]ZAP-Queries

---

---

$ solution: $

唉，数论实在有点烂了，昨天还会的，今天就不会了，周末刚证明的，今天全忘了，还不如早点写好题解。

这题首先我们可以列出来答案就是：

$ ans=sum_{i=1}^{a}{sum_{j=1}^{b}{[gcd(i,j)==d]}} $

我们发现后面那个东西（只有 $ gcd(i,j)d $ 时才为一）跟莫比乌斯很像，莫比乌斯是（只有 $ n $ 1 才为一），所以我们再尝试转化一下（把d转化成1）：

$ ans=sum_{i=1}^{frac{a}{d}}{sum_{j=1}^{frac{b}{d}}{[gcd(i,j)==1]}} $

于是我们就可以把后面那个东西用莫比乌斯函数的第一条性质转换成这样：

$ ans=sum_{i=1}^{frac{a}{d}}{sum_{j=1}^{frac{b}{d}}{sum_{k|gcd(i,j)}{mu(k)}}} $

但是这样显然还不够，我们想办法把莫比乌斯的式子挪到前面去：

$ ans=sum_{k}^{min(a,b)}{mu(k)}{sum_{i=1}^{frac{a}{d}}{sum_{j=1}^{frac{b}{d}}{[k|gcd(i,j)]}}} $

这个其实就相当于我们从小到大枚举k，但是我们是从上面那个式子转化过来的，所以必须满足 $ [k|gcd(i,j)] $ 这个条件。好了，现在我们肉眼观察一下，发现如下的东西：

$ {sum_{i=1}^{frac{a}{d}}{sum_{j=1}^{frac{b}{d}}{[k|gcd(i,j)]}}}=lfloor frac{lfloor frac{a}{d} floor}{k} floor imes lfloor frac{lfloor frac{b}{d} floor}{k} floor=lfloor frac{a}{d imes k} floor imes lfloor frac{b}{d imes k} floor $

$ ans=sum_{k}^{min(a,b)}{mu(k) imes lfloor frac{a}{d imes k} floor imes lfloor frac{b}{d imes k} floor} $

然后我们发现这样子的复杂度是 $ O(min(a,b)) $ 的，然而它的询问次数太多。于是出现了一个很奇妙的东西：整除分块（又叫数论分块）。举个栗子：

$ frac{10}{1}=10 $

$ frac{10}{2}=5 $

$ frac{10}{3}=3 $

$ frac{10}{4}=frac{10}{5}=2 $

$ frac{10}{6}=frac{10}{7}=frac{10}{8}=frac{10}{9}=frac{10}{10}=1 $

我们发现分子相同分母越大，则出现相同结果的概率越高，所以我们可以一次求出某一段相同结果的左端点和右端点（假设这一段的结果都为y，则这一段的最右端就是用分子除以y得到的值），从而使算法效率变高，这就是整除分块。

---

---

$ code: $

#include<iostream>
#include<cstdio>
#include<iomanip>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>

#define ll long long
#define db double
#define inf 0x7fffffff
#define rg register int

using namespace std;

int Q;
int pr[50005];
int mu[50005];
bool use[50005];

inline int min(const rg &x,const rg &y){
    if(x<y)return x; else return y;
}

inline int qr(){
    register char ch; register bool sign=0; rg res=0;
    while(!isdigit(ch=getchar())) if(ch=='-')sign=1;
    while(isdigit(ch)) res=res*10+(ch^48),ch=getchar();
    return sign?-res:res;
}

inline void get_mu(int x){
    rg t=0; mu[1]=1;
    for(rg i=2;i<=x;++i){
        if(!use[i])mu[i]=-1,pr[++t]=i;
        for(rg j=1;j<=t;++j){
            if(i*pr[j]>x)break;
            use[i*pr[j]]=1;
            if(!(i%pr[j]))break;
            else mu[i*pr[j]]=-mu[i];
        }
    }
    for(rg i=2;i<=x;++i) mu[i]+=mu[i-1];
}

int main(){
    //freopen(".in","r",stdin);
    //freopen(".out","w",stdout);
    Q=qr();
    get_mu(50000);
    while(Q--){
        rg a=qr(),b=qr(),k=qr();
        a/=k; b/=k;
        rg r,n=min(a,b),ans=0;
        for(rg l=1;l<=n;l=r+1){
            r=min(a/(a/l),b/(b/l));
            ans+=((a/l)*(b/l)*(mu[r]-mu[l-1]));
        }printf("%d
",ans);
    }
    return 0;
}