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

传送门

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

---

$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)}dsum_{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)}dsum_{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)}dsum_{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)}dsum_{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)}dsum_{T=1}^{min(frac n d,frac m d)}mu(T)T^2sum_{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^2sum_{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]<<'
';
}
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<<'
';
}