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

传送门

题意：求${\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j),n,m\le 1e7$

---

$Solution$

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

$$ans=\sum _{i=1}^n\sum _{j=1}^m\frac{i\ast 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^{\prime }=\frac{d}{i},j^{\prime }=\frac{d}{j}$(即从枚举$d$的倍数变成枚举是$d$的几倍)
则(为了方便仍用$i,j$表示$i^{\prime },j^{\prime }$)

$$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\ast 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^{\prime }=i\ast T,j^{\prime }=j\ast 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}}ij{T}^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\ast (n+1)/2\ast m\ast (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}\ast \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';
}