$ Ans = sum_i^Nsum_j^Mlcm(i,j)$
因为 $lcm(x,y)=xy/gcd(x,y)$
所以 $Ans = sum_i^Nsum_j^Mij/gcd(i,j)=sum_i^Nsum_j^M( i/gcd(i,j) )( j/gcd(i,j) ) cdot gcd(i,j)$
改成枚举所有可能的 $gcd(i,j)=t$，变成 $Ans=sum_tsum_i^{left lfloor frac{N}{t} ight floor}sum_j^{left lfloor frac{M}{t} ight floor}ijt[gcd(i,j)==1]$
显然莫比乌斯反演
设 $f[x]=sum_tsum_i^{left lfloor frac{N}{t} ight floor}sum_j^{left lfloor frac{M}{t} ight floor}ijt[gcd(i,j)==x]$，$F[x]=sum_tsum_i^{left lfloor frac{N}{t} ight floor}sum_j^{left lfloor frac{M}{t} ight floor}ijt[x|gcd(i,j)]$
那么 $Ans=f[1]$
那么有 $F[x]=sum_{x|d}f[d]$，直接反演
$f[1]=sum_{1|d}mu(d/1)sum_tsum_i^{left lfloor frac{N}{t} ight floor}sum_j^{left lfloor frac{M}{t} ight floor}ijt[d|gcd(i,j)]$
$=sum_dmu(d)sum_tsum_i^{left lfloor frac{N}{td} ight floor}sum_j^{left lfloor frac{M}{td} ight floor}ijtd^2$
枚举 $T=td$，则原式 $=sum_dmu(d)sum_{d|T}sum_i^{left lfloor frac{N}{T} ight floor}sum_j^{left lfloor frac{M}{T} ight floor}ijTd$
$=(sum_dmu(d)d)(sum_{d|T}T)(sum_i^{left lfloor frac{N}{T} ight floor}i)(sum_j^{left lfloor frac{M}{T} ight floor}j)$
改成先枚举 $T$ 再枚举 $d$
$=(sum_TT)(sum_{d|T}mu(d)d)(sum_i^{left lfloor frac{N}{T} ight floor}i)(sum_j^{left lfloor frac{M}{T} ight floor}j)$
显然后面的 $(sum_i^{left lfloor frac{N}{T} ight floor}i)$ 和 $(sum_j^{left lfloor frac{M}{T} ight floor}j)$ 可以预处理
考虑能不能预处理 $g[T]=sum_{d|T}mu(d)d$
发现好像可以线性筛，考虑假设此时已经知道了 $g[x]$，且 $y=xp$，$p$ 为 $y$ 唯一一个最小的质因数
那么 $g[xp]=sum_{d|xp}mu(d)d=(sum_{d|x}mu(d)d)+(sum_{d|x}mu(dp)dp)$
因为 $p$ 为 $y$ 唯一一个最小的质因数，所以 $gcd(d,p)=1$ ，所以 $mu(dp)=-mu(d)$
则 $g[xp]=(sum_{d|x}mu(d)d)-(sum_{d|x}mu(d)dp)=(sum_{d|x}mu(d)d)(1-p)=(1-p)g[x]$
如果 $y=xp$，$p$ 为 $y$ 最小的一个质因数但不唯一，则根据 $g$ 的定义发现 $p$ 对 $x$ 没有贡献，所以此时 $g[xp]=g[x]$
然后直接线性筛就行了，具体看代码
#include<iostream> #include<cstdio> #include<algorithm> #include<cmath> #include<cstring> using namespace std; typedef long long ll; inline int read() { int x=0,f=1; char ch=getchar(); while(ch<'0'||ch>'9') { if(ch=='-') f=-1; ch=getchar(); } while(ch>='0'&&ch<='9') { x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); } return x*f; } const int N=2e7+7,mo=20101009; inline int fk(int x) { return x>=mo ? x-mo : x; } int n,m,mi,pri[N],mu[N],g[N],tot; bool not_pri[N]; ll ans; void pre() { mu[1]=1; g[1]=1; for(int i=2;i<=mi;i++) { if(!not_pri[i]) pri[++tot]=i,mu[i]=mo-1,g[i]=mo+1-i; for(int j=1;j<=tot;j++) { ll t=1ll*i*pri[j]; if(t>mi) break; not_pri[t]=1; if(!(i%pri[j])) { g[t]=g[i]; break; } mu[t]=mo-mu[i]; g[t]= 1ll*g[i]*(mo+1-pri[j])%mo; } } } inline int sum(int x) { return 1ll*x*(x+1)/2%mo; } int main() { n=read(),m=read(); mi=min(n,m); pre(); for(int _d=1;_d<=mi;_d++) { ll t=g[_d],s=1ll*sum(n/_d)*sum(m/_d)%mo; ans=fk(ans+ ((t*s%mo)*_d%mo) ); } printf("%lld ",(ans%mo+mo)%mo); return 0; }