BZOJ2820 YY的GCD题解（Mobius反演+除法分块）

题目：BZOJ2820.
题目大意：求有多少对$x,y$满足$1\le x\le n,1\le y\le m$且$gcd(x,y)$为质数.
$1\le n,m\le 10^7$，数据组数$=10^4$.

设素数集为P，那么题目要求即为：
$$\sum _{p\in P}\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=p]$$

钦定$n\le m$，设$a=\lfloor \frac{n}{p}\rfloor ,b=\lfloor \frac{m}{p}\rfloor$，那么：
$$\begin{gathered} \sum _{p\in P}\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=p] \\ =\sum _{p\in P}\sum _{i=1}^a\sum _{j=1}^b[gcd(i,j)=1] \\ =\sum _{p\in P}\sum _{i=1}^a\sum _{j=1}^b\sum _{d|i\wedge d|j}\mu (d) \\ =\sum _{d=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{p\in P}\mu (d)\lfloor \frac{n}{pd}\rfloor \lfloor \frac{m}{pd}\rfloor \end{gathered}$$

设$k=pd$，那么：
$$\begin{gathered} \sum _{p\in P}\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=p] \\ =\sum _{k=1}^n\sum _{p\in P\wedge p|k}\mu (\frac{k}{p})\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor \end{gathered}$$

设一个函数：$f(k)={\sum }_{p\in P\wedge p|k}\mu (\frac{k}{p})$，那么：
$$\sum _{p\in P}\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=p]=\sum _{k=1}^{n}f(k)\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor$$

发现这个式子中，只要$f$函数能够用线性筛筛出即可用除法分块做到$O(\sqrt{n}+\sqrt{m})$求解，所以尝试推导这个函数在线性筛三种情况下的值：
1.当$i$为质数时，显然$f(i)=1$.
2.当$gcd(i,pr[j])=1$时，可以推导：
KaTeX parse error: Expected '}', got '\wdge' at position 27: …\\=\sum_{p\in P\̲w̲d̲g̲e̲ ̲p|(i*pr[j])}\mu…

3.当$pr[j]|i$时，可以推导：
$$\begin{gathered} f(i\ast pr[j]) \\ =\sum _{p\in P\wedge p|(i\ast pr[j])}\mu (\frac{i\ast pr[j]}{p}) \\ =\mu (\frac{i\ast pr[j]}{pr[j]})+\sum _{p\in P\wedge p|i}\mu (\frac{i\ast pr[j]}{p}) \\ =\mu (i) \end{gathered}$$

发现我们就这么把$f$函数筛出来了，所以除法分块一下就可以了，时间复杂度$O(m+T\sqrt{m})$.

代码如下：

#include<bits/stdc++.h>
  using namespace std;

#define Abigail inline void
typedef long long LL;

const int N=10000000;

int n,m;
int pr[N+9],tp,b[N+9];
LL mu[N+9],f[N+9];

void sieve(int n){
  for (int i=2;i<=n;++i) b[i]=1;
  mu[1]=1;
  for (int i=2;i<=n;++i){
    if (b[i]) pr[++tp]=i,mu[i]=-1,f[i]=1;
    for (int j=1;j<=tp&&pr[j]*i<=n;++j){
      b[i*pr[j]]=0;
      if (i%pr[j]==0){
	    mu[i*pr[j]]=0;
	    f[i*pr[j]]=mu[i];
		break;
	  }
	  mu[i*pr[j]]=-mu[i];
	  f[i*pr[j]]=mu[i]-f[i];
    }
  }
  for (int i=1;i<=n;++i)
    f[i]+=f[i-1];
}

LL Query(int n,int m){
  if (n>m) swap(n,m);
  LL ans=0; 
  for (int l=1,r;l<=n;l=r+1){
    r=min(n/(n/l),m/(m/l));
    ans+=(f[r]-f[l-1])*(n/l)*(m/l);
  }
  return ans;
}

Abigail start(){
  sieve(N);
}

Abigail into(){
  scanf("%d%d",&n,&m);
}

Abigail outo(){
  printf("%lld\n",Query(n,m));
}

int main(){
  start();
  int T;
  scanf("%d",&T);
  while (T--){
  	into();
  	outo();
  }
  return 0;
}