求
\[\sum^{n}_{i = 1}\sum^{m}_{j = 1}lcm(i, j) \]即
\[\sum^{n}_{i = 1}\sum^{m}_{j = 1}\dfrac{ij}{\gcd(i, j)} \]即
\[\sum^{\min(n, m)}_{k = 1}\sum^{n}_{i = 1}\sum^{m}_{j = 1}[\gcd(i,j) = k]\times\dfrac{ij}{k} \]即
\[\begin{aligned} &\sum^{\min(n, m)}_{k = 1}\sum^{\small[\dfrac{n}{k}]}_{i = 1}\sum^{\small[\dfrac{m}{k}]}_{j = 1}[\gcd(i,j) = 1]\dfrac{ik\times jk}{k}\\ &\sum^{\min(n, m)}_{k = 1}k\sum^{\small[\dfrac{n}{k}]}_{i = 1}\sum^{\small[\dfrac{m}{k}]}_{j = 1}[\gcd(i,j) = 1]\times ij\\ \end{aligned} \]设 \(H(n,m) = \sum^{n}_{i = 1}\sum^{m}_{j = 1}[\gcd(i,j) = 1]\times ij\).
即
\[\sum^{\min(n, m)}_{k = 1}k \times H([\dfrac{n}{k}], [\dfrac{m}{k}]) \]数论分块即可。
接下来求解 \(H(n,m) = \sum^{n}_{i = 1}\sum^{m}_{j = 1}[\gcd(i,j) = 1]\times ij\).
\[\begin{aligned} H(n,m) &= \sum^{n}_{i = 1}\sum^{m}_{j = 1}[\gcd(i,j) = 1]\times ij\\ &=\sum^{n}_{i = 1}\sum^{m}_{j = 1}\sum^{}_{d \mid \gcd(i,j)}\mu(d)\times ij\\ &=\sum^{min(n,m)}_{d}\sum^{n}_{i = 1}\sum^{m}_{j = 1}\mu(d)\times ij [d \mid \gcd(i,j)]\\ &=\sum^{min(n,m)}_{d}\sum^{\small\dfrac{n}{d}}_{i = 1}\sum^{\small\dfrac{m}{d}}_{j = 1}\mu(d)\times di \times dj\\ &=\sum^{min(n,m)}_{d}\mu(d)d^2 \times\sum^{\small\dfrac{n}{d}}_{i = 1}\sum^{\small\dfrac{m}{d}}_{j = 1}ij\\ &=\sum^{min(n,m)}_{d}\mu(d)d^2 \times\sum^{\small\dfrac{n}{d}}_{i = 1}i\sum^{\small\dfrac{m}{d}}_{j = 1}j\\ &=\sum^{min(n,m)}_{d}\mu(d)d^2 \times\dfrac{\small\dfrac{n}{d}(\dfrac{n}{d} + 1)}{2} \times \dfrac{\small\dfrac{m}{d}(\dfrac{m}{d} + 1)}{2}\\ \end{aligned} \]继续分块。
#include <algorithm>
#include <cstring>
#include <cmath>
#include <vector>
#include <cstdio>
#include <map>
#include <iostream>
using namespace std;
typedef long long ll;
const ll N = 1e7 + 11;
const ll p = 20101009;
vector<ll> prime;
ll mu[N], vis[N], S[N], SS[N];
ll H(ll, ll);
void Table(ll);
int main(){
ll n, m, Ans = 0;
scanf("%lld %lld", &n, &m);
if(n > m) swap(n, m);
Table(max(n,m));
for(ll l = 1, r; l <= n; l = r + 1){
r = min(n / (n/l), m / (m/l));
Ans += H(n/l, m/l) * ((SS[r] - SS[l - 1] + p) % p) % p;
Ans %= p;
}
printf("%lld\n", Ans);
return 0;
}
ll H(ll n, ll m){
if(n > m) swap(n, m);
ll Ans = 0;
for(ll l = 1, r; l <= n; l = r + 1){
r = min(n / (n/l), m / (m/l));
ll x = n/l, y = m/l;
ll rhs = ((x + 1) * x / 2 % p) * ((y + 1) * y / 2 % p) % p;
Ans += rhs * ((S[r] - S[l - 1] + p) % p) % p;
Ans %= p;
}
return Ans;
}
void Table(ll num){
mu[1] = 1;
for(ll i = 2; i <= num; i ++){
if(!vis[i]){
vis[i] = i;
mu[i] = -1;
prime.push_back(i);
}
for(auto j:prime){
if(i * j > num) break;
vis[i * j] = j;
if(i % j == 0){
mu[i * j] = 0;
break;
}else{
mu[i * j] = -mu[i];
}
}
}
for(ll i = 1; i <= num; i ++)
S[i] = (S[i - 1] + i*i%p*mu[i]%p) % p, SS[i] = (SS[i - 1] + i) % p;
}