LGOJ2257 YY的GCD

这题美好体验就是卡常

Description

link

求

[sum^{n}_ {i=1} sum^m_{j=1} [gcd(i,j)in prime] ]

其中：(1leq n,m leq 10^7)，多组询问

Solution

对于这种与(gcd)相关的反演题，有一个好的套路

设(f(d)=[gcd(i,j)=d]),(F(n))为(gcd(i,j)=d)和(d)的倍数的个数，即：

[f(d)=[gcd(i,j)=d] ]

[F(n)=sum_{d|n} f(d)= lfloor frac{N}{n} floor lfloor frac{M}{m} floor ]

[f(n)=sum_{n|d} mu(lfloor frac{d}{n} floor)F(d) ]

由着这个套路,我们开始化简这个式子

[sum_{p in prime} sum_{i=1}^n sum_{j=1}^m [gcd(i,j)==p] ]

将(f(p))带入：

[sum_{p in prime}f(p) ]

把(f(x))换成(F(x))

[sum_{p in prime}sum_{p|d} mu(lfloor frac{d}{p} floor)F(d) ]

我们枚举(lfloor frac{d}{p} floor)

[sum_{p in prime} sum_{d=1}^{min(lfloor frac{n}{p} floor,lfloor frac{m}{p} floor} mu(d)F(dp) ]

再把(F(dp))换成最终式：

[sum_{p in prime} sum_{d=1}^{min(lfloor frac{n}{p} floor,lfloor frac{m}{p} floor)} mu(d)lfloor frac{N}{n} floor lfloor frac{M}{m} floor ]

令(T=dp)，则有：

[sum^{min(n,m)}_ {T=1} sum_{t|T,tin prime} mu(lfloor frac{T}{t} floor)lfloor frac{n}{T} floorlfloor frac{m}{T} floor ]

[sum^{min(n,m)}_ {T=1} lfloor frac{n}{T} floorlfloor frac{m}{T} floor(sum_{t|T} mu(lfloor frac{T}{t} floor)) ]

推到这里，我们就都可以做了

(mu(space))可以线性筛，其他的可以整除分块

CODE

#include<bits/stdc++.h>
#define reg register
using namespace std;
namespace yspm{
	inline int read()
	{
		int res=0,f=1; char k;
		while(!isdigit(k=getchar())) if(k=='-') f=-1;
		while(isdigit(k)) res=res*10+k-'0',k=getchar();
		return res*f;
	}
	const int N=10000010;
	bool vis[N]; int pri[N],mu[N],g[N],cnt;
	#define ll long long
	ll sum[N];
	inline void prework()
	{
		mu[1]=1; 
		for(reg int i=2;i<N;++i)
		{
			if(!vis[i]){mu[i]=-1;pri[++cnt]=i;}
			for(reg int j=1;j<=cnt&&pri[j]*i<N;++j)
			{
				vis[i*pri[j]]=1;
				if(i%pri[j]==0) break;
				else mu[pri[j]*i]-=mu[i];
			}
		}
		for(reg int j=1;j<=cnt;++j)
		{
			for(reg int i=1;i*pri[j]<N;++i) g[i*pri[j]]+=mu[i];
		}
		for(reg int i=1;i<N;++i) sum[i]=sum[i-1]+(ll)g[i];
		return ;
	}
	inline void work()
	{
		int n=read(),m=read(); if(n>m) swap(n,m);
		ll ans=0;
		 for(reg int l=1,r;l<=n;l=r+1)
		 {
		 	r=min(n/(n/l),m/(m/l));
		 	ans+=1ll*(n/l)*(m/l)*(sum[r]-sum[l-1]); 
		 }printf("%lld
",ans);
		return ;
	}
	signed main()
	{
		prework(); int T=read(); while(T--) work();
		return 0;
	}
}
signed main(){return yspm::main();}