显然题目就是要我们求 \(\sum\limits_{i = 1}^n \sum\limits_{j = 1}^m \text{lcm}(i,j)\)。
根据小学奥数，就是 \(\sum\limits_{i = 1}^n \sum\limits_{j = 1}^m \dfrac{ij}{\gcd(i,j)}\)。
\[ \begin{aligned} & \sum\limits_{i = 1}^n \sum\limits_{j = 1}^m \dfrac{ij}{\gcd(i,j)} \\ = & \sum\limits_{i = 1}^{\min(n,m)} i \sum\limits_{xi \le \frac n i} \sum\limits_{yi \le \frac m i} [\gcd(x,y) = 1]xy \end{aligned} \]
敏锐地察觉到 \([n = 1]\) 这个东西可以替换成 \(\epsilon(n) = \sum\limits_{d | n} \mu(d)\)。
那么有
\[ \begin{aligned} & \sum\limits_{i = 1}^{\min(n,m)} i \sum\limits_{xi \le \frac n i} \sum\limits_{yi \le \frac m i} [\gcd(x,y) = 1]xy \\ = & \sum\limits_{i = 1}^{\min(n,m)} i \sum\limits_{xi \le \frac n i} \sum\limits_{yi \le \frac m i} xy \sum\limits_{d | \gcd(x,y)} \mu(d) \end{aligned} \]
根据套路，接着把 \(d\) 提出来。
\[ \begin{aligned} & \sum\limits_{i = 1}^{\min(n,m)} i \sum\limits_{xi \le \frac n i} \sum\limits_{yi \le \frac m i} xy \sum\limits_{d | \gcd(x,y)} \mu(d) \\ = & \sum\limits_{i = 1}^{\min(n,m)} i \sum\limits_{d = 1}^{\min(\lfloor \frac n i \rfloor,\lfloor \frac m i \rfloor)} \mu(d) d^2 \sum\limits_{xd \le \frac n i} \sum\limits_{yd \le \frac m i} xy \\ = & \sum\limits_{i = 1}^{\min(n,m)} i \sum\limits_{d = 1}^{\min(\lfloor \frac n i \rfloor,\lfloor \frac m i \rfloor)} \mu(d) d^2 \text{sum}(\lfloor \dfrac n {id} \rfloor,\lfloor \dfrac m {id} \rfloor) \end{aligned} \]
其中 \(\text{sum}(x,y) = \sum\limits_{i = 1}^x \sum\limits_{j = 1}^y ij = \dfrac{xy(x + 1)(y + 1)} 4\)。
那么预处理 \(\mu(i)i^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
using namespace std;
const int N = 1e7;
const long long mod = 20101009;
int n,m;
int cnt,prime[N + 5],vis[N + 5],sum[N + 5];
long long ans;
long long get_sum(int n,int m)
{
long long ret = 0;
if(n > m)
swap(n,m);
for(register int l = 1,r;l <= n;l = r + 1)
{
r = min(n / (n / l),m / (m / l));
ret = (ret + (long long)(sum[r] - sum[l - 1] + mod) * ((long long)(n / l) * (n / l + 1) / 2 % mod) % mod * ((long long)(m / l) * (m / l + 1) / 2 % mod) % mod);
}
return ret;
}
int main()
{
sum[1] = 1;
for(register int i = 2;i <= N;++i)
{
if(!vis[i])
prime[++cnt] = i,sum[i] = -1;
for(register int j = 1;j <= cnt && i * prime[j] <= N;++j)
{
vis[i * prime[j]] = 1;
if(!(i % prime[j]))
break;
sum[i * prime[j]] = -sum[i];
}
}
for(register int i = 1;i <= N;++i)
sum[i] = (sum[i - 1] + (long long)i * i % mod * (sum[i] + mod)) % mod;
scanf("%d%d",&n,&m);
if(n > m)
swap(n,m);
for(register int l = 1,r;l <= n;l = r + 1)
{
r = min(n / (n / l),m / (m / l));
ans = (ans + (long long)(l + r) * (r - l + 1) / 2 % mod * get_sum(n / l,m / l) % mod) % mod;
}
printf("%lld\n",ans);
}