求

∑i=1n∑j=1nlcm(i,j) \sum_{i=1}^{n}\sum_{j=1}^{n}\mathrm{lcm}(i,j)

n≤1010n\le 10^{10}

51nod 1238

Solution

注：以下所有（其实也包括之前的几篇文章）S(n)=∑i=1niS(n)=\sum_{i=1}^n i

∑i=1n∑j=1nlcm(i,j)=∑i=1n∑j=1nijgcd(i,j)=∑d=1n1d∑i=1n∑j=1nij[gcd(i,j)=d]=∑d=1nd∑i=1⌊nd⌋∑j=1⌊nd⌋ij[gcd(i,j)=1] \begin{aligned} &\sum_{i=1}^{n}\sum_{j=1}^{n}\mathrm{lcm}(i,j)\\ =&\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{ij}{\gcd(i,j)}\\ =&\sum_{d=1}^{n}\frac{1}{d}\sum_{i=1}^{n}\sum_{j=1}^{n}ij[\gcd(i,j)=d]\\ =&\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij[\gcd(i,j)=1] \end{aligned}

化到这一步之后，我硬刚用μ\mu推了两天时间，尝试了各种方法，最终也只能推出来各种比线性低一点点的玄学复杂度的做法，根本跑不过

因此事实证明，这道题必须用φ\varphi做，刚好也算是积累了一个φ\varphi的套路

这个套路长这样：

∑i=1ni[gcd(i,n)=1]=φ(n)n+[n=1]2 \sum_{i=1}^{n}i[\gcd(i,n)=1]=\frac{\varphi(n)n+[n=1]}{2}

感性证明：

不难发现，若gcd(i,n)=1\gcd(i,n)=1，则gcd(n−i,n)=1\gcd(n-i,n)=1（若n−in-i与nn有公共因子，则n−(n−i)=in-(n-i)=i与nn也有公共因子，矛盾）

那么任何一对这样的数对答案的贡献就是nn，而总共有φ2\frac{\varphi}{2}对，因此得证

Update 12.23

来一发理性证明

∑i=1ni[gcd(i,n)=1]=∑d∣nμ(d)∑i=1ndid=∑d∣nμ(d)n(nd+1)2=n2(∑d∣nμ(d)nd+∑d∣nμ(d))=n2([n==1]+φ(n)) \begin{aligned} &\sum_{i=1}^{n}i[gcd(i,n)=1]\\ =&\sum_{d|n}\mu(d)\sum_{i=1}^{\frac{n}{d}}id\\ =&\sum_{d|n}\mu(d)\frac{n(\frac{n}{d}+1)}{2}\\ =&\frac{n}{2}(\sum_{d|n}\mu(d)\frac{n}{d}+\sum_{d|n}\mu(d))\\ =&\frac{n}2([n==1]+\varphi(n)) \end{aligned}

有了这个套路，就可以接着上一步化式子了

ans=∑d=1nd{(2∑i=1⌊nd⌋i∑j=1ij[gcd(i,j)=1])−1}=∑d=1nd{(2∑i=1⌊nd⌋i∗i∗φ(i)+[i=1]2)−1}=∑d=1nd∑i=1⌊nd⌋i2φ(i)=∑i=1nφ(i)i2S(⌊ni⌋) \begin{aligned} ans&=\sum_{d=1}^{n}d\left\{(2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i\sum_{j=1}^{i}j[\gcd(i,j)=1])-1 \right\}\\ &=\sum_{d=1}^{n}d\left\{(2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i*\frac{i*\varphi(i)+[i=1]}{2})-1 \right\}\\ &=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i^2\varphi(i)\\ &=\sum_{i=1}^{n}\varphi(i)i^2S(\lfloor\frac{n}{i}\rfloor) \end{aligned}

tmd就和这道题几乎一模一样了

我tm还傻逼地用μ\mu搞了两天

Summary

这一道题带来收获还挺多的

