P2522 [HAOI2011]Problem b 莫比乌斯反演+除法优化

题目描述

　　对于给出的n个询问，每次求有多少个数对(x,y)，满足a≤x≤b，c≤y≤d，且gcd(x,y) = k，gcd(x,y)函数为x和y的最大公约数。
　　输入格式：　第一行一个整数n，接下来n行每行五个整数，分别表示a、b、c、d、k
　　输出格式：　共n行，每行一个整数表示满足要求的数对(x,y)的个数
　　输入样例#1：
　　2
　　2 5 1 5 1
　　1 5 1 5 2
　　输出样例#1：
　　14
　　3
说明：100%的数据满足：1≤n≤50000，1≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000
洛谷传送门：https://daniu.luogu.org/problem/show?pid=2522

解题报告：

　　　　莫比乌斯反演之经典题，同时此题也体现了除法优化的优美。
　　　　题目要求：
　

$$\sum_{i=a}^{b}\sum_{j=c}^{d}[gcd(i,j)==k]$$
　　　我们先不去管a和c，先来研究这个（下面均设n小于m）：
　 $$\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==k]$$
　 　我们设 $f(k)$表示 $1\le i \le n$&& $1\le j \le m$中满足 $gcd(i,j)==k$的数对 $(i,j)$的个数。显然 $f(k)$很难直接求。
　 　接着我们设 $g(k)$表示 $1\le i \le n$&& $1\le j \le m$中满足 $k|gcd(i,j)$的数对 $(i,j)$的个数。我们很显然能够得出下面这个式子：
　 　 $$g(k)=\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}f(i*k)= \lfloor\frac{n}{k}\rfloor * \lfloor\frac{m}{k}\rfloor$$
　 　那么我们将 $f(k)和g(k)$反演可得：
　 　 $$f(k)=\sum_{d=1}^{\lfloor\frac{n}{k}} \mu(d) g(k*d)=\sum_{d=1}^{\lfloor\frac{n}{k}} \mu(d)*\lfloor\frac{n}{k*d}\rfloor * \lfloor\frac{m}{k*d}\rfloor$$
　 　优美的式子已经出来了，但如果直接这样算的话单次询问是 $O(\lfloor\frac{n}{k}\rfloor)$的效率，只能拿70分。那么我们下面就要用到 除法优化了。
　 　
　 　我们发现不管是 $\lfloor\frac{n}{k*d}\rfloor$还是 $\lfloor\frac{m}{k*d}\rfloor$，都有一大段相同的，这就可以让我们不去重复计算。除法优化的具体写法如下（你先需要求个 $\mu(x)$的前缀和）：
　
接着我们再回过头来看原式子，我们发现可以用容斥的方法得出答案，即：
答案=ans(b,d)-ans(a-1,d)-ans(c-1,b)+(a-1,c-1)
那么我们就解决了。

代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
int mu[50005],prime[50005];
int size,n,a,b,c,d,k;
bool notprime[50005];
void init(){
    mu[1]=1;notprime[1]=true;
    for(int i=2;i<=50000;i++){
        if(!notprime[i]) 
            prime[++size]=i,mu[i]=-1;
        for(int j=1,k;j<=size;j++){
            k=prime[j]*i; 
            if(k>50000)break;
            notprime[k]=true;
            if(i%prime[j]==0){
                mu[k]=0;
                break;
            } 
            else 
                mu[k]=-mu[i];
        }
    }
    for(int i=2;i<=50000;i++)mu[i]+=mu[i-1];
}
long long work(int x,int y){  //除法优化
    x/=k,y/=k;
    long long ans=0;
    int i=1,j,lim=min(x,y);
    while(i<=lim){
        j=min(x/(x/i),y/(y/i));
        ans+=(long long)(mu[j]-mu[i-1])*(long long)(x/i)*(y/i);
        i=j+1;
    }
    return ans;
}
int main(){
    scanf("%d",&n);
    init();
    for(int i=1;i<=n;i++){
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        printf("%lld\n",work(b,d)-work(a-1,d)-work(c-1,b)+work(a-1,c-1));
    }
    return 0;
}