[可能有用科技]更多 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;