简要题意:
求
[sum_{i=1}^n sum_{j=1}^m operatorname{lcm}(i,j)
]
一言不合就推式子。
[sum_{i=1}^n sum_{j=1}^m operatorname{lcm}(i,j)
]
[= sum_{i=1}^n sum_{j=1}^m frac{ij}{gcd(i,j)}
]
[= sum_{d=1}^{min(n,m)} frac{1}{d} sum_{i=1}^m sum_{j=1}^n ij [gcd(i,j)==d]
]
[= sum_{d=1}^{min(n,m)} frac{1}{d} sum_{i=1}^{lfloor frac{n}{d}
floor} sum_{j=1}^{lfloor frac{m}{d}
floor} ijd^2 [gcd(i,j)==1]
]
[= sum_{d=1}^{min(n,m)} d sum_{i=1}^{lfloor frac{n}{d}
floor} sum_{j=1}^{lfloor frac{m}{d}
floor} [gcd(i,j)==1]
]
[= sum_{d=1}^{min(n,m)} d sum_{i=1}^{lfloor frac{n}{d}
floor} sum_{j=1}^{lfloor frac{m}{d}
floor} ij sum_{x | gcd(i,j)} mu_x
]
[= sum_{d=1}^{min(n,m)} d sum_{x=1}^{min(lfloor frac{n}{d}
floor,lfloor frac{m}{d}
floor)} x^2 mu_x s_{lfloor frac{n}{dx}
floor} s_{lfloor frac{m}{dx}
floor}
]
其中
[s_x = sum_{i=1}^x i = frac{x imes (x+1)}{2}
]
回到原式:
[= sum_{T=1}^n s_{lfloor frac{n}{T}
floor} s_{lfloor frac{m}{T}
floor} ig(T sum_{d|T} d mu_T ig)
]
括号内的东西我们预处理,括号外面的直接 整除分块。(实际上没有这个必要了)
那么如何预处理呢?单独搞出这个式子。
[T sum_{d|T} d mu_T
]
[T mu_T sum_{d|T} d
]
咦?是不是很熟悉? (T mu_T) 我们直接预处理就可以了。考虑,令:
[f_i = sum_{d|i} d
]
(其实 (f_i) 就是 (i) 的因数之和)
那么可得:
[f_{i imes j} = f_i imes f_j ig([gcd(i,j)]==1 ig)
]
即 (f) 为 积性函数。
所以用线性筛预处理即可。
时间复杂度:(O(n)).
实际得分:(100pts).
#pragma GCC optimize(2)
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e7+1;
const ll MOD=20101009;
inline int read(){char ch=getchar(); int f=1;while(ch<'0' || ch>'9') {if(ch=='-') f=-f; ch=getchar();}
int x=0;while(ch>='0' && ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();return x*f;}
int prime[N],mu[N],n,m;
int cnt=0; bool h[N];
ll ans=0;
inline void Euler(int n) {
mu[1]=1; for(register int i=2;i<=n;i++) {
if(!h[i]) mu[i]=MOD-i+1,prime[++cnt]=i;
for(register int j=1;j<=cnt && i*prime[j]<=n;j++) {
h[i*prime[j]]=1;
if(i%prime[j]==0) {mu[i*prime[j]]=mu[i];break;}
mu[i*prime[j]]=(1ll*mu[i]*mu[prime[j]])%MOD;
}
} for(register int i=1;i<=n;i++) mu[i]=1ll*mu[i]*i%MOD;
for(register int i=1;i<=n;i++) mu[i]=(mu[i]+mu[i-1])%MOD;
} //欧拉筛预处理
inline ll Sieve(ll n) {return (n*(n+1)>>1)%MOD;} //求 s
int main() {
Euler(N-1);
n=read(),m=read();
if(n>m) swap(n,m);
for(int l=1;l<=n;) {
int r=min(n/(n/l),m/(m/l));
ans=(ans+1ll*(mu[r]-mu[l-1]+MOD)*Sieve(n/l)%MOD*Sieve(m/l)%MOD)%MOD;
l=r+1; //整除分块
} printf("%lld
",ans);
return 0;
}