题解 P5435 【基于值域预处理的快速 GCD】
moongazer
2020-01-24 00:35:24
# 一种$O($值域$)$时间预处理$O(1)$时间求最大公约数($\gcd$)的算法
[$\mathfrak{View\space it\space on\space my\space Blog}$](https://blog.seniorious.cc/2020/quick-gcd/)
Update in 2020.08.07：修正了一处错误，感谢@Kinesis的指正
## 一些约定
1. $N$为询问的值域
2. $Prime$为全体素数集合
3. 集合$\{a_1,a_2,\cdots,a_m\}$是$n$的分解，当且仅当$a_1\times a_2\times\cdots\times a_m=n$
## 原理
### 定理一
#### 内容
可以将值域中的每个$x$分解成$\{a,b,c\}$，满足$a,b,c\leq\sqrt{x}$或$\in Prime$（定义这种分解为合法分解）
#### 证明
不妨设$a\leq b\leq c$若$c\notin Prime$且$c>\sqrt{x}$，则$c$可分解为$\{d,e\}$且$d\leq e$有$d\leq\sqrt{x}$,而$a\times b=\frac{x}{c}<\sqrt{x}$则有$n$的分解$\{d,a\times b,e\}$，若$e>\sqrt{x}$则可按该规律一直分解直到$e\in Prime$或$\leq\sqrt{x}$
### 定理二
#### 内容
对于询问$\gcd(x,y)$，分别考虑$a,b,c$对答案的贡献，$a$对答案的贡献为$\gcd(a,y)$，再将$y$除以$\gcd(a,y)$（这个因子已经被算过，不能重复计算），再对$b,c$干同样的事，三个贡献相乘即为$\gcd(x,y)$
#### 证明
> 易得引理若$r\mid p, r\mid q$则$\gcd(p,q)=r\times\gcd(\frac{p}{r},\frac{q}{r})$
分别代入$\left\{
\begin{aligned}
&p_1=a\times b\times c,q_1=y,r_1=\gcd(a,q_1) \\
&p_2=b\times c,q_2=\frac{q_1}{r_1},r_2=\gcd(b,q_2) \\
&p_3=c,q_3=\frac{q_2}{r_2},r_3=\gcd(c,q_3)
\end{aligned}
\right.
$即可
## 实现
我们发现实现的难点在于如何在$O(N)$时间内进行分解，询问部分较为容易
### 分解
对于$x\geq2$，找到$x$的最小质因子$p$以及$\frac{x}{p}$的合法分解$\{a_0,b_0,c_0\}$且$a_0\leq b_0\leq c_0$,$x$的一种合法分解即为$\{a_0\times p,b_0,c_0\}$的升序排序
考虑证明：
1. $x\in Prime$时显然成立,分解为$\{1,1,x\}$
2. 当$p\le\sqrt[4]{x}$时将$a_0\leq\sqrt[3]{\frac{x}{p}}$带入有$a_0\times p\le\sqrt{x}$
3. 考虑$p>\sqrt[4]{x}$的情况，$\left(1.\right)$ $a_0=1$，显然有$a_0\times p=p\le\sqrt{x}$；$\left(2.\right)$ $a\neq1$，由于$x$不是素数，$\frac{x}{p}$的最小质因子$q$即为$x$的第二小质因子，一定$\geq p$，而$a_0,b_0,c_0$都为$\frac{x}{p}$的非$1$因子，有$p\leq q\leq a_0\leq b_0\leq c_0$,$p\times a_0\times b_0\times c_0>(\sqrt[4]{x})^4=x$与$p\times a_0\times b_0\times c_0=x$相矛盾，故不存在此情况
所以只用跑一次线性筛，用最小质因子更新即可，然后预处理出$\sqrt{n}\times\sqrt{n}$的$\gcd$数组
代码如下
```cpp
// fac为合法分解，isp表示是否非质数，pri为质数数组，tot为pri的大小，pre为预处理的gcd数组，M为值域，T为sqrt(M)
void work() {
fac[1][0] = fac[1][1] = fac[1][2] = 1;
for (int i = 2; i <= M; ++i) {
if (!isp[i]) {
fac[i][0] = fac[i][1] = 1;
fac[i][2] = i;
pri[++tot] = i;
}
for (int j = 1; pri[j] * i <= M; ++j) {
int tmp = pri[j] * i;
isp[tmp] = true;
fac[tmp][0] = fac[i][0] * pri[j];
fac[tmp][1] = fac[i][1];
fac[tmp][2] = fac[i][2];
if (fac[tmp][0] > fac[tmp][1]) {
fac[tmp][0] ^= fac[tmp][1] ^= fac[tmp][0] ^= fac[tmp][1];
}
if (fac[tmp][1] > fac[tmp][2]) {
fac[tmp][1] ^= fac[tmp][2] ^= fac[tmp][1] ^= fac[tmp][2];
}
// 对于整数运算a ^= b ^= a ^= b等价于swap(a, b)这里就是手动进行length = 3的排序
if (i % pri[j] == 0) {
break;
}
}
}
for (int i = 0; i <= T; ++i) {
pre[0][i] = pre[i][0] = i;
}
for (int i = 1; i <= T; ++i) {
for (int j = 1; j <= i; ++j) {
pre[i][j] = pre[j][i] = pre[j][i % j];
}
}
}
```
### 查询
若当前枚举的为$a$，只需注意$a>\sqrt{N}$时分$a\mid y$和$a\nmid y$两种情况即可
以下为代码
```cpp
int gcd(int a, int b) {
// 不想写if-else所以用三目运算符代替并缩进了一下
int ans = 1;
for (int i = 0; i < 3; ++i) {
int tmp = (fac[a][i] > T) ?
(b % fac[a][i] ?
1
: fac[a][i]
)
: pre[fac[a][i]][b % fac[a][i]];
b /= tmp;
ans *= tmp;
}
return ans;
}
```
## 本题做法
基本上已经没了，就是求$\gcd$然后按题目给的算，放个代码吧
```cpp
const int N = 5000;
const int M = 1000000;
const int T = 1000;
const int Mod = 998244353;
void work();
int gcd(int, int);
int pre[T + 2][T + 2];
int a[N + 2], b[N + 2];
int fac[M + 2][3];
bool isp[M + 2];
int pri[M / 10], tot;
int n;
int main () {
work();
read(n);
for (int i = 1; i <= n; ++i) {
read(a[i]);
}
for (int i = 1; i <= n; ++i) {
read(b[i]);
}
for (int i = 1; i <= n; ++i) {
int ans = 0;
for (int j = 1, now = i; j <= n; ++j, now = now * 1ll * i % Mod) {
ans = (ans + int(now * 1ll * gcd(a[i], b[j]) % Mod)) % Mod;
}
write(ans), EL;
}
return 0;
}
void work() {
fac[1][0] = fac[1][1] = fac[1][2] = 1;
for (int i = 2; i <= M; ++i) {
if (!isp[i]) {
fac[i][0] = fac[i][1] = 1;
fac[i][2] = i;
pri[++tot] = i;
}
for (int j = 1; pri[j] * i <= M; ++j) {
int tmp = pri[j] * i;
isp[tmp] = true;
fac[tmp][0] = fac[i][0] * pri[j];
fac[tmp][1] = fac[i][1];
fac[tmp][2] = fac[i][2];
if (fac[tmp][0] > fac[tmp][1]) {
fac[tmp][0] ^= fac[tmp][1] ^= fac[tmp][0] ^= fac[tmp][1];
}
if (fac[tmp][1] > fac[tmp][2]) {
fac[tmp][1] ^= fac[tmp][2] ^= fac[tmp][1] ^= fac[tmp][2];
}
if (i % pri[j] == 0) {
break;
}
}
}
for (int i = 0; i <= T; ++i) {
pre[0][i] = pre[i][0] = i;
}
for (int i = 1; i <= T; ++i) {
for (int j = 1; j <= i; ++j) {
pre[i][j] = pre[j][i] = pre[j][i % j];
}
}
}
int gcd(int a, int b) {
int ans = 1;
for (int i = 0; i < 3; ++i) {
int tmp = (fac[a][i] > T) ?
(b % fac[a][i] ?
1
: fac[a][i]
)
: pre[fac[a][i]][b % fac[a][i]];
b /= tmp;
ans *= tmp;
}
return ans;
}
```