题目
P1829 [国家集训队]Crash的数字表格 / JZPTAB
求：
\[\sum_{i=1}^n\sum_{j=1}^mlcm(i,j)
\]
分析
先将式子化简：
\[\large
\begin{split}
\sum_{i=1}^n\sum_{j=1}^mlcm(i,j)
&=\sum_{i=1}^n\sum_{j=1}^{m}\sum_{d\vert i,d\vert j,\gcd(\frac id,\frac jd)=1}\frac {ij}{d}\\
&=\sum_{d=1}^nd\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}ij[\gcd(i,j)=1]\\
\end{split}
\]
这时我们单独计算后半部分，这里没有使用莫比乌斯反演而是直接用 \(\large \mu*I=\varepsilon\) 的性质：
\[\large
\begin{split}
S(n,m)
&=\sum_{i=1}^{n}\sum_{j=1}^{m}ij[\gcd(i,j)=1]\\
&=\sum_{i=1}^{n}\sum_{j=1}^{m}ij\sum_{d\vert \gcd(i,j)}\mu(d)\\
&=\sum_{d=1}^{n}\mu(d)\cdot d^2\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}ij\\
&=\sum_{d=1}^{n}\mu(d)\cdot d^2\sum_{i=1}^{\lfloor\frac nd\rfloor}i\sum_{j=1}^{\lfloor\frac md\rfloor}j\\
&=\sum_{d=1}^{n}\mu(d)\cdot d^2\cdot G(\lfloor\frac nd\rfloor)G(\lfloor\frac md\rfloor)
\end{split}
\]
其中，\(\large G(x)=\frac {x(x+1)}{2}\) ，表示等差数列求和。
很容易发现，对于 \(\large S\) 函数，可以通过整除分块来计算，按照 \(\large\lfloor\frac nd\rfloor\) 以及 \(\large \lfloor\frac md\rfloor\) 来分块即可，前面的部分就是预处理一下前缀和即可。
现在代回原式子：
\[\large \sum_{d=1}^nd\cdot S({\lfloor\frac nd\rfloor},{\lfloor\frac md\rfloor})
\]
然后我们可以发现，原式还需要一次除法分块，而第二次除法分块将借助第一次的作为上界进行分块，这就是二次除法分块，这样做的复杂度可以证明是 \(\large O(n^{\frac 34})\) 的。
扩展：貌似有单次除法分块的做法。
代码
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define rep(i,x,y) for(int i=(x);i<=(y);i++)
#define dep(i,y,x) for(int i=(y);i>=(x);i--)
const int N=1e7+5,M=2e5+5,MOD=20101009;
int n,m,k;
int prime[N],mu[N],cnt,pre[N];
ll val[N],Ans;
bool isprime[N];
inline int inc(int x,int y){x+=y;return x>=MOD?x-MOD:x;}
inline void incc(int &x,int y){x+=y;if(x>=MOD) x-=MOD;}
void GetPrimes(int v){
mu[1]=1;
for(int i=2;i<=v;i++){
if(!isprime[i]) prime[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&prime[j]*i<=v;j++){
isprime[i*prime[j]]=true;
if(i%prime[j]==0) break;
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=v;i++) pre[i]=(pre[i-1]+1ll*i*i%MOD*(mu[i]+MOD))%MOD;
return ;
}
inline int GetSum(int x){return (1ll*x*(x+1)/2)%MOD;}
inline int Sum(int x,int y){
int res=0;
if(x>y) swap(x,y);
for(int l=1,r;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
incc(res,1ll*(inc((pre[r]-pre[l-1])%MOD,MOD))*GetSum(x/l)%MOD*GetSum(y/l)%MOD);
}
return res;
}
inline int Calc(int x,int y){
int res=0;
for(int l=1,r;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
incc(res,1ll*inc((GetSum(r)-GetSum(l-1))%MOD,MOD)%MOD*inc(Sum(x/l,y/l)%MOD,MOD)%MOD);
}
return res;
}
signed main(){
ios::sync_with_stdio(false);
GetPrimes(1e7);
int t;cin>>t;
while(t--){
cin>>n>>m;
if(n>m) swap(n,m);
cout<<Calc(n,m)<<endl;
}
return 0;
}