BZOJ 2301: [HAOI2011]Problem b（莫比乌斯反演 + 容斥原理 + 分块优化）

传送门

Problem 2301. – [HAOI2011]Problem b

2301: [HAOI2011]Problem b

Time Limit: 50 Sec  Memory Limit: 256 MB
Submit: 3671  Solved: 1643
[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

解题思路：

首先我们可以看一下这个题目传送门 ，HDU这个题比BZOJ这个题简单，这个题目

首先利用容斥原理将这个答案拆分成四个询问，每次询问都是有多少个数对 (x,y) 满足

1≤x≤N,1≤y≤M 而且 GCD(x,y)==k,这个题目可以和上次一样转化为有

多少个数对 (x,y) 满足 1≤x≤floor(Nk),1≤y≤floor(Mk)，同时满足

GCD(x,y)=1,然后令 f(i) 为 1≤x≤N,1≤y≤M,而且GCD(x,y)==1

的数对个数，令 F(i)是GCD(x,y)==k的倍数的个数，显然有

F(i)=floor(ni)∗floor(mi)

经过莫比乌斯反演之后得到：

f(i)=∑i|dmu(di)∗F(d)=∑i|dmu(di)∗floor(nd)∗floor(md)

然后经过观察发现 floor(nd) 最多有 2n−−√ 个取值，然后枚举除法，写一个莫比乌斯的前缀和就行了。比如 n=100, 那么d 在 [34,50] 之间 nd 都是 2。 那么我们可以把连续的一段 d 一起来算（分块）:设 nd=x ，那么最后一个 nd=x 的 d=nx，所以这段连续的区间就是 [d,n/(n/d)] 结合 m/d，取个最小值就可以了。My Code：

/**2016 - 08 - 17 晚上Author: ITAKMotto:今日的我要超越昨日的我，明日的我要胜过今日的我，以创作出更好的代码为目标，不断地超越自己。**/#include <iostream>#include <cstdio>#include <cstring>#include <cstdlib>#include <cmath>#include <vector>#include <queue>#include <algorithm>#include <set>using namespace std;typedef long long LL;typedef unsigned long long ULL;const int INF = 1e9+5;const int MAXN = 1e5+5;const int MOD = 1e9+7;const double eps = 1e-7;const double PI = acos(-1);using namespace std;LL Scan_LL()///输入外挂{    LL res=0,ch,flag=0;    if((ch=getchar())=='-')        flag=1;    else if(ch>='0'&&ch<='9')        res=ch-'0';    while((ch=getchar())>='0'&&ch<='9')        res=res*10+ch-'0';    return flag?-res:res;}int Scan_Int()///输入外挂{    int res=0,ch,flag=0;    if((ch=getchar())=='-')        flag=1;    else if(ch>='0'&&ch<='9')        res=ch-'0';    while((ch=getchar())>='0'&&ch<='9')        res=res*10+ch-'0';    return flag?-res:res;}void Out(LL a)///输出外挂{    if(a>9)        Out(a/10);    putchar(a%10+'0');}LL mu[MAXN], p[MAXN], k;bool prime[MAXN];void Mobius(){    memset(prime, 0, sizeof(prime));    mu[1] = 1, mu[0] = 0;    k = 0;    for(int i=2; i<MAXN; i++)    {        if(!prime[i])        {            p[k++] = i;            mu[i] = -1;        }        for(int j=0; j<k&&i*p[j]<MAXN; j++)        {            prime[i*p[j]] = 1;            if(i%p[j] == 0)            {                mu[i*p[j]] = 0;                break;            }            mu[i*p[j]] = -mu[i];        }    }    for(int i=1; i<MAXN; i++)        mu[i] += mu[i-1];}int kk;LL Solve(int m, int n){    LL ret = 0;    m /= kk, n /= kk;    int tp = 0;    for(int i=1; i<=m&&i<=n; i=tp+1)    {        tp = min(n/(n/i),m/(m/i));        ret += (mu[tp]-mu[i-1])*(m/i)*(n/i);    }    return ret;}int main(){    Mobius();    int T, a, b, c, d;    scanf("%d",&T);    while(T--)    {        scanf("%d%d%d%d%d",&a,&b,&c,&d,&kk);        LL ans = Solve(b, d) - Solve(a-1, d) + Solve(a-1, c-1) - Solve(b, c-1);        printf("%lld\n",ans);    }    return 0;}