hdu4746 Mophues

莫反+数论分块+前缀和

4.HDU4746 Mophues

题目描述

给$Q(1\le Q\le 5000)$个询问，每个询问给出$n,m,p$求$\sum\limits_{i=1}^N\sum\limits_{j=1}^M[h(gcd(i,j))<=p]$。

题目分析

因为把数据范围内的数字分解，p最多为19，那么$18\le p$时直接输出$n*m$，对于另一种情况，令$f(x)=\sum\limits_{i=1}^N\sum\limits_{j=1}^M[gcd(i,j)=x]$,那么答案就是$\sum\limits_{1\le x\le n,h(x)\le p}f(x)$,要求$f(x)$，考虑莫比乌斯反演，定义$g(x)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^m[x|gcd(i,j)]$，易得所有的$i,j$为x的倍数，那就很容易求出来$g(x)=\frac{n}{x}\frac{m}{x}$，同时有$g(x)=\sum\limits_{x|y}f(y)$,通过倍数形式的莫反得到$f(x)=\sum\limits_{x|y}\mu(y/x)*g(y)=\sum\limits_{x|y}\mu(y/x)(n/y)*(m/y)$.那么答案为

$\sum\limits_{1\le x\le n,h(x)\le p}\sum\limits_{x|y}\mu(y/x)(n/y)(m/y) \\=\sum\limits_{1\le x\le n,h(x)\le p}\ \sum\limits_{1\le k\le n/x}\mu(k)(n/xk)(m/xk) \\=\sum\limits_{1\le x\le n,h(x)\le p}\ \sum\limits_{1\le k\le n/x}\mu(k)(n/xk)(m/xk) \\交换求和顺序 \\=\sum\limits_{1\le T\le n}(n/T)(m/T) \sum\limits_{x|T}\mu(T/x)[h(x)\le p]$

然后通过一些妙哇的筛法求出来后面的一坨东西。

代码

//
// Created by mrx on 2022/9/24.
//

#include <vector>
#include <algorithm>
#include <iostream>
#include <array>

using ll = long long;

int main() {
	std::ios::sync_with_stdio(false);
	std::cin.tie(nullptr);
	std::cout.tie(nullptr);

	std::vector<int> prime, vis(5e5 + 10), mu(5e5 + 10), h(5e5 + 10);
	mu[1] = 1, h[1] = 0;
	for (int i = 2; i <= 5e5; ++i) {
		if (!vis[i])prime.push_back(i), mu[i] = -1, h[i] = 1;
		for (int j = 0; j < prime.size() && i * prime[j] <= 5e5; ++j) {
			vis[i * prime[j]] = 1;
			h[prime[j] * i] = h[i] + 1;
			if (i % prime[j]) {
				mu[i * prime[j]] = -mu[i];
			} else {
				mu[i * prime[j]] = 0;
				break;
			}
		}
	}

	std::vector<std::array<ll, 20>> pref(5e5 + 10);
	for (int i = 1; i <= 5e5; ++i) {
		for (int k = 1; i * k <= 5e5; k++) {
			pref[i * k][h[i]] += mu[k];
		}
	}

	for (int i = 0; i <= 5e5; ++i) for (int j = 1; j < 20; ++j)pref[i][j] += pref[i][j - 1];
	for (int i = 1; i < 5e5; ++i)for (int j = 0; j < 20; ++j)pref[i][j] += pref[i - 1][j];

	int q;
	std::cin >> q;
	while (q--) {
		ll n, m, p;
		std::cin >> n >> m >> p;
		if (p >= 18) {
			std::cout << n * m << '\n';
		} else {
			if (n > m)std::swap(n, m);
			ll ans = 0;
			for (ll l = 1, r; l <= n; l = r + 1) {
				r = std::min(n / (n / l), m / (m / l));
				ans += (pref[r][p] - pref[l - 1][p]) * (n / l) * (m / l);
			}
			std::cout << ans << '\n';

		}
	}
	return 0;
}