BZOJ2301[HAOI2011]problem b

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

---

算是写的第一道反演题。题目要求的是一段区间，通常的做法是求出1到右边界，然后再用容斥来求解区间。

我们考虑如何求1到右边界，设右边界为$n$和$m$，则要求的式子是$\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=k]$，我们来推下式子。

$\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=k]=\sum_{i=1}^{\frac{n}{k}}\sum_{j=1}^{\frac{m}{k}}[gcd(i,j)=1]=\sum_{i=1}^{\frac{n}{k}}\sum_{j=1}^{\frac{m}{k}}\sum_{x|gcd(i,j)}\mu(x)$(根据性质1)

接下来要做最重要的一步转换，大家要好好理解：$\sum_{i=1}^{\frac{n}{k}}\sum_{j=1}^{\frac{m}{k}}\sum_{x|gcd(i,j)}\mu(x)=\sum_{x=1}^{\frac{min(n,m)}{k}}\mu(x) [\frac{n}{kx}][\frac{m}{kx}]$

这步大家可能不太理解，我刚看的时候也是一脸懵逼，简单的说就是将$x$提到外面，而右边$[\frac{n}{kx}][\frac{m}{kx}]$表示的是$x$的倍数的个数。

我们发现$[\frac{n}{kx}][\frac{m}{kx}]$可以用数论分块，所以直接分块枚举x即可。

---

这道题还有另外一种推法，需要构造。

我们令$f(x)=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=x]$,$F(x)=\sum_{i=1}^{n}\sum_{j=1}^{m}[x|gcd(i,j)]$

那么显然$F(x)=\sum_{x|d}f(d)$,那么根据莫比乌斯反演$f(x)=\sum_{x|d}\mu(\frac{d}{x})F(d)$

我们考虑$F(x)$是什么,$F(x)$必须要求$i,j$都是$x$的倍数，所以$F(x)=[\frac{n}{x}][\frac{m}{x}]$

所以答案就是$f(k)=\sum_{k|x}\mu(\frac{x}{k})[\frac{n}{x}][\frac{m}{x}]=\sum_{x=1}^{\frac{min(n,m)}{k}}\mu(x) [\frac{n}{kx}][\frac{m}{kx}]$

与上面的做法殊途同归。

#include<bits/stdc++.h>
using namespace std;
int read(){
    char c;int x;while(c=getchar(),c<'0'||c>'9');x=c-'0';
    while(c=getchar(),c>='0'&&c<='9') x=x*10+c-'0';return x;
}
void print(int x){
    if(x/10) print(x/10);
    putchar(x%10+'0');
}
int T,a,b,c,d,k,top,ans,mu[50005],sum[50005],pri[50005],vis[50005];
int calc(int n,int m){
    if(n>m) swap(n,m);
    int l=1,r=0,res=0;n=n/k;m=m/k;
    while(l<=n){
        r=min(n/(n/l),m/(m/l));r=min(r,n);
        res+=(n/l)*(m/l)*(sum[r]-sum[l-1]);
        l=r+1;
    }
    return res;
}
int main()
{
    T=read();mu[1]=sum[1]=1;
    for(int i=2;i<=50000;i++){
        if(!vis[i]) pri[++top]=i,mu[i]=-1;
        for(int j=1;j<=top&&pri[j]*i<=50000;j++){
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) break;
            mu[i*pri[j]]=-mu[i];
        }
        sum[i]=sum[i-1]+mu[i];
    }
    while(T--){
        a=read();b=read();c=read();d=read();k=read();
        ans=calc(b,d)-calc(a-1,d)-calc(b,c-1)+calc(a-1,c-1);
        print(ans);puts("");
    }
    return 0;
}