●BZOJ 3309 DZY Loves Math

题链：

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

题解：

莫比乌斯反演，线筛

化一化式子：

f(x)表示x的质因子分解中的最大幂指数

$sum_{i=1}^n sum_{j=1}^m f(gcd(i,j))$

$quadquad=sum_{g=1}^{n}f(g)sum_{d=1}^{lfloor frac{n}{g} floor} mu(d)lfloor frac{n}{gd} floorlfloor frac{m}{gd} floor$

$quadquad=sum_{D=gd=1}^{n}(lfloor frac{n}{D} floorlfloor frac{m}{D} floor)sum_{g|D} f(g)u(frac{D}{g})$

令 $w[D]=sum_{g|D} f(g)u(frac{D}{g})$

然后如果能够预处理出w[D]，那么这个题的每个询问就可以在$O(sqrt N)$的复杂度内解决。

虽然w[D]不是积性函数，但仍可以在线筛时求出，详见BZOJ 3309: DZY Loves Math [莫比乌斯反演 线性筛]。

代码：

#include<cstdio>
#include<cstring>
#include<iostream>
#define MAXN 10000007
using namespace std;
int g[MAXN];
void Sieve(){
	static bool np[MAXN];
	static int prime[MAXN],idx[MAXN],hav[MAXN],pnt;
	for(int i=2,tmp,d;i<=10000000;i++){
		if(!np[i]) prime[++pnt]=i,hav[i]=1,idx[i]=1,g[i]=1;
		for(int j=1;j<=pnt&&i<=10000000/prime[j];j++){
			np[i*prime[j]]=1; hav[i*prime[j]]=hav[i]+(i%prime[j]!=0);

			d=1; tmp=i; while(tmp%prime[j]==0) d++,tmp/=prime[j];
			if(idx[tmp]==d||tmp==1) idx[i*prime[j]]=d;
			if(tmp==1) g[i*prime[j]]=1;
			else if(idx[i*prime[j]]) g[i*prime[j]]=-1*(hav[i*prime[j]]&1?-1:1);
			
			if(i%prime[j]==0) break;
		}
	}
	for(int i=1;i<=10000000;i++) g[i]+=g[i-1];
}
int main(){
	Sieve();
	int Case,n,m,mini; long long ans;
	scanf("%d",&Case);
	while(Case--){
		scanf("%d%d",&n,&m); 
		ans=0; mini=min(n,m);
		for(int D=1,last;D<=mini;D=last+1){
			last=min(n/(n/D),m/(m/D));
			ans+=1ll*(g[last]-g[D-1])*(n/D)*(m/D);
		}
		printf("%lld
",ans);
	}
	return 0;
}