(color{#0066ff}{ 题目描述 })
今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时整除a和b的最小正整数。例如，LCM(6, 8) = 24。
回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张NM的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个45的表格如下：
1 2 3 4 5
2 2 6 4 10
3 6 3 12 15
4 4 12 4 20
看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod20101009的值。
(color{#0066ff}{输入格式})
输入的第一行包含两个正整数，分别表示N和M。
(color{#0066ff}{输出格式})
输出一个正整数，表示表格中所有数的和mod20101009的值。
(color{#0066ff}{输入样例})
4 5
(color{#0066ff}{输出样例})
122
(color{#0066ff}{数据范围与提示})
30%的数据满足(N, M≤ 10^3。)
70%的数据满足(N, M≤ 10^5。)
100%的数据满足(N, M≤ 10^7。)
(color{#0066ff}{ 题解 })
题目就是要求
[sum_{i=1}^n sum_{j=1}^m lcm(i,j)
]
显然要转为gcd的东西
[sum_{i=1}^n sum_{j=1}^m frac {i*j} {gcd(i,j)}
]
gcd套路，考虑枚举gcd
[sum_{d=1}^{min(n,m)} sum_{i=1}^n sum_{j=1}^m [gcd(i,j)==d] frac {i*j} {d}
]
老套路，把d除上去，注意后面的(i*j)，各除去了一个d，所以要乘回来(d^2)
[sum_{d=1}^{min(n,m)} sum_{i=1}^{lfloorfrac{n}{d}
floor} sum_{j=1}^{lfloorfrac{m}{d}
floor} [gcd(i,j)==1] i*j*d
]
然后发现d可以提出来
[sum_{d=1}^{min(n,m)} dsum_{i=1}^{lfloorfrac{n}{d}
floor} sum_{j=1}^{lfloorfrac{m}{d}
floor} [gcd(i,j)==1] i*j
]
我们令(f(x, y))等于后面那一坨，即
[f(a,b)=sum_{i=1}^{a} sum_{j=1}^{b} [gcd(i,j)==1] i*j
]
看到这种形式，继续利用(mu * i =e)，来化简(f)
[sum_{i=1}^{a} sum_{j=1}^{b} sum_{d|gcd(i,j)} mu(d)* i*j
]
然后枚举gcd
[sum_{d=1}^{min(a,b)} mu(d) sum_{d|i}^{a} sum_{d|j}^b i*j
]
后面的显然不好枚举，考虑(i=p_1d,j=p_2d)，枚举(p_1,p_2)
[sum_{d=1}^{min(a,b)} mu(d) sum_{i=1}^{lfloorfrac a k
floor} sum_{j=1}^{lfloorfrac b k
floor} i*d*j*d
]
然后我们又可以开心的把d提出来了
[sum_{d=1}^{min(a,b)} mu(d) * d^2 sum_{i=1}^{lfloorfrac a k
floor} sum_{j=1}^{lfloorfrac b k
floor} i*j
]
后面的东西好像很好求呢
[令g(a,b)=sum_{i=1}^asum_{j=1}^bi*j=frac{a*(a+1)}{2}*frac{b*(b+1)}{2}
]
嗯，(O(1))就出来了
至此，本题就完了
总结一下
[sum_{i=1}^n sum_{j=1}^m lcm(i,j)=sum_{d=1}^{min(n,m)} d *f(lfloorfrac{n}{d}
floor,lfloorfrac{m}{d}
floor)
]
[f(a,b)=sum_{d=1}^{min(a,b)} mu(d) * d^2* g(lfloorfrac a k
floor,lfloorfrac b k
floor)
]
[g(a,b)=frac{a*(a+1)}{2}*frac{b*(b+1)}{2}
]
最后一个(O(1))求，前两个数列分块即可
#include<bits/stdc++.h>
#define LL long long
LL in() {
char ch; LL x = 0, f = 1;
while(!isdigit(ch = getchar()))(ch == '-') && (f = -f);
for(x = ch ^ 48; isdigit(ch = getchar()); x = (x << 1) + (x << 3) + (ch ^ 48));
return x * f;
}
const int mod = 20101009;
const int maxn = 1e7 + 10;
LL pri[maxn], tot, mu[maxn], s[maxn];
bool vis[maxn];
void predoit() {
mu[1] = 1;
for(int i = 2; i < maxn; i++) {
if(!vis[i]) pri[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && (LL)i * pri[j] < maxn; j++) {
vis[i * pri[j]] = true;
if(i % pri[j] == 0) break;
else mu[i * pri[j]] = -mu[i];
}
}
for(int i = 1; i < maxn; i++) s[i] = (s[i - 1] + 1LL * mu[i] * i * i) % mod;
}
LL getsum1(LL l, LL r) {
return ((l + r) * (r - l + 1) >> 1LL) % mod;
}
LL getsum2(LL a, LL b) {
return (((a * (a + 1) >> 1LL) % mod) * ((b * (b + 1) >> 1LL) % mod)) % mod;
}
LL getans(LL a, LL b) {
LL ans = 0;
for(LL l = 1, r; l <= std::min(a, b); l = r + 1) {
r = std::min(a / (a / l), b / (b / l));
(ans += ((((s[r] - s[l - 1]) % mod + mod) % mod) * getsum2(a / l, b / l) % mod)) %= mod;
}
return ans;
}
LL work(LL a, LL b) {
LL ans = 0;
for(LL l = 1, r; l <= std::min(a, b); l = r + 1) {
r = std::min(a / (a / l), b / (b / l));
(ans += getsum1(l, r) * getans(a / l, b / l) % mod) %= mod;
}
return ans;
}
int main() {
predoit();
int n = in(), m = in();
printf("%lld
", work(n, m));
return 0;
}