【bzoj2301】[HAOI2011]Problem b

2301: [HAOI2011]Problem b

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

Source

[ Submit][ Status][ Discuss]



心塞。。

这题就是莫比乌斯反演

我们设g(i)为gcd(x,y)为i 倍数的对数

那么显然有

g(i) = ([b / i] - [(a - 1) / i]) * ([d / i] - [(c - 1) / i])

再设f(i) 为gcd(x,y)为i的对数

显然有

f[i] = Σμ(D) * g(D / i)        (d为i的倍数)

这是一个O (n ^ 2)的算法，显然会超时

那么我们考虑优化

展开一下f[i]

有

f[i] = Σμ(D) *  ([b / D] - [(a - 1) / D]) * ([d / D] - [(c - 1) / D])

     = Σμ(D) *  ([b / D] * [d / D] - [(a - 1) / D] * [d / D] - [b / D] * [(c - 1) / D] + [(a - 1) / D] * [(c - 1) / D])

会发现有一堆向下取整的东西

然后联想一下向下取整的性质

如果说对于一个[n / i]来说

对于固定的n来说，[n / i]的个数不会超过2 * n ^ 0.5个

那么[n / i] * [m / i]的个数不会超过2 * n ^ 0.5 + 2 * m ^ 0.5个

于是可以把[n / i] * [m / i]相同的部分利用前缀和合在一起算

值得注意的是，如何求出[n / i]相同部分的最右端呢？

我们考虑当前左端点i,则n / [n / i]就是符合条件的最右端

然后一通瞎搞就好了

代码：

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<cstdlib>
#include<vector>
#include<queue>
#include<stack>
using namespace std;

typedef long long LL;

const LL INF = 1000000000000000000LL;
const int maxn = 50100;

LL n,m,a,b,c,d,k,q;
LL sum[maxn],prime[maxn],tot,mu[maxn];
bool mark[maxn];

inline LL getint()
{
	LL ret = 0,f = 1;
	char c = getchar();
	while (c < '0' || c > '9')
	{
		if (c == '-') f = -1;
		c = getchar();
	}
	while (c >= '0' && c <= '9')
		ret = ret * 10 + c - '0',c = getchar();
	return ret * f;
}

inline LL cal(LL n,LL m,LL k)
{
	n /= k; m /= k;
	LL last,ret = 0;
	for (int i = 1; i <= min(n,m); i = last + 1)
	{
		last = min(n / (n / i),m / (m / i));
		ret += (n / i) * (m / i) * (sum[last] - sum[i - 1]);
	}
	return ret;
}

int main()
{
	mu[1] = 1;
	for (int i = 2; i <= maxn; i++)
	{
		if (!mark[i])
		{
			mu[i] = -1;
			prime[++tot] = i;
		}
		for (int j = 1; j <= tot; j++)
		{
			if (i * prime[j] > maxn) break;
			mark[i * prime[j]] = 1;
			if (i % prime[j]) mu[i * prime[j]] = -mu[i];
			else
			{
				mu[i * prime[j]] = 0;
				break;
			}
		}
	}
	for (int i = 1; i <= maxn; i++) sum[i] = sum[i - 1] + mu[i];
	q = getint();
	while (q--)
	{
		LL ans = 0;
		a = getint(); b = getint(); c = getint(); d = getint(); k = getint();
		printf("%lld\n",cal(b,d,k) - cal(a - 1,d,k) - cal(b,c - 1,k) + cal(a - 1,c - 1,k));
	}
	return 0;
}