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

2021-11-07 09:00:44

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 5 1 5 1

1 5 1 5 2

Sample Output

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+;

inline int read(){

    char c=getchar();int x=,f=;

    while(c<''||c>''){if(c=='-')f=-; c=getchar();}

    while(c>=''&&c<=''){x=x*+c-''; c=getchar();}

    return x*f;

}

int n,a,b,c,d,k;

bool notp[N];

ll p[N],mu[N];

void sieve(){

    mu[]=;

    for(int i=;i<=N-;i++){

        if(!notp[i]) p[++p[]]=i,mu[i]=-;

        for(int j=;j<=p[]&&i*p[j]<=N-;j++){

            int t=i*p[j];

            notp[t]=;

            if(i%p[j]==){

                mu[t]=;

                break;

            }

            mu[t]=-mu[i];

        }

    }

    for(int i=;i<=N-;i++) mu[i]+=mu[i-];

}

ll cal(int n,int m,int k){

    n/=k;m/=k;

    if(n>m) swap(n,m);

    ll ans=;int r=;

    for(int i=;i<=n;i=r+){

        r=min(n/(n/i),m/(m/i));

        ans+=(mu[r]-mu[i-])*(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-,d,k)-cal(b,c-,k)+cal(a-,c-,k));

    }

    return ;

}