BZOJ2301 [HAOI2011]Problem b

本文版权归ljh2000和博客园共有，欢迎转载，但须保留此声明，并给出原文链接，谢谢合作。

  

本文作者：ljh2000
作者博客：http://www.cnblogs.com/ljh2000-jump/
转载请注明出处，侵权必究，保留最终解释权！

 

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

 

 

正解：线性筛+莫比乌斯函数+分块

解题报告：

　　 这道题是一道莫比乌斯函数的入门题，我从昨晚一直看题解和分析，直到今天才完全弄懂。系统性地写一份解题报告吧。

　　 ps：我在学习的时候参考了这些论文与博客，贴出来以供参考吧：

　　 http://wenku.baidu.com/view/fbe263d384254b35eefd34eb.html

　　 http://wenku.baidu.com/view/fbec9c63ba1aa8114431d9ac.html

　　 http://www.cnblogs.com/zhsl/p/3269288.html

 

　　 可以发现题目要求的是∑ ∑ [gcd(i,j)==k] (1<=i<=n && 1<=j<=m)，这样一个东西。（至于区间的话我们很快可以发现根据容斥原理[a,b]和[c,d]组合起来的答案就等于[1,b]组合[1,d] - [1,a-1]组合[1][d] - [1][b]组合[1][c-1] + [1][a-1] 组合 [1][c-1]，问题从而转化。）但是直接求的复杂度很高，还是多组询问，肯定会TLE。

　　题目所求的可以变成

　　∑ ∑ [gcd(i,j)==1] (1<=i<=n/k && 1<=j<=m/k) ，相当于是同时除以k；

　　考虑再把式子变形：

　　∑ ∑ ∑ μ(d) (1<=i<=n && 1<=j<=m && d|gcd(i,j))；

　　

      根据莫比乌斯函数的性质，当n=1时，∑μ(i)=1(1<=i<=n)；否则∑μ(i)=0(1<=i<=n)，所以实际上只有在gcd(i,j)==1时上式才有贡献。可是复杂度还是没有改变啊。我看了一篇解题报告，里面对于这个问题讲的很详细，大概是说根据容斥原理，可以把1<=i<=n && 1<=j<=m且满足gcd(i,j)==k的方案数转化成1<=i<=n && 1<=j<=m且满足gcd(i,j)>=k的方案数。二者的联系就在于莫比乌斯函数。论文里面画了一个图，一目了然地展示了二者的联系。

 

　　可以根据莫比乌斯函数的性质很快反应过来二者的联系。

 　  所以得到新式：

    　　 ans=∑μ(i)*[n/i]*[m/i] (1<=i<=min(n,m),[n/i]表示n/i下取整)

　　容易发现这个式子可以使我的复杂度降低到O(n)，但是询问次数也是5w级别，还是会TLE。考虑继续优化上式，考虑到[n/i]、[m/i]都会有大量的完全相等的部分，我们可以把[n/i]、[m/i]都相等的部分放在一起算，也就是一个分块的思想。容易发值为[n/i]的起点是i，终点是n/([n/i])，而取值不会超过根号n种，也就是说我只要算根号次。但是因为是两个都要相等，在数轴上画线段可以发现，这样的话总取值不会超过2×根号n种（论文里面也有图说明）。所以复杂度降低到可以接受的范围内。

　　值得一提的是，由于我们是分块计算的，所以我们必须要预处理出莫比乌斯函数和他的前缀和。预处理的方式是根据莫比乌斯函数的性质，进行一次线性筛。在筛出质数的同时随便求一下莫比乌斯函数，比较巧妙，然后直接构前缀和就可以了。

 

//It is made by ljh2000
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <ctime>
#include <vector>
#include <queue>
#include <map>
#include <set>
using namespace std;
typedef long long LL;
const int inf = (1<<30);
const int MAXN = 50011;
const int MAXL = 50000;
LL a,b,c,d,k,ans;
int prime[MAXN],cnt;
bool ok[MAXN];
int mobius[MAXN],sum[MAXN];//莫比乌斯函数及其前缀和

inline int getint()
{
    int w=0,q=0; char c=getchar();
    while((c<'0' || c>'9') && c!='-') c=getchar(); if(c=='-') q=1,c=getchar(); 
    while (c>='0' && c<='9') w=w*10+c-'0', c=getchar(); return q ? -w : w;
}

inline void init(){//递推得到莫比乌斯函数
    //1的莫比乌斯函数是1；质数的莫比乌斯函数为1；含有相同质因子的数莫比乌斯函数为0；
    //不含有相同质因子的数的莫比乌斯函数函数为（-1）^k，k为质因子个数
    mobius[1]=1;
    for(int i=2;i<=MAXL;i++) {
    if(!ok[i]) {
        mobius[i]=-1;
        prime[++cnt]=i;
    }
    for(int j=1;j<=cnt && i*prime[j]<=MAXL;j++) {
        ok[i*prime[j]]=1;//标记合数
        if(i%prime[j]) mobius[i*prime[j]]=-mobius[i];//互质的两个数乘起来得到一个不含有相同质因子的数，质因子个数奇偶性改变，莫比乌斯函数变号
        else { mobius[i*prime[j]]=0; break; }//留到后面再筛，此处已经可以break
    }
    }
    for(int i=1;i<=MAXL;i++) sum[i]=sum[i-1]+mobius[i];//求莫比乌斯函数的前缀和
}

inline LL solve(LL n,LL m){//计算a在[1,n]且b在[1,m]中的gcd(a,b)==1的数目
    n/=k; m/=k;
    if(n>m) swap(n,m); if(n==0) return 0;
    LL i,next1,next2,next;//把相等的部分直接分块一起计算
    LL tot=0;
    for(i=1;i<=n;i=next) {
    next1=n/(n/i); next2=m/(m/i);
    next=min(next1,next2);
    tot+=(n/i)*(m/i)*(sum[next]-sum[i-1]);
    next++;
    }
    return tot;
}

inline void work(){
    int T=getint(); init();
    while(T--) {
    a=getint(); b=getint(); c=getint(); d=getint(); k=getint();
    ans=solve(b,d)-solve(a-1,d)-solve(b,c-1)+solve(a-1,c-1);//容斥原理
    printf("%lld\n",ans);
    }
}

int main()
{
    work();
    return 0;
}