BZOJ2301: [HAOI2011]Problem b[莫比乌斯反演容斥原理]【学习笔记】

2301: [HAOI2011]Problem b

Time Limit: 50 Sec Memory Limit: 256 MB
Submit: 4032 Solved: 1817
[Submit][Status][Discuss]

Description

对于给出的n个询问，每次求有多少个数对(x,y)，满足a≤x≤b，c≤y≤d，且gcd(x,y) = k，gcd(x,y)函数为x和y的最大公约数。

Input

第一行一个整数n，接下来n行每行五个整数，分别表示a、b、c、d、k

Output

共n行，每行一个整数表示满足要求的数对(x,y)的个数

Sample Input

2

2 5 1 5 1

1 5 1 5 2

Sample Output

14

3

HINT

100%的数据满足：1≤n≤50000，1≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000

研究了好长时间差不多明白了，第一道莫比乌斯反演，好多值得学习的东西

首先，由容斥原理易得答案为

cal(b,d,k)-cal(a-1,d,k)-cal(b,c-1,k)+cal(a-1,c-1,k)

这个问题等价于询问有多少个数对(x,y)满足1<=x<=[n/k],1<=y<=[m/k]且x与y互质
考虑莫比乌斯反演，
f(i)为1<=x<=n,1<=y<=m且gcd(x,y)=i的数对(x,y)的个数
F(i)为1<=x<=n,1<=y<=m且i|gcd(x,y)的数对(x,y)的个数
显然，F(i)=(n/i)*(m/i) 整除，并且F(i)=Σ{i|d} f(d) 是倍数和
反演后，f(i)=Σ{i|d} miu(d/i)*F(d)=Σ{i|d} miu(d/i)*(n/d)*(m/d)
但这样每个询问复杂度是O(n)
观察式子，发现[n/d] 最多有2sqrt(n) 个取值(整除....一段相同参考链接)
那么 (n/d)*(m/d)就至多有2sqrt(n)+2sqrt(m)个取值 (当然不是乘起来，因为对于一个n只有一个值而不是2sqrt(n)个)
计算每个询问时枚举这2sqrt(n)+2sqrt(m)个取值，因为一个取值是一段，要乘一段miu的和，所以对莫比乌斯函数维护一个前缀和，就可以在sqrt(n)时间内出解

【WT1(WT是从小新那里学来的....发现竟然是问题的首字母)：】

f(k)=Σ{k|d} miu(d/k)*(n/d)*(m/d)这个式子怎么计算？

d是k的倍数，取值k,2*k,3*k,...,t*k

f(k)=Σ{i=1..n/k} miu(i)*(n/(k*i))*(m/(k*i)) //注意，【整除满足 x/a/b=a/(a*b)】

更一般的：

f(k)=Σ{k|d} miu(d/k)*F(d)

--> f(k)=Σ{i=1..n/k}miu(i)*F(i*k)

【WT2 】如何按照整除取值相同分段？

当前除法为n/i，与它相同的上界到n/(n/i)

为什么？我想了好久，最后的方法是

考虑n是一段区间，n=p*i+q，被分成p段长为i的

i每增加1 q就减少p，（这时候整除的取值没有改变），最多能减少q/p个，那么此时i=i+q/p=(i*p+q)/p=n/(n/i)

注意：miu的区间和*(n/i)*(m/i)可能会溢出，对拍都没有发现.......

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N=5e4+5;
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1; c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0'; c=getchar();}
    return x*f;
}
int n,a,b,c,d,k;
bool notp[N];
ll p[N],mu[N];
void sieve(){
    mu[1]=1;
    for(int i=2;i<=N-1;i++){
        if(!notp[i]) p[++p[0]]=i,mu[i]=-1;
        for(int j=1;j<=p[0]&&i*p[j]<=N-1;j++){
            int t=i*p[j];
            notp[t]=1;
            if(i%p[j]==0){
                mu[t]=0;
                break;
            }
            mu[t]=-mu[i];
        }
    }
    for(int i=1;i<=N-1;i++) mu[i]+=mu[i-1];
}
ll cal(int n,int m,int k){
    n/=k;m/=k;
    if(n>m) swap(n,m);
    ll ans=0;int r=0;
    for(int i=1;i<=n;i=r+1){
        r=min(n/(n/i),m/(m/i));
        ans+=(mu[r]-mu[i-1])*(n/i)*(m/i);
    }
    return ans;
}
int main(int argc, const char * argv[]) {
    //freopen("in.txt","r",stdin);
    //freopen("1.out","w",stdout);
    sieve();
    int T=read();
    while(T--){
        a=read();b=read();c=read();d=read();k=read();
        printf("%lld\n",cal(b,d,k)-cal(a-1,d,k)-cal(b,c-1,k)+cal(a-1,c-1,k));
    }
    
    return 0;
}

附：还有另一种思考的角度，从莫比乌斯函数的角度考虑，殊途同归

复制鏼爷的题解

推导：

令 $a'=\left \lfloor a/d \right \rfloor,b'=\left \lfloor b/d \right \rfloor$
$\sum_{i=1}^a\sum_{j=1}^b[\text{gcd}(i,j)==d]$
$=\sum_{i=1}^{a'}\sum_{j=1}^{b'}[\text{gcd}(i,j)==1]$
用莫比乌斯函数的性质把求和的式子换掉，
$=\sum_{i=1}^{a'}\sum_{j=1}^{b'}\sum_{d|\text{gcd}(i,j)}\mu(d)$
其中 $d|\text{gcd}(i,j) \Leftrightarrow d|i \,\text{and} \,d|j$ ，更换求和指标，
$=\sum_{d=1}^{a'}\mu(d)\left \lfloor a'/d \right \rfloor \left \lfloor b'/d \right \rfloor$
容易知道 $\left\lfloor a'/d\right \rfloor$ 单调不上升，且最多有 $2\sqrt{a'}$ 种不同的取值。所以按取值分成 $O(\sqrt{n})$ 个段分别处理，一个连续段内的和可以用预处理出的莫比乌斯函数前缀和求出