[可能有用科技]更多 GCD 算法

前言#

重写旧文 & 增加新内容。

欧几里得算法#

$$\large \gcd(a,a) = a\\ \gcd(a,b) = \gcd(b,a \bmod b)$$

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;
}

更相减损法#

可半者半之，不可半者，副置分母、子之数，以少减多，更相减损，求其等也。以等数约之。

白话 ：

1. 将被 $2$ 整除的一直除以 $2$ 直到因数不包含 $2$。
2. 用较大的数减去较小的，直到减数和差相等，得到的就是最大公因子。
3. 同时除了 $k$ 次 $2$ 的情况，要在结尾乘以 $2^k$

$$\large \gcd(a,b) = \gcd(\min(a,b),|a - b|)$$

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}$，我们之前得到 ：

$$\large a \le \sqrt[3]{\frac{n}{p}}\\ pa \le p\sqrt[3]{\frac{n}{p}}$$

然后

$$\large p\sqrt[3]{\frac{n}{p}} = \sqrt[3]{p^3\frac{n}{p}} = \sqrt[3]{p^2n}$$

并且 $\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$。

最后回顾一下全过程 ：

1. 预处理 $\sqrt{\mathbb{V}}$ 以内的 $\gcd$，复杂度 $\mathcal{O}(\mathbb{V})$。
2. 预处理 $\mathbb{V}$ 以内的所有数的分解方案，类似线性筛的方式，复杂度 $\mathcal{O}(\mathbb{V})$。
3. 询问的时候需要的信息都预处理出来了，复杂度 $\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;