BZOJ2820 YY的GCD <莫比乌斯反演>

Problem

YY的GCD

$\mathrm{Time\;Limit:\;10\;Sec}$
$\mathrm{Memory\;Limit:\;512\;MB}$

Description

神犇$\mathrm{YY}$虐完数论后给$\mathrm{傻叉kAc}$出了一题：
给定$N,M$，求$1\le x\le N,\;1\le y\le M$且$\gcd(x, y)$为质数的$(x, y)$有多少对
$\mathrm{kAc}$这种$\mathrm{傻叉}$必然不会了，于是向你来请教。
多组输入。

Input

第一行一个整数$T$表述数据组数。
接下来$T$行，每行两个正整数，表示$N,M$。

Output

$T$行，每行一个整数表示第$i$组数据的结果。

Sample Input

2
10 10
100 100

Sample Output

30
2791

HINT

$T=10^4$
$N,M\le10^7$

标签：莫比乌斯反演

Solution

套路反演+积性函数预处理。

先套路推一波反演式：

$$\begin{aligned} Ans&=\sum_{p}\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=p]\\ &=\sum_{p}\sum_{i=1}^{\lfloor n/p\rfloor}\sum_{j=1}^{\lfloor m/p\rfloor}\sum_{d|i,j}\mu(d)\\ &=\sum_{p}\sum_{d=1}^{\lfloor\frac{\min(n,m)}{p}\rfloor}\mu(d)\times\lfloor\frac{n}{d\cdot p}\rfloor\times\lfloor\frac{m}{d\cdot p}\rfloor\\ &=\sum_{T=1}^{n}\sum_{p|T}\mu(\frac{T}{p})\times\lfloor\frac{n}{T}\rfloor\times\lfloor\frac{m}{T}\rfloor\\ &=\sum_{T=1}^{n}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{p|T}\mu(\frac{T}{p})\\ \\ Let\;f(x)&=\sum_{p|T}\mu(\frac{T}{p}),\\ then\;Ans&=\sum_{T=1}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor f(T) \end{aligned}$$

那么我们需要用线筛预处理出所有$f$值。
对于当前预处理到的数$x$和枚举到的质数$p$，有

$$f(x)=\sum_{d|x,\;d\in\;pri}\mu(\frac{x}{d})\\ f(x\cdot p)=\sum_{d|x\cdot p,\;d\in\;pri}\mu(\frac{x\cdot p}{d})\\$$

下面的式子的一部分可以化为上面的式子，剩余部分直接加上去。
分类讨论：

• 若$p|x$，则下式中$d\ne p$时，$\frac{x\cdot p}{d}$的分解式中$p$的指数一定大于$1$，于是只有$d=p$时会对答案产生$\mu(\frac{x\cdot p}{p})=\mu(x)$的贡献，所以$f(x\cdot p)=\mu(x)$。
• 若$p\nmid x$，则下式中对于任意$d|x$，其贡献都为上式中对应项的贡献乘$\mu(p)=-1$，即$d|x$的贡献为$f(x)$。当$d\nmid x$时，只存在$d=p$，此时贡献为$\mu(x)$。因此$f(x\cdot p)=\mu(x)-f(x)$。

由此，可以线筛预处理出所有$f$值，然后根号分块计算答案即可。

Code

#include <bits/stdc++.h>
#define MAX_N 100000
using namespace std;
typedef long long lnt;
template <class T>
inline void read(T &x) {
    x = 0; int c = getchar(), f = 1;
    for (; !isdigit(c); c = getchar()) if (c == 45) f = -1;
    for (; isdigit(c); c = getchar()) (x *= 10) += f*(c-'0');
}
int pri[MAX_N+5], mu[MAX_N+5], fac[MAX_N+5], cnt;    bool NotPri[MAX_N+5];
void PriS() {
    mu[1] = 1;
    for (int i = 2; i <= MAX_N; i++) {
        if (!NotPri[i]) pri[cnt++] = i, mu[i] = -1;
        for (int j = 0, x; j < cnt; j++) {
            if ((x = i*pri[j]) > MAX_N) break;
            NotPri[x] = true;
            if (i%pri[j]) mu[x] = -mu[i];
            else {mu[x] = 0;	break;}
        }
    }
}
int main() {
    int n;	read(n), PriS();	lnt ans = 0LL;
    for (lnt i = 1, l, r; i < sqrt(n); i++) {
        cnt = 0; for (int j = 1; j <= sqrt(i); j++) if (i%j == 0) fac[cnt++] = j;
        for (l = 1; l < i; l = r+1) {
            lnt val = n/(i*(i+l)); if (!val) break;	r = min(n/val/i-i, i-1);
            for (lnt k = 0, j; k < cnt; k++) {
                j = fac[k], ans += mu[j]*(r/j-(l-1)/j)*val;	lnt t = i/j;
                if (i%t == 0 && (j^t)) ans += mu[t]*(r/t-(l-1)/t)*val;
            }
        }
    }
    return printf("%lld\n", ans), 0;
}