题解 P5176 【公约数】
Nemlit
2019-09-27 22:28:20
~~看到这题目就不想做了系列,出题人是不是都不知道杜教筛是什么东西啊,他家杜教筛可以预处理优化到~~$O(n^\frac{2}{3})$
~~先吓唬你一下,我们要求:~~
$$\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p\gcd(i\cdot j,i\cdot k,j\cdot k)\times \gcd(i,j,k)\times \left(\frac{\gcd(i,j)}{\gcd(i,k)\times \gcd(j,k)}+\frac{\gcd(i,k)}{\gcd(i,j)\times \gcd(j,k)}+\frac{\gcd(j,k)}{\gcd(i,j)\times \gcd(i,k)}\right)$$
我们令$x = gcd(i, j), y = gcd(i, k), z = gcd(j, k)$
于是我们有:
$$\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p\gcd(i\cdot j,i\cdot k,j\cdot k)\times \gcd(i,j,k)\times \left(\frac{x}{y\times z}+\frac{y}{x\times z}+\frac{z}{x\times y}\right)$$
$$\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p\gcd(i\cdot j,i\cdot k,j\cdot k)\times \gcd(i,j,k) \times \frac{x ^ 2 + y ^ 2 + z ^ 2}{x \times y \times z}$$
展开第一陀恶心人的柿子,可以得到:
$$gcd(i\cdot j,i\cdot k,j\cdot k) = \frac{x \times y\times z}{gcd(i, j, k)}$$
(证明可以通过**唯一分解定理**展开)
于是我们惊喜地发现,我们只需要求$($还原$x, y, z)$:
$$\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p gcd(i, j)^2 +gcd(j, k) ^ 2 + gcd(i, k) ^ 2$$
我们把柿子分开写,得到:
$$\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p gcd(i, j)^2+\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^pgcd(j, k) ^ 2 +\sum_{i=1}^n\sum_{j=1}^m\sum_{k=1}^p gcd(i, k) ^ 2$$
因为$gcd(i, j)$与k无关(其他同理),于是我们可以提出k,变成p个$\sum_{i=1}^n\sum_{j=1}^m gcd(i, j)^2$相加,故有:
$$p \times \sum_{i=1}^n\sum_{j=1}^m gcd(i, j)^2+n \times \sum_{j=1}^m\sum_{k=1}^pgcd(j, k) ^ 2 +m \times \sum_{i=1}^n\sum_{k=1}^p gcd(i, k) ^ 2$$
设$F(n, m) = \sum_{i=1}^n\sum_{j=1}^m gcd(i, j)^2$
原式$=p \times F(n, m) + n \times F(p, m) + m \times F(n, p)$
所以问题转化为如何快速的求出$F(n, m)$(不妨设$n ≤ m$)
$$F(n, m) = \sum_{d = 1}^n\sum_{i=1}^n\sum_{j=1}^m d^2 \times [gcd(i, j) == d]$$
$i, j$同时除以d得:
$$F(n, m) = \sum_{d = 1} ^ n\sum^{n / d}_{i = 1}\sum^{m / d}_{j = 1} d ^ 2 \times [gcd(i, j) == 1]$$
$$F(n, m) = \sum_{d = 1} ^ nd ^ 2\sum^{n / d}_{i = 1}\sum^{m / d}_{j = 1}\sum_{k|gcd(i, j)} \mu(k)$$
提前枚举K得:
$$F(n, m) = \sum_{d = 1} ^ n d ^ 2\sum_{k = 1}^{n / d}\mu(k)\sum^{n / (d * k)}_{i = 1}\sum^{m / (d * k)}_{j = 1} $$
$$F(n, m) = \sum_{d = 1} ^ n d ^ 2\sum_{k = 1}^{n / d}\mu(k)\left \lfloor \frac{n}{d * k} \right \rfloor \times \left \lfloor \frac{m}{d * k} \right \rfloor$$
跟据反演套路,我们令$T = d * k$
$$F(n, m) = \sum_{T = 1}^n \left \lfloor \frac{n}{T} \right \rfloor \times \left \lfloor \frac{m}{T} \right \rfloor\sum_{d | T} d ^ 2\times \mu(\frac{T}{d})$$
令$g(x) = x * x, mu(x) = \mu(x)$
于是$F(n, m) = \sum_{T = 1}^n \left \lfloor \frac{n}{T} \right \rfloor \times \left \lfloor \frac{m}{T} \right \rfloor\times(g(x) \times mu(x))$
发现后面那一堆就是一个卷积的形式,而且是两个积性函数的卷积,于是我们可以用线性筛筛出,筛法见[于神之怒加强版](https://www.luogu.org/problem/P4449)
总结一下求法:
我们记$prim[i]$为第i个质数,$low[i]$为i质因数分解后最小的质因子的指数次幂,$f[x]$为$g(x)$卷$mu(x)$,$mu[x]$为$\mu(x)$
类似于线性筛,我们分三种情况讨论:
$case1:$ $i$是质数
我们可以直接算出来:$low[i] = i,\;\ f[i] = i \times i - 1,\;\ mu[i] = -1$
$case2:i \% prim[j] != 0$:这种情况是满足积性函数定义的,即$i, prim[j]$互质,令$p = i \times prim[j]$,所以我们有$low[p] = prim[j],\;\ mu[p] = -mu[i],\;\ f[i] = f[i] \times f[prim[j]]$
$case3:i \% prim[j] == 0$:这种情况不满足积性函数的定义,于是我们考虑用我们刚刚维护的$low[i]$
考虑到$i / low[i],prim[j]$是互质的,那我们可不可以类比$case2$通过$f[i / low[i]], f[prim[j]]$来推出$f[p]$呢?
假设$a, b$互质,我们有:$f[a \times b] = f[a] \times f[b]$
所以我们也有:$f[p] = f[i] \times f[prim[j]] = f[i\ /\ low[i]] \times f[prim[j] \times low[i]]$
根据定义,$i\ /\ low[i]$没有$prim[j]$这个约数,所以$i / low[i],\;\ prim[j] \times low[i]$互质
i肯定有$prim[j]$这个质因数,而且根据线性筛的性质,$prim[j]$肯定是i的最小质因数,故我们有$mu[p] = 0,\;\ low[p] = low[i] \times prim[j]$
如果$low[i] < i$,证明$low[i] \times prim[j] <= i$,故我们可以直接用$f[p] = f[i\ /\ low[i]] \times f[prim[j] \times low[i]]$
如果$low[i] == i$,证明i以前肯定没有被筛过,若$low[i] = prim[j] ^ k$,则k肯定为质数,那么我们就可以直接用:$f[p] = (i * prim[j]) ^ 2 - i ^ 2$
### $ps:$
写完后看题解才发现:这道题根于神之怒不一样,并不需要记录$low[i]$,对于$case3$,$f[p] = f[i] \times prim[j] ^2$即可
# $Code:$
```
#include<bits/stdc++.h>
using namespace std;
#define il inline
#define re register
il int read() {
re int x = 0, f = 1; re char c = getchar();
while(c < '0' || c > '9') { if(c == '-') f = -1; c = getchar();}
while(c >= '0' && c <= '9') x = x * 10 + c - 48, c = getchar();
return x * f;
}
#define rep(i, s, t) for(re int i = s; i <= t; ++ i)
#define mod 1000000007
#define maxn 20000007
il int Pow(int x) {return 1ll * x * x % mod;}
int cnt, prim[maxn], f[maxn], mu[maxn], low[maxn];
bool vis[maxn];
il void init() {
f[1] = mu[1] = 1;
rep(i, 2, 2e7) {
if(!vis[i]) f[i] = Pow(i) - 1, mu[i] = -1, low[i] = i, prim[++ cnt] = i;
for(re int j = 1; j <= cnt; ++ j) {
if(i * prim[j] > 2e7) break;
int p = i * prim[j];
vis[p] = 1;
if(i % prim[j] == 0) {
low[p] = low[i] * prim[j];
if(low[i] == i) f[p] = (Pow(1ll * i * prim[j] % mod) - Pow(i) + mod) % mod;
else f[p] = 1ll * f[i / low[i]] * f[low[i] * prim[j]] % mod;
break;
}
mu[p] = -mu[i], low[p] = prim[j], f[p] = 1ll * f[i] * f[prim[j]] % mod;
}
}
rep(i, 1, 2e7) f[i] = (f[i - 1] + f[i]) % mod;
}
il int F(int n, int m) {
if(n > m) swap(n, m);
int ans = 0;
for(re int l = 1, r; l <= n; l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans = (ans + 1ll * (1ll * (n / l) * (m / l) % mod) * (f[r] - f[l - 1]) % mod) % mod;
}
return (ans + mod) % mod;
}
il int get(int n, int m, int p) {
return ((1ll * p * F(n, m) % mod + 1ll * n * F(p, m) % mod) % mod + 1ll * m * F(n, p)) % mod;
}
il void outit() {
for(re int i = 1, T = read(); i <= T; ++ i) printf("%d\n", get(read(), read(), read()));
}
int main() {
init(), outit();
return 0;
}
```