BZOJ 2301 [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

题目链接

题意为求出 $T$T 组 ${\sum }_{i=a}^b{\sum }_{i=c}^d\left[gcd\right(i,j)=k]$i=a∑b​i=c∑d​[gcd(i,j)=k]

根据容斥原理可以将算式转化为

$\sum _{i=1}^b\sum _{j=1}^d\left[gcd\right(i,j)=k]-\sum _{i=1}^b\sum _{j=1}^{c-1}\left[gcd\right(i,j)=k]-\sum _{i=1}^{a-1}\sum _{j=1}^d\left[gcd\right(i,j)=k]+\sum _{i=1}^{a-1}\sum _{j=1}^{c-1}\left[gcd\right(i,j)=k]$i=1∑b​j=1∑d​[gcd(i,j)=k]−i=1∑b​j=1∑c−1​[gcd(i,j)=k]−i=1∑a−1​j=1∑d​[gcd(i,j)=k]+i=1∑a−1​j=1∑c−1​[gcd(i,j)=k]

其中每一块算式都为

$\sum _{i=1}^n\sum _{j=1}^m\left[gcd\right(i,j)=k](n\&lt;m)$i=1∑n​j=1∑m​[gcd(i,j)=k](n<m)

化简该式

$\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }\left[gcd\right(i,j)=1]$i=1∑⌊kn​⌋​j=1∑⌊km​⌋​[gcd(i,j)=1]

利用莫比乌斯函数性质可得

$\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }\sum _{d\mid gcd(i,j)}\mu \left(d\right)$i=1∑⌊kn​⌋​j=1∑⌊km​⌋​d∣gcd(i,j)∑​μ(d)

因为 $d\mid gcd(i,j)$d∣gcd(i,j) 所以 $d\mid i$d∣i 且 $d\mid j$d∣j ，考虑枚举d得

$\sum _{d=1}^{\lfloor \frac{n}{k}\rfloor }\mu \left(d\right)\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\left[d\mid i\right]\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }\left[d\mid j\right]$d=1∑⌊kn​⌋​μ(d)i=1∑⌊kn​⌋​[d∣i]j=1∑⌊km​⌋​[d∣j]

其中 ${\sum }_{i=1}^{\lfloor \frac{n}{k}\rfloor }\left[d\mid i\right]$i=1∑⌊kn​⌋​[d∣i] 为 $[1,\lfloor \frac{n}{k}\rfloor ]$[1,⌊kn​⌋] 中 $d$d 的倍数， ${\sum }_{j=1}^{\lfloor \frac{m}{k}\rfloor }\left[d\mid j\right]$j=1∑⌊km​⌋​[d∣j] 为 $[1,\lfloor \frac{m}{k}\rfloor ]$[1,⌊km​⌋] 中 $d$d 的倍数，两项可化简为 $\lfloor \frac{n}{kd}\rfloor$⌊kdn​⌋ 和 $\lfloor \frac{m}{kd}\rfloor$⌊kdm​⌋

这样就可推导出 ${\sum }_{i=1}^n{\sum }_{j=1}^m\left[gcd\right(i,j)=k]={\sum }_{d=1}^{\lfloor \frac{n}{k}\rfloor }\mu \left(d\right)\lfloor \frac{n}{kd}\rfloor \lfloor \frac{m}{kd}\rfloor (n\&lt;m)$i=1∑n​j=1∑m​[gcd(i,j)=k]=d=1∑⌊kn​⌋​μ(d)⌊kdn​⌋⌊kdm​⌋(n<m)

其中莫比乌斯函数 $\mu$μ 使用线性筛筛出，这样算式就可以用一个函数计算（ $N$N 、 $M$M 传参之前已除 $K$K ）

int Solve(int N, int M) {
    if (N > M) swap(N, M);
    int Res = 0;
    for (int d = 1; d <= N; ++d) {
        Res += Mu[d] * (N / d) * (M / d);
    }
    return Res;
}

但是这道题目有 $T$T 组询问，而我们每次四个算式的计算都会有 $O\left(n\right)$O(n) 的复杂度，这样会导致超时

所以需要用到分块来进行优化将代码改为这样

int Solve(int N, int M) {
	if (N > M) swap(N, M);
	int Res = 0;
	for (int d = 1, tp; d <= N; d = tp + 1) {
		tp = min(N / (N / d), M / (M / d));
		Res += (PrefixMu[tp] - PrefixMu[d - 1]) * (N / d) * (M / d);
	}
	return Res;
}

这样每次计算的时间复杂度就为 $O\left(\sqrt{n}\right)$O(n ​)

分块的原理为对于一些 $i$i 其 $\lfloor \frac{n}{i}\rfloor$⌊in​⌋ 的值是相等的（因为下取整），所以可以将其合并以求 ${\sum }_{i=1}^n\lfloor \frac{n}{i}\rfloor$i=1∑n​⌊in​⌋ 那么 $\lfloor \frac{n}{kd}\rfloor \lfloor \frac{m}{kd}\rfloor$⌊kdn​⌋⌊kdm​⌋ 部分便可以使用分块优化，而在分块优化的同时另一部分 $\mu \left(d\right)$μ(d) 需要使用前缀和预处理

分块上界的计算小技巧为 $\lfloor \frac{n}{\lfloor \frac{n}{l}\rfloor }\rfloor$⌊⌊ln​⌋n​⌋

AC代码:

#include <bits/stdc++.h>
using namespace std;

const int maxn = (int)(5e4 + 5);

bool IsPrime[maxn];
int Tot;
int Prime[maxn];
int Mu[maxn];
int PrefixMu[maxn];

void Moblus() {
	for (int i = 0; i < maxn; ++i) IsPrime[i] = true;
	Mu[1] = 1;
	for (int i = 2; i < maxn; ++i) {
		if (IsPrime[i]) {
			Prime[Tot++] = i;
			Mu[i] = -1;
		}
		for (int j = 0; j < Tot && Prime[j] * i < maxn; ++j) {
			IsPrime[i * Prime[j]] = false;
			if (i % Prime[j] == 0) {
				Mu[i * Prime[j]] = 0;
				break;
			}
			Mu[i * Prime[j]] = -Mu[i];
		}
	}
	for (int i = 1; i < maxn; ++i) PrefixMu[i] = PrefixMu[i - 1] + Mu[i];
}

int T;
int A, B, C, D, K;
int Ans;

int Solve(int N, int M) {
	if (N > M) swap(N, M);
	int Res = 0;
	for (int d = 1, tp; d <= N; d = tp + 1) {
		tp = min(N / (N / d), M / (M / d));
		Res += (PrefixMu[tp] - PrefixMu[d - 1]) * (N / d) * (M / d);
	}
	return Res;
}

int main(int argc, char *argv[]) {
	Moblus();
	scanf("%d", &T);
	for (int Case = 1; Case <= T; ++Case) {
		scanf("%d%d%d%d%d", &A, &B, &C, &D, &K);
		Ans = 0;
		Ans += Solve(B / K, D / K);
		Ans -= Solve((B / K), (C - 1) / K);
		Ans -= Solve((A - 1) / K, D / K);
		Ans += Solve((A - 1) / K, (C - 1) / K);
		printf("%d\n", Ans);
	}
	return 0;
}