洛谷 P2257 - YY的GCD（莫比乌斯反演+整除分块）
题意：
- 求满足 \(1 \leq x \leq n\)，\(1 \leq y \leq m\)，\(\gcd(x,y)\) 为质数的数对 \((x,y)\) 的个数。
- \(T\) 组询问。
- \(1 \leq T \leq 10^4\)，\(1 \leq n,m \leq 10^7\)。
今天终于学会了莫比乌斯反演反演~~，就写篇博客加深下印象吧。
要说这莫比乌斯反演有多么博大精深，就不得不从莫比乌斯函数 \(\mu(x)\) 说起。
我们定义 \(\mu(x)\) 为：
当然，\(\mu(x)\) 有不少有趣的性质，其中最重要且最常用的一条是：
为什么？其实很好理解。
当 \(x=1\) 时显然得证。
当 \(x \neq 1\) 时，设 \(n=\prod\limits_{i=1}^kp_i^{c_i},n'=\prod\limits_{i=1}^kp_i\)。
那么显然有 \(\sum\limits_{d|n}\mu(d)=\sum\limits_{d|n'}\mu(d)\)。
而 \(\sum\limits_{d|n'}\mu(d)=\sum\limits_{i=0}^k\dbinom{k}{i}\times (-1)^i\)。
由二项式定理易得原式 \(=0\)。故原命题成立。
那这形式诡异的莫比乌斯函数对我们解题能有什么帮助呢？
这就要引入一个更加神的知识：莫比乌斯反演。
何谓反演？相信不少刚学莫比乌斯反演的人都会有这样的疑问。
回忆起圆的反演。如果存在圆 \(O\) 与点 \(P\) 与 \(P'\) 使得 \(O,P,P'\) 共线且 \(|OP| \times |OP'|=r^2\)，那么就称“\(P\) 与 \(P'\) 互为圆 \(O\) 的反演点”。知道了 \(P\) 的位置，\(P'\) 的位置也确定下来了。反之亦然。
莫比乌斯反演也类似。假设我们有两个函数 \(f(n),F(n)\)，满足 \(F(n)=\sum\limits_{d|n}f(d)\)，那么 \(f(n)=\sum\limits_{d|n}\mu(d)F(\dfrac{n}{d})\)。
怎么证明这个定理呢？最简单的方法是利用狄利克雷卷积来证明。可惜我不会， 5 天后咱再回头来看用狄利克雷卷积的解法。这里介绍一种纯推式子的方法：
同时，还有第二种形式的莫比乌斯反演。如果 \(F(n)=\sum\limits_{n|d}f(n)\)，那么 \(f(n)=\sum\limits_{n|d}\mu(\frac{d}{n})F(d)\)。证明方法与第一种形式的莫比乌斯反演几乎一样，这里就不再赘述了。
说了这么多，相信大家对莫比乌斯反演有了一定的了解。那么莫比乌斯函数及莫比乌斯反演怎样真正应用在题目中呢？我们就以这道题为例来看一看。
很自然地想到枚举 \(\gcd\)，即 \(\sum\limits_{p \in prime}\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=p]\)。
下面就是套路了。令 \(f(d)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=d]\)，\(F(d)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[d|\gcd(i,j)]\)。
显然 \(F(p)=\sum\limits_{p|d}f(d)\)。
套用第二种形式的莫比乌斯反演得 \(f(p)=\sum\limits_{p|d}\mu(\frac{d}{p})F(d)\)。
这样转化有什么好处呢？不难发现要求 \(f(d)\) 非常棘手，但是要求 \(F(d)\) 非常容易————\(F(d)=\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor\)
一不做二不休，把 \(F(d)\) 带回原式得 \(f(p)=\sum\limits_{p|d}\mu(\frac{d}{p})\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor\)
最后计算答案。\(ans=\sum\limits_{p \in prime}f(p)=\sum\limits_{p \in prime}\sum\limits_{p|d}\mu(\frac{d}{p})\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor=\sum\limits_{d=1}^n\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor\sum\limits_{p|d \operatorname{and} p \in prime}\mu(\frac{d}{p})\)
记 \(t_d=\sum\limits_{p|d \operatorname{and} p \in prime}\mu(\frac{d}{p})\)，显然可以在 \(\mathcal O(n \log n)\) 的时间复杂度内求出所有 \(t_d\)。还是把 \(t_d\) 带回去，\(ans=\sum\limits_{d=1}^n\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor t_d\)。
如果单次询问那么的确一遍循环就可以了。可是本题有 \(10^4\) 组数据，出题人显然不会放过 \(\mathcal O(nT)\) 解法。如何进行优化呢？
考虑一个特别经典的问题，求 \(\sum\limits_{i=1}^n\lfloor \frac{n}{i} \rfloor\)。
显然，暴力的 \(\mathcal O(n)\) 的解法是不能被接受的。
不过仔细观察就会发现一个性质：有不少重复的 \(\lfloor \frac{n}{i} \rfloor\)，而且 \(\lfloor \frac{n}{i} \rfloor\) 相同的 \(i\) 是成段分部的。
故我们可以找出 \(\lfloor \frac{n}{i} \rfloor\) 相同的段然后计算贡献。最多有 \(2 \sqrt{n}\) 个这样的段。
这就是著名的整除分块的思想。
回到本题来。使用类似的套路，找出 \(\lfloor \frac{n}{i} \rfloor\) 与 \(\lfloor \frac{m}{i} \rfloor\) 都相同的 \(i\)——它们也一定是成段分部的，使用前缀和计算贡献就完事了。
时间复杂度 \(\mathcal O(T\sqrt{n}+n \log n)\)。
/*
Contest: -
Problem: P2257
Author: tzc_wk
Time: 2020.8.31
*/
#include <bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define pb push_back
#define fz(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
#define foreach(it,v) for(__typeof(v.begin()) it=v.begin();it!=v.end();it++)
#define all(a) a.begin(),a.end()
#define fill0(a) memset(a,0,sizeof(a))
#define fill1(a) memset(a,-1,sizeof(a))
#define fillbig(a) memset(a,0x3f,sizeof(a))
#define y1 y1010101010101
#define y0 y0101010101010
typedef pair<int,int> pii;
typedef long long ll;
inline int read(){
int x=0,neg=1;char c=getchar();
while(!isdigit(c)){
if(c=='-') neg=-1;
c=getchar();
}
while(isdigit(c)) x=x*10+c-'0',c=getchar();
return x*neg;
}
int T=read();
bool vis[10000005];
ll mu[10000005],pri[5000005],sum[10000005],t[10000005],pcnt=0;
inline void get_mu(){
mu[1]=1;
for(int i=2;i<=10000000;i++){
if(!vis[i]){mu[i]=-1;pri[++pcnt]=i;}
for(int j=1;j<=pcnt&&pri[j]*i<=10000000;j++){
vis[i*pri[j]]=1;
if(i%pri[j]==0) break;
else mu[pri[j]*i]=-mu[i];
}
}
for(int i=1;i<=pcnt;i++) for(int j=1;j*pri[i]<=10000000;j++) t[j*pri[i]]+=mu[j];
for(int i=1;i<=10000000;i++) sum[i]=sum[i-1]+t[i];
}
int main(){
get_mu();
while(T--){
int n=read(),m=read();
ll ans=0;
for(int l=1,r;l<=min(n,m);l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=1ll*(n/l)*(m/l)*(sum[r]-sum[l-1]);
}
printf("%lld\n",ans);
}
return 0;
}
