YY的GCD(莫比乌斯反演)

YY的GCD

Statement

求${\sum }_{i=1}^N{\sum }_{j=1}^M[\gcd (i,j)=Prime]$.

Solution

设: $f(d)={\sum }_{i=1}^N{\sum }_{j=1}^M[\gcd (i,j)=d]$.

设: $g(d)={\sum }_{d|d^{\prime }}f(d^{\prime })={\sum }_{i=1}^N{\sum }_{j=1}^M[d|i][d|j]={\sum }_{i=1}^N[d|i]{\sum }_{i=1}^M[d|j]={\sum }_{i=1}^N[d|i]\lfloor \frac{M}{d}\rfloor =\lfloor \frac{N}{d}\rfloor \cdot \lfloor \frac{M}{d}\rfloor$.

则: $f(d)={\sum }_{d|d^{\prime }}\mu (\frac{d^{\prime }}{d})g(d^{\prime })$.

由 $Ans={\sum }_{p\in prime}{\sum }_{i=1}^N{\sum }_{j=1}^N[\gcd (i,j)=p]$.

可得: $Ans={\sum }_{p\in prime}f(p)$.
A n s = ∑ p ∈ p r i m e ∑ i = 1 N ∑ j = 1 N [ gcd ⁡ ( i , j ) = p ] = ∑ p ∈ p r i m e f ( p ) = ∑ p ∈ p r i m e ∑ p ∣ d min ⁡ { N , M } μ ( d p ) ⌊ N d ⌋ ⌊ M d ⌋ = ∑ d = 1 min ⁡ { N , M } ∑ p ∣ d 且 p ∈ p r i m e μ ( d p ) ⌊ N d ⌋ ⌊ M d ⌋ = ∑ d = 1 min ⁡ { N , M } ⌊ N d ⌋ ⌊ M d ⌋ ∑ p ∣ d 且 p ∈ p r i e m μ ( d p ) \begin{aligned}Ans&=\sum_{p\in prime}\sum_{i=1}^N\sum_{j=1}^N[\gcd(i,j)=p]\\ &=\sum_{p\in prime}f(p)\\ &=\sum_{p\in prime}\sum_{p|d}^{\min\{N,M\}}\mu(\frac{d}{p})\lfloor\frac{N}{d}\rfloor\lfloor\frac{M}{d}\rfloor\\ &=\sum_{d=1}^{\min\{N,M\}}\sum_{p|d且p\in prime}\mu(\frac{d}{p})\lfloor\frac{N}{d}\rfloor\lfloor\frac{M}{d}\rfloor\\ &=\sum_{d=1}^{\min\{N,M\}}\lfloor\frac{N}{d}\rfloor\lfloor\frac{M}{d}\rfloor\sum_{p|d且p\in priem}\mu(\frac{d}{p}) \end{aligned} Ans​=p∈prime∑​i=1∑N​j=1∑N​[gcd(i,j)=p]=p∈prime∑​f(p)=p∈prime∑​p∣d∑min{N,M}​μ(pd​)⌊dN​⌋⌊dM​⌋=d=1∑min{N,M}​p∣d且p∈prime∑​μ(pd​)⌊dN​⌋⌊dM​⌋=d=1∑min{N,M}​⌊dN​⌋⌊dM​⌋p∣d且p∈priem∑​μ(pd​)​
我们发现$\sum _{p|d且p\in prime}\mu (\frac{d}{p})$是可以预处理的, 该部分我们枚举$p\in prime$即可, 时间复杂度是$O(n\log \log n)$的和$\text{Eratosthenes}$筛法(埃拉托特尼筛法, 简称埃氏筛法)该部分时间复杂度和写法相似.

然后一次询问我们只需要用数论分块就可以将复杂度做到$O(N\log \log N+T\sqrt{N})$.

Code

# define Fast_IO std::ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
# include "unordered_map"
# include "algorithm"
# include "iostream"
# include "cstdlib"
# include "cstring"
# include "cstdio"
# include "vector"
# include "bitset"
# include "queue"
# include "cmath"
# include "ctime"
# include "map"
# include "set"
# define ll long long
# define ld long double
# define rep1(i,a,b) for(ll i=(a);i<=(b);i++)
# define rep2(i,a,b) for(ll i=(b);i>=(a);i--)
# define pii pair<int,int>
# define pll pair<ll,ll>
# define ph push_back
# define pb pop_back
# define eb emplace_back
# define vi vector<int>
# define vll vector<ll>
# define vpi vector<pii >
# define vpll vector<pll >
# define ri(x) scanf("%d",&x)
# define rf(x) scanf("%f",&x)
# define rl(x) scanf("%lld",&x)
# define rd(x) scanf("%lf",&x)
# define rs(s) scanf("%s",s+1)
# define wi(x) printf("%d",x)
# define wl(x) printf("%lld",x)
# define ws(s) printf("%s",s+1)
# define resize(v,x) v.resize(x)
# define all(v) v.begin(),v.end()
# define reverse(v) reverse(all(v))
# define fi first
# define se second
# define lowbit(x) ((x)&(-(x)))
# define repauto(Name,v) for(auto Name:v)
# define Endl "\n"
using namespace std;
template<class I> inline I GCD(I A,I B){return B?GCD(B,A%B):A;}
template<class I> inline I LCM(I A,I B){return A/GCD(A,B)*B;}
template<class I> I Sqrt(I N){
	I sqrtN=sqrt(N)-1;
	while(sqrtN+1<=N/(sqrtN+1))sqrtN++;
	return sqrtN;
}
template<class I> I Pow(I X,I Y,__int128 Mod1=998244353){
	static __int128 Ans; Ans=1;
	for(;Y;Y>>=1,X=(__int128)X*X%Mod1) if(Y&1) Ans=Ans*X%Mod1;
	return Ans;
}

namespace Mu_Class{
	# define maxm 10000010 
	int Prime_Cnt;
	int Prime[maxm];
	int Mu[maxm];
	int Visit[maxm];
	long long F[maxm];
	inline void Init(int N){
		static int i,j;
		for(Mu[1]=1,i=2;i<=N;i++){
			if(!Visit[i]){
				Visit[i]=i,Mu[i]=-1;
				Prime[++Prime_Cnt]=i;
			}for(j=1;j<=Prime_Cnt && Visit[i]>=Prime[j] && Prime[j]<=N/i;j++){
				Mu[i*Prime[j]]=Visit[i]==Prime[j]?0:Mu[i]*(-1);
				Visit[i*Prime[j]]=Prime[j];
			}
		}for(i=1;i<=Prime_Cnt;i++) for(j=1;Prime[i]<=N/j;j++) F[j*Prime[i]]+=Mu[j];
		for(i=1;i<=N;i++) F[i]+=F[i-1];
		return;
	}
	# undef maxm
}using namespace Mu_Class;

int T,N,M;
long long Ans;

int main(){
# ifdef LH_Frank
    freopen("1.in","r",stdin);
	freopen("1.out","w",stdout);
# endif
	Init(10000000);
	ri(T);
	while(T--){
		Ans=0;
		ri(N),ri(M);
		for(int i=1;i<=min(N,M);){
			int j=min(N/(N/i),M/(M/i));
			Ans+=1LL*(N/i)*(M/i)*(F[j]-F[i-1]);
			i=j+1;
		}printf("%lld\n",Ans);
	}
	return 0;
}

LInk

[1] Luogu P2257 YY的GCD