D. [国家集训队] Crash的数字表格 / JZPTAB
\[\begin{aligned}
&\sum_{i=1}^{n}\sum_{j=1}^{m}\operatorname{lcm}(i,j)\\
=&\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{i\times j}{\gcd(i,j)}\\
=&\sum_{d=1}^{\min(n,m)}\frac{1}{d}\sum_{i=1}^{n}\sum_{j=1}^{m}i\times j\times[\gcd(i,j)=d]\\
\end{aligned}
\]
不妨令 \(n\le m\)，即
\[\sum_{d=1}^{n}\frac{1}{d}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}(i\times d)\times(j\times d)\times[\gcd(i,j)=1]
\]
莫比乌斯反演，得
\[\begin{aligned}
&\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}i\times j\sum_{k\mid\gcd(i,j)}\mu(k)\\
=&\sum_{d=1}^{n}d\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum_{i=1}^{n}i\times[k\mid i]\sum_{j=1}^{m}j\times[k\mid j]\\
=&\sum_{d=1}^{n}d\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\times(k\times\frac{\lfloor\frac{n}{d\times k}\rfloor\times(\lfloor\frac{n}{d\times k}\rfloor+1)}{2})\times(k\times\frac{\lfloor\frac{m}{d\times k}\rfloor\times(\lfloor\frac{m}{d\times k}\rfloor+1)}{2})
\end{aligned}
\]
于是可以调和级数求。预处理 \(i^2\) 和 \(\frac{i\times(i+1)}{2}\) 减少取模即可。
Code
#include<bits/stdc++.h>
#define ll long long
#define il inline
using namespace std;
namespace asbt{
const int maxn=1e7+5,mod=20101009;
int n,m,prn,prm[maxn],mu[maxn],f[maxn],g[maxn];
bool npr[maxn];
il void euler(int n=1e7){mu[1]=npr[1]=1;for(int i=2;i<=n;i++){if(!npr[i]){prm[++prn]=i,mu[i]=-1;}for(int j=1;j<=prn&&i*1ll*prm[j]<=n;j++){npr[i*prm[j]]=1;if(i%prm[j]==0){mu[i*prm[j]]=0;break;}mu[i*prm[j]]=-mu[i];}}
}
il void init(int n=1e7){for(int i=1;i<=n;i++){f[i]=i*1ll*i%mod;g[i]=i*1ll*(i+1)/2%mod;}
}
int main(){ios::sync_with_stdio(0),cin.tie(0);euler(),init();cin>>n>>m;if(n>m){swap(n,m);}int ans=0;for(int i=1;i<=n;i++){int res=0;for(int j=1;j<=n/i;j++){int tmp=g[n/i/j]*1ll*g[m/i/j]%mod*f[j]%mod;if(mu[j]==1){(res+=tmp)%=mod;}else if(mu[j]==-1){(res+=mod-tmp)%=mod;}}ans=(res*1ll*i+ans)%mod;}cout<<ans;return 0;
}
}
int main(){return asbt::main();}
然而这还不够优秀。观察这个式子：
\[\begin{aligned}
&\sum_{d=1}^{n}d\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\times(k\times\frac{\lfloor\frac{n}{d\times k}\rfloor\times(\lfloor\frac{n}{d\times k}\rfloor+1)}{2})\times(k\times\frac{\lfloor\frac{m}{d\times k}\rfloor\times(\lfloor\frac{m}{d\times k}\rfloor+1)}{2})\\
=&\sum_{d=1}^{n}d\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\times k^2\times\frac{\lfloor\frac{n}{d\times k}\rfloor\times(\lfloor\frac{n}{d\times k}\rfloor+1)}{2}\times\frac{\lfloor\frac{m}{d\times k}\rfloor\times(\lfloor\frac{m}{d\times k}\rfloor+1)}{2}
\end{aligned}
\]
考虑设 \(\operatorname{sum}(n,m)=\sum_{k=1}^{n}\mu(k)\times k^2\times\frac{\lfloor\frac{n}{k}\rfloor\times(\lfloor\frac{n}{k}\rfloor+1)}{2}\times\frac{\lfloor\frac{m}{k}\rfloor\times(\lfloor\frac{m}{k}\rfloor+1)}{2}\)，给 \(\mu(k)\times k^2\) 做前缀和，可以数论分块做。于是原式变成：
\[\sum_{d=1}^{n}d\times\operatorname{sum}(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)
\]
这个东西还可以数论分块，于是时间复杂度线性。
Code
#include<bits/stdc++.h>
#define ll long long
#define il inline
using namespace std;
namespace asbt{
const int maxn=1e7+5,mod=20101009;
int n,m,prn,prm[maxn],mu[maxn],f[maxn],g[maxn],h[maxn];
bool npr[maxn];
il void euler(int n=1e7){mu[1]=npr[1]=1;for(int i=2;i<=n;i++){if(!npr[i]){prm[++prn]=i,mu[i]=-1;}for(int j=1;j<=prn&&i*1ll*prm[j]<=n;j++){npr[i*prm[j]]=1;if(i%prm[j]==0){mu[i*prm[j]]=0;break;}mu[i*prm[j]]=-mu[i];}}
}
il void init(int n=1e7){for(int i=1;i<=n;i++){f[i]=(f[i-1]+i)%mod;g[i]=((g[i-1]+mu[i]*1ll*i*i)%mod+mod)%mod;h[i]=i*1ll*(i+1)/2%mod;}
}
il int calc(int n,int m){int res=0;for(int i=1;i<=n;i++){int j=min(n/(n/i),m/(m/i));res=((g[j]-g[i-1]+mod)*1ll*h[n/i]%mod*h[m/i]+res)%mod;i=j;}return res;
}
int main(){ios::sync_with_stdio(0),cin.tie(0);euler(),init();cin>>n>>m;if(n>m){swap(n,m);}int ans=0;for(int i=1;i<=n;i++){int j=min(n/(n/i),m/(m/i));ans=((f[j]-f[i-1]+mod)*1ll*calc(n/i,m/i)+ans)%mod;i=j;}cout<<ans;return 0;
}
}
int main(){return asbt::main();}