我膜拜PoPoQQQ
http://www.lydsy.com/JudgeOnline/problem.php?id=2154
∑ni=1∑mj=1lcm(i,j)=∑ni=1∑mj=1i×jgcd(i,j)
令
gcd(i,j)=d
那么
∑ni=1∑mj=1i×jgcd(i,j)=∑min(n,m)d=1d2×g(n,m,d)d
这里令
g(n,m,d)
表示为在
g(n,m,d)=∑⌊nd⌋i=1∑⌊md⌋j=1[gcd(i,j)==1]i×j=∑⌊nd⌋i=1∑⌊md⌋j=1∑t|i,t|jμ(t)i×j
那么就有
∑min(n,m)d=1d2×g(n,m,d)d=∑min(n,m)d=1d×g(n,m,d)
那么令
Sum(a,b)=∑ai=1∑bj=1i×j=(1+2....+a)×(1+2.....+b)=a(a+1)2×b(b+1)2
因为
g(n,m,d)=∑⌊nd⌋i=1∑⌊md⌋j=1∑t|i,t|jμ(t)i×j=∑min(n,m)dt=1μ(t)×t2×Sum(⌊ntd⌋,⌊mtd⌋)
#include <cstdio>
#include <iostream>
#include <cstring>
#include <climits>
#include <algorithm>
using namespace std;
#define mod 20101009LL
const int MAXN = 10000000;
bool notprime[MAXN+10];
int prime[MAXN+10], mu[MAXN+10], sum[MAXN+10];
long long Max;
void Init(){
mu[1] = 1;
int tmp;
for(int i=2;i<=Max;i++){
if(!notprime[i]){
mu[i] = -1;
prime[++prime[0]] = i;
}
for(int j=1;j<=prime[0]&&(tmp=prime[j]*i)<=Max;j++){
notprime[tmp] = true;
if(i%prime[j] == 0){
mu[tmp] = 0;
break;
}
mu[tmp] = -mu[i];
}
}
for(long long i=1;i<=Max;i++)
sum[i]=(sum[i-1]+(i*i*mu[i])%mod)%mod;
}
long long Sum(long long n, long long m){
return (n*(n+1)/2%mod)*(m*(m+1)/2%mod)%mod;
}
long long F(long long n, long long m){
long long ret = 0, last;
for(long long i=1;i<=n;i=last+1){
last = min(n/(n/i), m/(m/i));
ret += (sum[last] - sum[i-1]) * Sum(n/i, m/i)%mod;
ret %= mod;
}
return ret;
}
long long solve(int n, int m){
long long ret = 0, last;
for(long long i=1;i<=n;i=last+1){
last = min(n/(n/i), m/(m/i));
ret += (i+last)*(last-i+1)/2 * F(int(n/i), int(m/i))%mod;
ret %= mod;
}
return ret;
}
int main(){
long long a, b;
cin>>a>>b;
if(a > b) swap(a, b);
Max = a;
Init();
cout<<(solve(a, b)%mod+mod)%mod<<endl;
return 0;
}