BZOJ2301 Problem B 【莫比乌斯反演】【分块】

题目链接：https://www.lydsy.com/JudgeOnline/problem.php?id=2301

题解：
对于 $a \le x \le b，c \le y \le d$，这个条件，我们发现比较难以处理，这时候我们可以利用二维前缀和的思想，记 $x \le b，y \le d$ 时的答案为 $A[b][d]$，那么$Ans=A[b][d]-A[b][c-1]-A[a-1][d]+A[a-1][c-1]$，这样就把问题简化成了两个条件的问题。
下面推一下式子：
$Ans=\Sigma_{i=1}^{b}\ \Sigma_{j=1}^{d}\ (gcd(i,j)==k)$

$=\Sigma_{i=1}^{b}\ \Sigma_{j=1}^{d}\ (gcd(i/k,j/k)==1)$

我们发现 $gcd(X,Y)==1$，可以用卷积的单位元来做。

$=\Sigma_{i=1}^{b}\ \Sigma_{j=1}^{d}\ (e(gcd(i/k,j/k))==1)$

$=\Sigma_{i=1}^{b}\ \Sigma_{j=1}^{d}\ \Sigma_{T|(i/k),T|(j/k)}\mu(T)$

$=\Sigma_{T|(i/k),T|(j/k)}\ \mu(T)\ \Sigma_{i=1}^{b/(k*T)}\ \Sigma_{j=1}^{d/(k*T)}*1$

$=\Sigma_{T|(i/k),T|(j,k)}\ \mu(T)*(b/(k*T))*(d/(k*T))$

对于$(b/(k*T))*(d/(k*T))$，对于$(b/k)，(d/k)$，将其看做$n、m$，发现可以分块即可。

原因：$(b/(k*T))=((b/k)/T)$ –> 将 $b/k$ 看做 $n$，一项即为 $n/T$，
可以从小到大枚举 $T$，相同的块可以直接跳过。
现在我们的任务变成了：如何求一个块的最大元素。这个严格证明参考了某大神的blog。
设这个块的 $T$ 为 $i$
设 $n=p*i+q$，那么对于 $i$，若 $i++$，则 $q$ 会相应减少 $p$，而 $q$ 最多减少 $q/p$个$p$，
所以此时的 $i'=i+q/p=(pi+q)/p$，即为 $n/(n/i)$，每次分块做即可
PS：小心爆 int，如果不分块的话极限数据是50000^2级别的，会TLE。
关于求$\mu()$函数，由于其为积性函数，所以可以用线性筛解决，详见代码。

代码：

// by Balloons
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define mpr make_pair
#define debug() puts("okkkkkkkk")
#define rep(i,a,b) for(int (i)=(a);(i)<=(b);(i)++)

using namespace std;

typedef long long LL;

const int inf = 1 << 30;

const int maxn=50000;
int notpm[maxn+5],pm[maxn+5],mu[maxn+5],pcnt=0;
int tes;
int a,b,c,d,k;

/*
推导之后，我们发现就是一个式子：
sigma T=1 to min(x/k,y/k) miu(T) * (b/(k*T)) * (d/(k*T))
对于后面两项，我们发现可以分块。
原因：(b/(k*T))=((b/k)/T) --> 将 b/k 看做 n，一项即为 n/T，
可以从小到大枚举 T，相同的块可以直接跳过。
现在我们的任务变成了：如何求一个块的最大元素。设这个块的 T为i 
设 n=pi+q，那么对于 i，若 i++，则 q 会相应减少 p，而 q 最多减少 q/p个p，
所以此时的 i'=i+q/p=(pi+q)/p，即为 n/(n/i) 
*/

void xxs(){
    notpm[1]=1;
    mu[1]=1;
    for(int i=2;i<=maxn;i++){
        if(!notpm[i]){
            pm[++pcnt]=i;mu[i]=-1;
        }
        for(int j=1;j<=pcnt&&pm[j]*i<=maxn;j++){
            notpm[pm[j]*i]=1;
            if(i%pm[j]==0){
                mu[i*pm[j]]=0;
                break;
            }
            mu[i*pm[j]]=-mu[i];
        } 
        mu[i]+=mu[i-1]; 
    }
} 

LL solve(int x,int y){
    x/=k;y/=k;  // 推导公式需要 
    LL tmp=0;
    for(int i=1,j;i<=min(x,y);i=j+1){
        j=min(x/(x/i),y/(y/i)); // 详见上 
        tmp+=1ll*(mu[j]-mu[i-1])*(x/i)*(y/i);
    }
    return tmp;
}

int main(){
//  freopen("2301.in","r",stdin);
//  freopen("2301.out","w",stdout);
    scanf("%d",&tes);
    xxs(); 
    while(tes--){
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        printf("%lld\n",solve(b,d)-solve(b,c-1)-solve(a-1,d)+solve(a-1,c-1));
    }

    return 0;
}