求

∑i=1n∑j=1md(ij) \sum_{i=1}^{n}\sum_{j=1}^m\mathrm{d}(ij)

其中d(n)\mathrm{d}(n)表示nn的约数个数

多组数据，数据组数≤50000,n,m≤50000\le 50000, n,m\le 50000

Luogu P3327

Solution

首先有一个结论

d(ij)=∑x∣i∑y∣j[gcd(x,y)==1] \mathrm{d}(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]

（这个结论的感性理解见数论函数初探）

直接代入，得到

∑i=1n∑j=1md(ij)=∑i=1n∑j=1m∑x∣i∑y∣j[gcd(x,y)==1]=∑x=1n∑y=1m⌊nx⌋⌊my⌋[gcd(x,y)=1]=∑x=1n∑y=1m⌊nx⌋⌊my⌋∑d∣gcd(x,y)μ(d)=∑d=1min(n,m)μ(d)∑x=1⌊nd⌋⌊ndx⌋∑y=1⌊md⌋⌊mdy⌋ \begin{aligned} &\sum_{i=1}^{n}\sum_{j=1}^{m}\mathrm{d}(ij)\\ =&\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]\\ =&\sum_{x=1}^{n}\sum_{y=1}^m\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor[\gcd(x,y)=1]\\ =&\sum_{x=1}^{n}\sum_{y=1}^m\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor\sum_{d|\gcd(x,y)}\mu(d)\\ =&\sum_{d=1}^{\min(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{dx}\rfloor\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{m}{dy}\rfloor \end{aligned}

后面那一坨，我们可以设f(n)=∑i=1n⌊ni⌋f(n)=\sum_{i=1}^{n}\lfloor\frac{n}{i}\rfloor，那么有

ans=∑d=1min(n,m)μ(d)f(⌊nd⌋)f(⌊md⌋) ans=\sum_{d=1}^{\min(n,m)}\mu(d)f(\lfloor\frac{n}{d}\rfloor)f(\lfloor\frac{m}{d}\rfloor)

那么O(nn)O(n\sqrt n)整除分块预处理f(n)f(n)后，即可单次O(n)O(\sqrt n)整除分块算答案了

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#include <bits/stdc++.h>

#define x first
#define y second
#define y1 Y1
#define y2 Y2
#define mp make_pair
#define pb push_back

using namespace std;

typedef long long LL;
typedef pair<int, int> pii;

template <typename T> inline int Chkmax (T &a, T b) { return a < b ? a = b, 1 : 0; }
template <typename T> inline int Chkmin (T &a, T b) { return a > b ? a = b, 1 : 0; }

inline void proc_status()
{
	ifstream t ("/proc/self/status");
	cerr << string (istreambuf_iterator <char> (t), istreambuf_iterator <char> ()) <<endl;
}

inline int read ()
{
	int sum = 0, fl = 1; char ch = getchar();
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') fl = -1;
	for (; isdigit(ch); ch = getchar()) sum = (sum << 3) + (sum << 1) + ch - '0';
	return sum * fl;
}

const int Maxn = 5e4 + 100;

int N, M;
int Prime[Maxn], Not_Prime[Maxn], prime_cnt, Mu[Maxn], Sum_Mu[Maxn];
LL f[Maxn];

inline void Solve ()
{
	LL ans = 0;
	for (int l = 1, r; l <= min(N, M); l = r + 1)
	{
		r = min(N / (N / l), M / (M / l));
		ans += 1ll * (Sum_Mu[r] - Sum_Mu[l - 1]) * f[N / l] * f[M / l];
	}
	printf("%lld\n", ans);
}

inline void Input ()
{
	N = read(), M = read();
}

inline void Init ()
{
	Mu[1] = Not_Prime[1] = 1;
	for (int i = 2; i <= Maxn - 5; ++i)
	{
		if (!Not_Prime[i])
		{
			Prime[++prime_cnt] = i;
			Mu[i] = -1;
		}
		for (int j = 1; j <= prime_cnt && i * Prime[j] <= Maxn - 5; ++j)
		{
			Not_Prime[i * Prime[j]] = 1;
			if (!(i % Prime[j])) { Mu[i * Prime[j]] = 0; break; }
			else Mu[i * Prime[j]] = -Mu[i];
		}
	}

	for (int i = 1; i <= Maxn - 5; ++i)
		for (int l = 1, r; l <= i; l = r + 1)
		{
			r = i / (i / l);
			f[i] += 1ll * (r - l + 1) * (i / l);
		}

	for (int i = 1; i <= Maxn - 5; ++i) Sum_Mu[i] = Sum_Mu[i - 1] + Mu[i];
}

int main()
{
#ifndef ONLINE_JUDGE
	freopen("A.in", "r", stdin);
	freopen("A.out", "w", stdout);
#endif
	Init();
	int Test_Case = read();
	while (Test_Case--)
	{
		Input();
		Solve();
	}
	return 0;
}