[可能有用科技]更多 GCD 算法
前言
重写旧文 & 增加新内容。
欧几里得算法
int gcd(int a,int b) {
return b ? gcd(b,a % b) : a;
}
非递归 ：
int gcd(int a,int b) {
int p = a % b;
while(p)
a = b,b = p,p = a % b;
return b;
}
更相减损法
可半者半之，不可半者，副置分母、子之数，以少减多，更相减损，求其等也。以等数约之。
白话 ：
- 将被 \(2\) 整除的一直除以 \(2\) 直到因数不包含 \(2\)。
- 用较大的数减去较小的，直到减数和差相等，得到的就是最大公因子。
- 同时除了 \(k\) 次 \(2\) 的情况，要在结尾乘以 \(2^k\)
int gcd(int x,int y) {
if(!x || !y) return x | y;
if(!(x & 1) && !(y & 1)) return gcd(x >> 1,y >> 1) << 1;
else if(!(y & 1)) return gcd(x,y >> 1);
else if(!(x & 1)) return gcd(x >> 1,y);
return gcd(abs(x - y),std::min(x,y));
}
非递归 ：
int gcd(int x,int y) {
int p = 0;
while(!(x & 1) && !(y & 1)) {
x >>= 1,y>>= 1;
p++;
}
while(!(x & 1)) x >>= 1;
while(!(y & 1)) y >>= 1;
while(x && y) {
if(x > y) std::swap(x,y);
y = y - x;
}
return (1 << p) * (x | y);
}
能不能再给力一点？
尝试快速判断 \(a,b\) 能除以 \(2\) 几次,众所周知整除 \(2\) 就是二进制末位为 \(1\) 的时候二进制右移一位,那么求出 \(a\) 按位或 \(b\) 的二进制末尾连续 \(0\) 个数即可,这个过程通过跑得飞快的内置函数 __builtin_ctz() 实现.
inline int gcd(int a,int b) {
int k = __builtin_ctz(a | b);
a >>= __builtin_ctz(a);
b >>= __builtin_ctz(b);
while(a && b) {
if(a > b) std::swap(a,b);
b -= a;
}
return a << k;
}
当然可以把这个玩意和欧几里得算法结合一下 ：
inline int gcd(int a,int b) {
int k = __builtin_ctz(a | b);
a >>= __builtin_ctz(a);
b >>= __builtin_ctz(b);
int p = a % b;
while(p)
a = b,b = p,p = a % b;
return b << k;
}
能不能再给力一点？
值域预处理
将值域表示为 \(\mathbb{V}\)，质数集表示为 \(\mathbb{P}\)。
对于值域较小的情况，可以做到 \(\mathcal{O}(\mathbb{V})\) 预处理，\(\mathcal{O}(1)\) 查询。
Theorem 1
\(\forall n \in \mathbf{V}\) 可以被分解为一个三元组 \(\{a,b,c\}\ (a \le b \le c)\) 满足 \(a,b \le \sqrt{n}\)，\(c \le \sqrt{n}\) 或 \(c \in \mathbb{P}\)。
Proof 1
归纳法 ：
基础 ： 对于 \(n = 1\) 有 \(\{1,1,1\}\)。
归纳 ： 假设现在对于一个数 \(n\)，\(x \in [0,n)\) 都有解，尝试证明对于 \(n\) 有解。
首先求 \(n\) 的最小质因数 \(p\)，可知现在 \(\dfrac{n}{p}\) 有解，记为 \(\{a,b,c\}\)。
根据定义 \(n = pabc,\dfrac{n}{p} = abc\)，那么 \(a,b \le \sqrt[3]{\dfrac{n}{p}}\)。
如果 \(p \le \sqrt[4]{n}\)，我们之前得到 ：
然后
并且 \(\sqrt[3]{p^2n} \le \sqrt{n}\)
那么得到 \(pa \le \sqrt{n}\)，有一组合法的分解为 \(\{pa,b,c\}\)
然后是 \(p > \sqrt[4]{n}\) 的情况。
我们知道 \(n = pabc\)，那么显然 \(p,a,b,c > 1\) 时 \(n > p^4\)。
于是这时候 \(a\) 显然是 \(1\)，一组可行的分解为 \(\{p,b,c\}\)。
证毕。
可以发现分解一个数 \(n\) 成三个因数的方式就是首先求其最小质因子 \(p\) 然后求 \(\dfrac{n}{p}\) 的分解方案 \(\{a,b,c\}(a \le b \le c)\)，那么 \(\{pa,b,c\}\) 就是 \(n\) 的分解方案了。
然后线性筛的时候一个数会被其最小质因子筛到，那么显然这个分解可以直接线性预处理的。
最后我们回到如何求 \(\gcd\)。
对于一个 \(\gcd(n,m)\)，将 \(n\) 拆分为 \(\{a,b,c\}\)。
首先算出 \(\gcd(a,m)\) 然后 \(m \leftarrow \dfrac{m}{\gcd(a,m)}\)
然后对 \(b,c\) 重复这个过程即可。
具体令当前细化到求一个 \(\gcd(p,q)\) ：
如果 \(p \le \sqrt{\mathbb{V}}\) 就求 \(\gcd(p,q \bmod p)\) 使得要求出的两个数都在 \(\sqrt{\mathbb{V}}\) 范围内，这些 \(\gcd\) 可以 \(\sqrt{\mathbb{V}} \times \sqrt{\mathbb{V}} = \mathcal{O}(\mathbb{V})\) 预处理出来。
然后如果 \(p > \sqrt{\mathbb{V}}\)，如果 \(p | q\) 那么结果就是 \(p\)，如果 \(p\not|\ q\) 那就是 \(1\)。
最后回顾一下全过程 ：
- 预处理 \(\sqrt{\mathbb{V}}\) 以内的 \(\gcd\)，复杂度 \(\mathcal{O}(\mathbb{V})\)。
- 预处理 \(\mathbb{V}\) 以内的所有数的分解方案，类似线性筛的方式，复杂度 \(\mathcal{O}(\mathbb{V})\)。
- 询问的时候需要的信息都预处理出来了，复杂度 \(\mathcal{O}(1)\)。
于是就有 \(<\mathcal{O}(\mathbb{V}),\mathcal{O}(1)>\) 的 \(\gcd\) 了。
模板 ： \(\mathtt{Link.}\)
Code :
为啥我数学相关的代码都巨大常数啊……
V 是值域常量，SV 是值域根号。
struct CommonFac {
bool vis[V + 5];
int pri[V],pcnt;
int fac[V + 5][3];
int G[SV + 5][SV + 5];
void init() {
rep(i,1,SV) {
G[i][i] = G[i][0] = G[0][i] = i;
repl(j,1,i)
G[i][j] = G[j][i] = G[j][i % j];
}
fac[1][0] = fac[1][1] = fac[1][2] = 1;
rep(i,2,V) {
if(!vis[i]) {
fac[i][0] = fac[i][1] = 1;
fac[i][2] = i;
pri[++pcnt] = i;
}
for(int j = 1;j <= pcnt && i * pri[j] <= V;++j) {
int t = i * pri[j];vis[t] = 1;
fac[t][0] = fac[i][0] * pri[j];
fac[t][1] = fac[i][1],fac[t][2] = fac[i][2];
if(fac[t][0] > fac[t][1]) std::swap(fac[t][0],fac[t][1]);
if(fac[t][1] > fac[t][2]) std::swap(fac[t][1],fac[t][2]);
if(!i % pri[j]) break;
}
}
}
inline int gcd(int a,int b) {
if(a <= SV && b <= SV) return G[a][b];
int res = 1;
repl(i,0,3) {
int t;
if(fac[a][i] > SV)
t = (b % fac[a][i]) ? 1 : fac[a][i];
else
t = G[fac[a][i]][b % fac[a][i]];
b /= t,res *= t;
}
return res;
}
}F;
