【BZOJ2154】—Crash的数字表格（莫比乌斯反演+整除分块）

传送门

题意：求$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j),n,m\le1e7$

---

$Solution$

考虑到$lcm$无法处理，我们先变成$gcd$的形式

$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}\frac {i*j}{gcd(i,j)}$

考虑枚举$gcd$
则

$ans=\sum_{d=1}^{min(n,m)}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d]\frac {ij}d$$=\sum_{d=1}^{min(n,m)}\sum_{d|i}^{n}\sum_{d|j}^{m}[gcd(i,j)=d]\frac {ij}d$

令$i'=\frac d i,j'=\frac d j$(即从枚举$d$的倍数变成枚举是$d$的几倍)
则(为了方便仍用$i,j$表示$i',j'$)

$ans=\sum_{d=1}^{min(n,m)}\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac m d}[gcd(di,dj)=d] \frac{id*jd}{d}$$=\sum_{d=1}^{min(n,m)}d\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac m d}[gcd(i,j)=1]ij$

考虑有个莫比乌斯函数的简单结论$\sum_{d|n}\mu(d)=[n==1]$

考虑的式子中有个$[gcd(i,j)=1]$
则：

$ans=\sum_{d=1}^{min(n,m)}d\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac m d}\sum_{T|gcd(i,j)}\mu(T)ij$

考虑把$T$提到前面来：

$ans=\sum_{d=1}^{min(n,m)}d\sum_{T=1}^{min(\frac n d,\frac m d)}\mu(T)\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac m d}ij[T|gcd(i,j)]$

考虑令$i'=i*T,j'=j*T$
则(同样为了方便直接用$i,j$表示了)

$ans=\sum_{d=1}^{min(n,m)}d\sum_{T=1}^{min(\frac n d,\frac m d)}\mu(T)\sum_{i=1}^{\frac {n}{dT}}\sum_{j=1}^{\frac{m}{dT}}ijT^2$
$=\sum_{d=1}^{min(n,m)}d\sum_{T=1}^{min(\frac n d,\frac m d)}\mu(T)T^2\sum_{i=1}^{\frac {n}{dT}}\sum_{j=1}^{\frac{m}{dT}}ij$

我们发现$\sum_{i=1}^{\frac {n}{dT}}\sum_{j=1}^{\frac{m}{dT}}ij$这一团是可以$O(1)$求的
（$\sum_{i=1}^{n}\sum_{j=1}^{m}ij=n*(n+1)/2*m*(m+1)/2$）

又$\mu(T)T^2$可以线性筛$\mu$时求出

就可以先整除分块一下$d$，得到$\frac n d,\frac m d$的值后在里面套一个整除分块求

$\sum_{T=1}^{min(\frac n d,\frac m d)}\mu(T)T^2\sum_{i=1}^{\frac {n}{dT}}\sum_{j=1}^{\frac{m}{dT}}ij$这一团

总复杂度$O(\sqrt n *\sqrt n)=O(n)$

#include<bits/stdc++.h>
using namespace std;
#define int long long
inline int read(){
	char ch=getchar();
	int res=0,f=1;
	while(!isdigit(ch)){if(ch=='-')f=-f;ch=getchar();}
	while(isdigit(ch))res=(res<<3)+(res<<1)+(ch^48),ch=getchar();
	return res*f;
}
const int N=10000005;
const int mod=20101009;
int mu[N],pr[N],vis[N],sum[N],tot,ans;
inline void init(){
	mu[1]=1;
	for(int i=2;i<N;i++){
		if(!vis[i])pr[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot&&i*pr[j]<N;j++){
			vis[pr[j]*i]=1;
			if(i%pr[j]==0)break;
			mu[i*pr[j]]=-mu[i];
		}
	}
	for(int i=1;i<N;i++){
		sum[i]=(sum[i-1]+(i*i)*mu[i])%mod;
	}
	//for(int i=20;i<=30;i++)cout<<sum[i]<<'\n';
}
inline int t(int n,int m){
	return (((n*(n+1)/2)%mod)*((m*(m+1)/2)%mod)%mod);
}
inline int calc(int n,int m){
	int p=min(n,m),res=0;
	for(int i=1,nxt;i<=p;i=nxt+1){
		nxt=min(n/(n/i),m/(m/i));
		res=(res+((sum[nxt]-sum[i-1]+mod)%mod)*t(n/i,m/i)%mod)%mod;
	}
	return res;
}
signed main(){
	init();
	int n=read(),m=read();
	int p=min(n,m);
	for(int i=1,nxt;i<=p;i=nxt+1){
		nxt=min(n/(n/i),m/(m/i));
		ans=(ans+(((nxt-i+1)*(nxt+i)/2)%mod)*calc(n/i,m/i)%mod)%mod;
	}
	cout<<ans<<'\n';
}