首先是那个关于φ\varphi的套路：∑i=1ni[gcd(i,n)=1]=φ(n)n+[n=1]2\sum_{i=1}^{n}i[\gcd(i,n)=1]=\frac{\varphi(n)n+[n=1]}{2}

其次就是在推假的μ\mu做法的时候，发现的一些东西：

对于一个这样的式子∑d=1n...∑g=1⌊nd⌋...\sum_{d=1}^{n}...\sum_{g=1}^{\lfloor\frac{n}{d}\rfloor}...

在化式子的时候，既可以考虑把枚举倍数转化成枚举约数，也可以直接交换dd和gg的位置

这就取决于里面的东西是与dgdg相关的比较多还是与dd或gg相关的比较多

还有一个小收获是证明了个(Id2μ)∗Id2=Id(\mathrm{Id}^2\mu)*\mathrm{Id}^2=\mathrm{Id}，不知道以后有没有可能用得到

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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>

#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 LL read ()
{
	LL 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 = 1e7 + 100, Mod = 1e9 + 7;

LL N, inv2, inv6;
int Phi[Maxn], Not_Prime[Maxn], Prime[Maxn], prime_cnt;
__gnu_pbds :: gp_hash_table<LL, int> F;
LL FF[Maxn];

inline void Add (LL &a, LL b) { a += b; if (a >= Mod) a -= Mod; }

inline LL Pow (LL a, int b)
{
	LL ans = 1;
	while (b) { if (b & 1) ans = 1ll * ans * a % Mod; a = 1ll * a * a % Mod; b >>= 1; }
	return ans;
}

inline LL Calc_S (LL x)
{
	x %= Mod;
	LL sum = 1ll * x * (x + 1) % Mod * inv2 % Mod;
	return sum;
}

inline LL Calc_G (LL x)
{
	x %= Mod;
	return x * (x + 1) % Mod * (2 * x + 1) % Mod * inv6 % Mod;
}

inline LL Calc_F (LL n)
{
	if (n <= Maxn - 5) return FF[n];
	if (F[n]) return F[n];
	LL ans = 1ll * Calc_S(n) * Calc_S(n) % Mod, sum1, sum2, sum;
	for (LL l = 2, r; l <= n; l = r + 1)
	{
		r = n / (n / l);
		sum1 = Calc_G(l - 1), sum2 = Calc_G(r);
		sum = Calc_F(n / l);
		Add (ans, Mod - sum * (sum2 - sum1 + Mod) % Mod);
	}
	return F[n] = ans;
}

inline void Solve ()
{
	LL sum1, sum2, sum, ans = 0;
	for (LL l = 1, r; l <= N; l = r + 1)
	{
		r = N / (N / l);
		sum1 = Calc_F(l - 1), sum2 = Calc_F(r);
		sum = Calc_S(N / l);
		Add (ans, sum * (sum2 - sum1 + Mod) % Mod);
	}
	cout<<ans<<endl;
}

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

inline void Init ()
{
	inv2 = Pow(2, Mod - 2), inv6 = Pow(6, Mod - 2);

	Phi[1] = 1;
	for (int i = 2; i <= min(N, 1ll * Maxn - 5); ++i)
	{
		if (!Not_Prime[i])
		{
			Prime[++prime_cnt] = i;
			Phi[i] = i - 1;
		}
		for (int j = 1; j <= prime_cnt && Prime[j] * i <= min(N, 1ll * Maxn - 5); ++j)
		{
			Not_Prime[i * Prime[j]] = 1;
			if (!(i % Prime[j])) { Phi[i * Prime[j]] = Phi[i] * Prime[j]; break; }
			else Phi[i * Prime[j]] = Phi[i] * (Prime[j] - 1);
		}
	}
	
	for (LL i = 1; i <= min(N, 1ll * Maxn - 5); ++i) FF[i] = Phi[i] * i % Mod * i % Mod;
	for (LL i = 1; i <= min(N, 1ll * Maxn - 5); ++i) Add (FF[i], FF[i - 1]);
}

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