【题目描述】
今天的数学课上,Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数$a$和$b$,$\operatorname{lcm}(a,b)$表示能同时被$a$和$b$整除的最小正整数。例如,$\operatorname{lcm}(6,8)$=$24$。回到家后,Crash还在想着课上学的东西,为了研究最小公倍数,他画了一张$N\times M$的表格。每个格子里写了一个数字,其中第i行第j列的那个格子里写着数为$\operatorname{lcm}(i,j)$。一个$4\times 5$的表格如下:
$\begin{matrix}1&2&3&4&5\2&2&6&4&10\3&6&3&12&15\4&4&12&4&20\ \end{matrix}$
看着这个表格,Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题:这个表格中所有数的和是多少。当N和M很大时,Crash就束手无策了,因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大,Crash只想知道表格里所有数的和$\operatorname{mod} 20101009$的值。

【输入格式】
输入的第一行包含两个正整数,分别表示$N$和$M$。

【输出格式】
输出一个正整数,表示表格中所有数的和$\operatorname{mod} 20101009$的值。

$n \le 10^7$
万幸20101009是个质数

题解

不妨设$n\le m$

那就直接开始推式子吧。。。

懵逼警告

需要求的是$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\operatorname{lcm}(i,j)$

$=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{i*j}{\gcd(i,j)}$

枚举$\gcd(i,j)$的值

$=\sum\limits_{g=1}^{n}\frac{1}{g}*\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j)=g]*i*j$

把后半部分提出来一个$g^2$

$=\sum\limits_{g=1}^{n}g*\sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{g}\rfloor}[\gcd(i,j)=1]*i*j$

套用公式$[n=1]=\sum\limits_{d|n}\mu(d)$

$=\sum\limits_{g=1}^{n}g*\sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{g}\rfloor}i*j*\sum\limits_{d|\gcd(i,j)}\mu(d)$

如法炮制 枚举$d$的值

$=\sum\limits_{g=1}^{n}g*\sum\limits_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)*\sum\limits_{i=1}^{\lfloor\frac{n}{g*d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{g*d}\rfloor}d^2*i*j$

提出$d^2$

$=\sum\limits_{g=1}^{n}g*\sum\limits_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)*d^2*\sum\limits_{i=1}^{\lfloor\frac{n}{g*d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{g*d}\rfloor}i*j$

后面的两个$\sum$可以用等差数列求和的方式表示

$=\sum\limits_{g=1}^{n}g*\sum\limits_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)*d^2*\frac{\lfloor\frac{n}{g*d}\rfloor*(\lfloor\frac{n}{g*d}\rfloor+1)}{2}*\frac{\lfloor\frac{m}{g*d}\rfloor*(\lfloor\frac{m}{g*d}\rfloor+1)}{2}$

好了 这个式子就是最终的形式了

观察一下 枚举$g$我们只需要知道$\frac{n}{g}$及$\frac{m}{g}$是多少,$\sum\limits_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)*d^2*\frac{\lfloor\frac{n}{g*d}\rfloor*(\lfloor\frac{n}{g*d}\rfloor+1)}{2}*\frac{\lfloor\frac{m}{g*d}\rfloor*(\lfloor\frac{m}{g*d}\rfloor+1)}{2}$这一串就能算出来

所以用数论分块来枚举$g$就可以了 记$\lfloor\frac{n}{g}\rfloor =n_0,\lfloor\frac{m}{g}\rfloor =m_0$

同理 只需要知道$\frac{n_0}{d}$和$\frac{m_0}{d}$, 里面那串$\frac{\lfloor\frac{n}{g*d}\rfloor*(\lfloor\frac{n}{g*d}\rfloor+1)}{2}*\frac{\lfloor\frac{m}{g*d}\rfloor*(\lfloor\frac{m}{g*d}\rfloor+1)}{2}$也就能算出来了

所以$d$也可以用数论分块枚举 处理一下$\mu(d)*d^2$这个东西的前缀和就可以了

时间复杂度是$O(\sqrt{n}^2)=O(n)$,很完美

注意处理负数$\mu(d)*d^2$的前缀和会有负数

推导很冗长 代码很简洁

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
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

inline ll read() {
ll x = 0, f = 1; char ch = getchar();
for (; ch > '9' || ch < '0'; ch = getchar()) if (ch == '-') f = -1;
for (; ch <= '9' && ch >= '0'; ch = getchar()) x = (x << 1) + (x << 3) + (ch ^ '0');
return x * f;
}

const ll mod = 20101009;
const ll div2 = mod / 2 + 1;
ll n, m, ans;

bool np[10000005];
ll prime[1000005], mob[10000005], sum[10000005], tot;

inline void getprime() {
mob[1] = 1;
for (ll i = 2; i <= max(n, m); i++) {
if (!np[i]) prime[++tot] = i, mob[i] = -1;
for (ll j = 1; j <= tot && i * prime[j] <= max(n, m); j++) {
np[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mob[i * prime[j]] = 0;
break;
} else mob[i * prime[j]] = -mob[i];
}
}
for (ll i = 1; i <= max(n, m); i++) {
sum[i] = (sum[i-1] + i * i % mod * mob[i] % mod) % mod;
}
}

inline ll calc(ll n0, ll m0) {
ll ret = 0;
for (ll l = 1, r; l <= min(n0, m0); l = r + 1) {
r = min(n0 / (n0 / l), m0 / (m0 / l));
ll s = (sum[r] - sum[l-1] + mod) % mod;
ll s1 = (n0 / l) * (n0 / l + 1) % mod * div2 % mod, s2 = (m0 / l) * (m0 / l + 1) % mod * div2 % mod;
ret = ((ret + s * s1 % mod * s2 % mod) % mod + mod) % mod;
}
return ret;
}

int main() {
n = read(); m = read();
getprime();
for (ll l = 1, r; l <= min(n, m); l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans = (((r - l + 1) * (r + l) % mod * div2 % mod * calc(n / l, m / l) % mod + ans) % mod + mod) % mod;
}
printf("%lld\n", ans);
return 0;
}

评论