题目大意：
给定$n,m(n,m\leq10^7)$，求$\displaystyle\sum_{x=1}^n\sum_{y=1}^m{\rm lcm}(x,y)$。思路：
令$d=\gcd(x,y),x=ad,y=bd$，则：$$\begin{align*} 原式&=\sum_{d=1}^{\min(n,m)}d\sum_{a=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{d}\rfloor}ab[\gcd(ab)=1]\\ &=\sum_{d=1}^{\min(n,m)}d\sum_{a=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{d}\rfloor}ab\sum_{p|\gcd(a,b)}\mu(p)\end{align*}$$ 令$p=\gcd(a,b),a=jp,b=kp$，则：$$\begin{align*} 原式&=\sum_{d=1}^{\min(n,m)}d\sum_{p=1}^{\min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\mu(p)p^2\sum_{j=1}^{\lfloor\frac{n}{dp}\rfloor}\sum_{k=1}^{\lfloor\frac{m}{dp}\rfloor}jk\\ &=\sum_{d=1}^{\min(n,m)}d\sum_{p=1}^{\min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\mu(p)p^2\frac{(\lfloor\frac{n}{dp}\rfloor+1)\lfloor\frac{n}{dp}\rfloor}{2}\cdot\frac{(\lfloor\frac{m}{dp}\rfloor+1)\lfloor\frac{m}{dp}\rfloor}{2}\\\end{align*}$$ 预处理$\mu(p)p^2$的前缀和，暴力枚举$d$。对于$\lfloor\frac{n}{dp}\rfloor$和$\lfloor\frac{m}{dp}\rfloor$分别相等的$p$可以分块计算。1 #include2 #include 3 #include 4 typedef long long int64; 5 inline int getint() { 6 register char ch; 7 while(!isdigit(ch=getchar())); 8 register int x=ch^'0'; 9 while(isdigit(ch=getchar())) x=(((x<<2)+x)<<1)+(ch^'0');10 return x;11 }12 const int N=10000001,M=664580,mod=20101009;13 bool vis[N];14 int mu[N],p[M],sum[N];15 inline void sieve(const int lim) {16 mu[1]=1;17 for(register int i=2;i<=lim;i++) {18 if(!vis[i]) {19 p[++p[0]]=i;20 mu[i]=-1;21 }22 for(register int j=1;j<=p[0]&&i*p[j]<=lim;j++) {23 vis[i*p[j]]=true;24 if(i%p[j]==0) {25 mu[i*p[j]]=0;26 break;27 }28 mu[i*p[j]]=-mu[i];29 }30 }31 for(register int i=1;i<=lim;i++) {32 sum[i]=(sum[i-1]+(int64)mu[i]*i%mod*i%mod)%mod;33 }34 }35 int main() {36 const int n=getint(),m=getint(),lim=std::min(n,m);37 sieve(lim);38 int ans=0;39 for(register int d=1;d<=lim;d++) {40 int tmp=0;41 const int x=n/d,y=m/d,lim=std::min(x,y);42 for(register int i=1,j;i<=lim;i=j+1) {43 j=std::min(x/(x/i),y/(y/i));44 (tmp+=(sum[j]-sum[i-1]+mod)%mod*((int64)(x/i+1)*(x/i)/2%mod)%mod*((int64)(y/i+1)*(y/i)/2%mod)%mod)%=mod;45 }46 (ans+=((int64)tmp*d)%mod)%=mod;47 }48 printf("%d\n",ans);49 return 0;50 }