题目要求，求：
[sum_{i=1}^nsum_{j=1}^mlcm(i,j)
]
先转化为gcd的形式，然后枚举gcd。
[sum_{i=1}^nsum_{j=1}^msum_{d=1}^nfrac{ij}{d}[gcd(i,j) = d]
]
把d除进去，套用莫比乌斯函数的性质：
[sum_{d=1}^ndsum_{i=1}^{frac{n}{d}}isum_{j=1}^{frac{m}{d}}jsum_{p|i.p|j}mu(p)
]
继续替换得到：
[sum_{d=1}^ndsum_{p=1}^{frac{n}{d}}p^2mu(p)s(frac{n}{dp})s(frac{m}{dp})
]
其中s(n)表示(sum_{i=1}^ni)
这个其实已经可以做了，直接枚举d，然后里面使用整除分块完成。这个看起来复杂度是(O(nsqrt{n}))的，但是实际上它每次没有跑满，复杂度是(O(n))左右的。
不过这个是弱化版，加强版还要求处理多组数据，这时候上面的做法就不好使了。继续推导，设(T=dp)
[sum_{T=1}^nTsum_{d|T}dmu(d)s(frac{n}{T})s({frac{m}{T}})
]
设(h(T) = sum_{d|T}dmu(d))问题在于怎么能快速求出(h(T))
这并不是一个积性函数，但是我们仍然能线性把它筛出来。首先考虑T为质数的时候，这时候显然(h(T) = 1 - T)。如果现在加入一个已经出现在T中的质因子p，那么所有T原来的因子在乘上这个p之后，p的指数必然大于1，也就是说其莫比乌斯函数的值是0，原来的因子不变，所以(h(Tp) = h(T)).再考虑新加入一个因子的情况。加入之后，原来所有的因子其莫比乌斯函数的值变成其相反数，而且因为前面还有乘p，所以(h(Tp) = (1-p)h(T)),即(h(Tp) = h(p)h(T))
所以我们把它线性筛出来，之后整除分块做即可。单次复杂度(O(sqrt{n}))。注意因为此题有取模，所以要注意前缀和相减的时候，先变成正数。
弱化版代码：
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#include<set>
#include<vector>
#include<map>
#include<queue>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('
')
#define fr friend inline
#define y1 poj
#define mp make_pair
#define pr pair<int,int>
#define fi first
#define sc second
#define pb push_back
using namespace std;
typedef long long ll;
const int M = 10000005;
const int INF = 1000000009;
const ll mod = 20101009;
const double eps = 1e-7;
int read()
{
int ans = 0,op = 1;char ch = getchar();
while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
return ans * op;
}
ll k,p[M],mu[M],tot,n,m;
ll sum[M],ans,pre[M];
bool np[M];
void euler()
{
np[1] = 1,mu[1] = 1;
rep(i,2,M-2)
{
if(!np[i]) p[++tot] = i,mu[i] = -1;
for(int j = 1;i * p[j] <= M-2;j++)
{
np[i * p[j]] = 1;
if(!(i % p[j])) break;
mu[i * p[j]] = -mu[i];
}
}
rep(i,1,M-2) pre[i] = (pre[i-1] + (ll)i * (ll)i % mod * mu[i]) % mod;
rep(i,1,M-2) sum[i] = (sum[i-1] + (ll)i) % mod;
}
int main()
{
euler();
n = read(),m = read();
int lim = min(n,m);
rep(d,1,lim)
{
int a = n / d,b = m / d,c = min(a,b);
ll cur = 0;
for(int i = 1,j;i <= c;i = j + 1)
{
j = min(a / (a / i),b / (b / i));
cur += ((pre[j] - pre[i-1] + mod) % mod * sum[a / i] % mod * sum[b / i] % mod),cur %= mod;
}
ans += (cur * (ll)d % mod),ans %= mod;
}
printf("%lld
",ans);
return 0;
}
强化版代码：
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#include<set>
#include<vector>
#include<map>
#include<queue>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('
')
#define fr friend inline
#define y1 poj
#define mp make_pair
#define pr pair<int,int>
#define fi first
#define sc second
#define pb push_back
using namespace std;
typedef long long ll;
const int M = 10000005;
const int INF = 1000000009;
const ll mod = 20101009;
const double eps = 1e-7;
int read()
{
int ans = 0,op = 1;char ch = getchar();
while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
return ans * op;
}
ll k,p[M],mu[M],tot,n,m;
ll sum[M],ans,pre[M],f[M];
bool np[M];
void euler()
{
np[1] = 1,f[1] = 1;
rep(i,2,M-2)
{
if(!np[i]) p[++tot] = i,f[i] = (1 - i + mod) % mod;
for(int j = 1;i * p[j] <= M-2 && j <= tot;j++)
{
np[i * p[j]] = 1;
if(!(i % p[j])) {f[i * p[j]] = f[i];break;}
f[i * p[j]] = f[i] * f[p[j]] % mod;
}
}
rep(i,1,M-2) pre[i] = (pre[i-1] + ((ll)i * f[i] % mod)) % mod;
rep(i,1,M-2) sum[i] = (sum[i-1] + (ll)i) % mod;
}
int main()
{
euler();
n = read(),m = read();
int lim = min(n,m);
for(int i = 1,j;i <= lim;i = j + 1)
{
j = min(n / (n / i),m / (m / i));
ans += ((pre[j] - pre[i-1] + mod) % mod * sum[n / i] % mod * sum[m / i] % mod),ans %= mod;
}
printf("%lld
",ans);
return 0;
}