洛谷P2522思路分享（莫比乌斯反演，数论分块，容斥）

https://www.luogu.com.cn/problem/P2522

题意概述

给定 \(a,b,c,d,k\)，求满足 \(a \le x \le b\)，\(c \le y \le d\)，且 \(\gcd(x,y) = k\) 的数对数量。

多测，\(1\le T \le 5 \times 10^4\)，\(1 \le a \le b \le 5 \times 10^4\)，\(1 \le c \le d \le 5 \times 10^4\)。

思路

考虑容斥。记 \(f(n,m)\) 为满足 \(1\le x \le n\)，\(1 \le y \le m\)，且 \(\gcd(x,y) = k\) 的数对数量，题目所求为 \(f(b,d)-f(a-1,d)-f(b,c-1)+f(a-1,c-1)\)。

接下来求 \(f(n,m)\)，即：

\[\sum_{x=1}^{n}{\sum_{y=1}^{m}{[\gcd(x,y)=k]}} \]

令 \(N = n/k\)，\(M = m/k\)，得：

\[\sum_{x=1}^{N}{\sum_{y=1}^{M}{[\gcd(x,y)=1]}} \]

利用莫比乌斯函数的性质，得：

\[\sum_{x=1}^{N}{\sum_{y=1}^{M}{\sum_{d}{[d\mid x][d\mid y] \mu(d)}}} \]

交换求和，得：

\[\sum_{d}{\mu(d) \left(\sum_{x=1}^{N}{[d \mid x]}\right) \left(\sum_{y=1}^{M}{[d \mid y]}\right) } \]

也就是：

\[\sum_{d}{\mu(d) \left\lfloor \frac{N}{d} \right\rfloor \left\lfloor \frac{M}{d} \right\rfloor} \]

使用数论分块求解即可。

时间复杂度 \(\mathcal{O}(N+T \sqrt{N})\)。\(N\) 为 \(a,b,c,d\) 的上界。

代码

//author:kzssCCC

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

class sieve{
public:
	int n;
	vector<int> prime;
	vector<bool> f;
	vector<int> phi;
	vector<int> mu,premu;
	
	sieve(int _n){
		n = _n;

		phi = vector<int>(n+1);
		phi[1] = 1;

		mu = premu = vector<int>(n+1);
		mu[1] = premu[1] = 1;

		f = vector<bool>(n+1,true);
		f[0] = f[1] = false;
		
		for (int i=2;i<=n;i++){
			if (f[i]){
				prime.push_back(i);
				phi[i] = i-1;
				mu[i] = -1;
			}
			
			for (auto& v:prime){
				if ((ll)v*i>n) break;
				f[v*i] = false;

				if (i%v==0){
					phi[v*i] = phi[i]*v;
					mu[v*i] = 0;
					break;
				}

				phi[v*i] = phi[i]*phi[v];
				mu[v*i] = mu[v]*mu[i];
			}

			premu[i] = premu[i-1]+mu[i];
		}
	}
	
	
	vector<int> factorize(ll x){
		vector<int> res;
		
		for (auto& v:prime){
			if ((ll)v*v>x) break;
			if (x%v==0){
				res.push_back(v);
				while (x%v==0) x/=v;
			}
		}
		if (x>1) res.push_back(x);
		
		return res;
	}
	
	vector<pair<int,int>> factorize_cnt(ll x){
		vector<pair<int,int>> res;
		
		for (auto& v:prime){
			if ((ll)v*v>x) break;
			if (x%v==0){
				pair<int,int> pr{v,0};
				while (x%v==0){
					pr.second++;
					x/=v;
				}
				res.push_back(pr);
			}
		}
		if (x>1) res.emplace_back(x,1);
		return res;
	}
};

sieve se(5e4+5);

void solve(){
	int a,b,c,d,k;
	cin >> a >> b >> c >> d >> k;

	auto f = [&](int n,int m){
		n /= k;
		m /= k;
		ll res = 0;

		int l = 1;
		while (l<=min(n,m)){
			int r = min({min(n,m),n/(n/l),m/(m/l)});

			res += (ll)(n/l)*(m/l)*(se.premu[r]-se.premu[l-1]);
			l = r+1;
		}

		return res;
	};

	cout << f(b,d)-f(a-1,d)-f(b,c-1)+f(a-1,c-1) << '\n';
}

int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);
	
	int t;
	cin >> t;
	while (t--) solve();

	return 0;
}