[bzoj2310][HAOI2011]Problem b

2301: [HAOI2011]Problem b

Time Limit: 50 Sec Memory Limit: 256 MB

对于给出的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

发现这题很gcd很反演，所以我们用反演来做。

题目要求输入a,b,c,d,k求
$\sum_a^b\sum_c^dgcd(x,y)=k$

=$\sum_{a/k}^{b/k}\sum_{c/k}^{d/k}[gcd(x,y)=1]$

=$\sum_{a/k}^{b/k}\sum_{c/k}^{d/k}\sum_{z|gcd(x,y)}\mu(z)$
所以发现如果gcd(x,y)为z的倍数的时候产生$\mu(z)$的贡献。
改变枚举顺序，在这里的时候我们先把a-1,c-1,原因等会就明白了。。

=$\sum_z^{min(n,m)}(\lfloor \frac{b}{z} \rfloor-\lfloor \frac{a}{z} \rfloor)*(\lfloor \frac{d}{z} \rfloor-\lfloor \frac{c}{z} \rfloor)*\mu(z)$

发现直接写肯定是T飞了。所以分段求和优化一下变成了数据组数*sqrt的复杂度

#include<cstdio>
#include<algorithm>
#include<string>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<iostream>
#define ll long long
using namespace std;
const int N=51000;
int n,m,a,b,c,d,k,j;
int mu[N],prime[N];
bool is[N];
ll ans;
inline int read()
{
    char c;
    int res,flag=0;
    while((c=getchar())>'9'||c<'0') if(c=='-')flag=1;
    res=c-'0';
    while((c=getchar())>='0'&&c<='9') res=(res<<3)+(res<<1)+c-'0';
    return flag?-res:res;
}
inline int get(int x,int y)
{
    return !(x%y)?x/y:x/y+1;
}
const int inf=2147483646;
inline void solve()
{
    ans=0;
    a=read();b=read();
    c=read();d=read();k=read();
    a=get(a,k)-1;b=b/k;c=get(c,k)-1;d=d/k;
    if(b>d) swap(b,d),swap(a,c);
    for(int i=1;i<=b;i=j+1)
    {
        j=min(b/(b/i),min((a/i)?a/(a/i):inf,min((c/i)?c/(c/i):inf,d/(d/i))));
        ans+=(ll)(b/i-a/i)*(d/i-c/i)*(mu[j]-mu[i-1]);
    }

/*  for(int i=1;i<=b;++i)
        ans+=(ll)(b/i-a/i)*(d/i-c/i)*(mu[i]-mu[i-1]);*/
    printf("%lld\n",ans);
}
int main()
{
//  freopen("2301.txt","r",stdin);
    int T=read();
    mu[1]=1;
    for(int i=2;i<N;++i)
    {
        if(!is[i]) mu[i]=-1,prime[++prime[0]]=i;
        for(j=1;j<=prime[0]&&prime[j]*i<N;++j)
        {
            is[prime[j]*i]=1;
            if(!(i%prime[j]))
            {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
        mu[i]+=mu[i-1];
    }
    while(T--)
        solve();
}