BZOJ2154: Crash的数字表格【莫比乌斯反演】

2154: Crash的数字表格

${\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)$

${\sum }_{i=1}^n{\sum }_{j=1}^m\frac{i\ast j}{(i,j)}$

设d=(i,j)

$F(x,y)={\sum }_{i\in [1,x]}{\sum }_{j\in [1,y]}i\ast j\ast [(i,j)==1]$

则有$Ans={\sum }_{d=1}^{min(n,m)}d\ast F(\frac{n}{d},\frac{m}{d})$

另$Sum(x,y)={\sum }_{i\in [1,x]}{\sum }_{y\in [1,y]}i\ast j=\frac{x(x+1)\ast y(y+1)}{4}$

$F(x,y)={\sum }_{i=1}^{min(x,y)}{i}^2\ast \mu (i)\ast Sum(\frac{x}{i},\frac{y}{i})$

F可以数论分块，Ans也可以，复杂度$O(n)$

#include<cstdio>
#include<algorithm>
using namespace std;
const int MAXN=1e7+5,MOD=20101009;
typedef long long LL;
int p[MAXN/10],n,m,mu[MAXN];bool vis[MAXN];
void make_p(){
	vis[0]=vis[1]=1;mu[1]=1;mu[0]=0;
	for(int i=2;i<=m;i++){
		if(!vis[i]) mu[i]=-1,p[++p[0]]=i;
		for(int j=1;j<=p[0]&&i*p[j]<=m;j++)
		if(vis[i*p[j]]=1,i%p[j]==0){mu[i*p[j]]=0;break;}
		else mu[i*p[j]]-=mu[i];
	}
}
int Sum(int x){return 1ll*(x+1)*x/2%MOD;}
int Solve(int x,int y){
	if(x>y) swap(x,y);int Ans=0;
	for(int i=1,j;i<=x;i=j+1){
		j=min(x/(x/i),y/(y/i));
		Ans=(Ans+1ll*((mu[j]-mu[i-1])%MOD+MOD)%MOD*Sum(x/i)%MOD*Sum(y/i)%MOD)%MOD;
	}
	return Ans;
}
int main(){
	scanf("%d%d",&n,&m);
	if(n>m) swap(n,m);make_p();
	for(int i=1;i<=1e7;i++) mu[i]=1ll*mu[i]*i%MOD*i%MOD;
	for(int i=1;i<=1e7;i++) mu[i]=(mu[i]+mu[i-1])%MOD;
	int Ans=0;
	for(int i=1,j;i<=n;i=j+1){
		j=min(n/(n/i),m/(m/i));
		Ans=(Ans+1ll*(i+j)*(j-i+1)/2%MOD*Solve(n/i,m/i)%MOD)%MOD;
	}
	printf("%d\n",Ans);
	return 0;
}