●BZOJ 2694 Lcm

题链：

http://www.lydsy.com/JudgeOnline/problem.php?id=2694

题解：

莫比乌斯反演

不难看出，造成贡献的(i,j)满足gcd(i,j)无平方因子。

其实也就是$mu(gcd(i,j))!=0$

先列出求ANS的式子

$egin{align*}ANS&=sum_{a=1}^{A}sum_{b=1}^{B} lcm(a,b)mu(gcd(a,b))^2;(同样的，先枚举gcd的值g)\&=sum_{g=1}^{min(A,B)} mu(g)^2 imes gsum_{d=1}^{min(frac{A}{g},frac{B}{g})}mu(d) d^2 imes sum(lfloor frac{A}{gd} floor)sum(lfloor frac{B}{gd} floor)\&(sum(x)=frac{(1+n)n}{2})end{align*}$

上式的$gsum_{d=1}^{min(frac{A}{g},frac{B}{g})}mu(d)d^2 imes sum(lfloor frac{A}{gd} floor)sum(lfloor frac{B}{gd} floor)$是求满足gcd(i,j)=g的lcm(i,j)之和，详见●BZOJ 2154 Crash的数字表格。

我们继续：

$egin{align*}ANS&=sum_{g=1}^{min(A,B)} mu(g)^2 imes gsum_{d=1}^{min(frac{A}{g},frac{B}{g})}mu(d) d^2 imes sum(lfloor frac{A}{gd} floor)sum(lfloor frac{B}{gd} floor)\&=sum_{D=gd=1}^{min(A,B)}sum(lfloor frac{A}{D} floor)sum(lfloor frac{B}{D} floor)sum_{g|D}mu(g)^2gcdotmu(frac{D}{g})(frac{D}{g})^2end{align*}$

令$egin{align*}w(D)=sum_{g|D}mu(g)^2gcdotmu(frac{D}{g})(frac{D}{g})^2end{align*}$

现在，如果能够求出w(D)，那么每个询问就可以在$O(sqrt N)$里完成。

由于$y=mu(x)，y=x$是积性函数，那么由狄利克雷乘积的性质可知，

w(D)也是一个积性函数，所以线筛就可以求出w(D)。

代码：

#include<cstdio>
#include<cstring>
#include<iostream>
#define MAXN 4000050
using namespace std;
const int mod=1<<30;
int w[MAXN];
void Sieve(){
	static bool np[MAXN];
	static int prime[MAXN],pnt;
	w[1]=1;
	for(int i=2,tmp,d;i<=4000000;i++){
		if(!np[i]) prime[++pnt]=i,w[i]=(1ll*i-1ll*i*i%mod+mod)%mod;
		for(int j=1;j<=pnt&&i<=4000000/prime[j];j++){
			np[i*prime[j]]=1; tmp=i; d=prime[j];
			while(tmp%prime[j]==0) tmp/=prime[j],d*=prime[j];
			if(tmp!=1) w[tmp*d]=1ll*w[tmp]*w[d]%mod;
			else if(1ll*d==1ll*prime[j]*prime[j]) w[d]=(-1ll*prime[j]*prime[j]%mod*prime[j]%mod+mod)%mod;
			if(i%prime[j]==0) break;
		}
	}
	for(int i=2;i<=4000000;i++) w[i]=(1ll*w[i]+w[i-1])%mod;
}
int sum(int n){
	return 1ll*(1+n)*n/2%mod;
}
int main(){
	Sieve();
	int Case,n,m,mini,ans;
	scanf("%d",&Case);
	while(Case--){
		scanf("%d%d",&n,&m);
		mini=min(n,m); ans=0;
		for(int D=1,last;D<=mini;D=last+1){
			last=min(n/(n/D),m/(m/D));
			ans=(1ll*ans+1ll*((w[last]-w[D-1]+mod)%mod)*sum(n/D)*sum(m/D))%mod;
		}
		printf("%d
",ans);
	}

	return 0;
}