chap 1-4 数论算法-欧几里得算法

1. 辗转相除法(Euclid’s algorithm)

有了模运算，我们可以利用它解决不少关于数字的问题。其中，最为著名的问题当属求解两个数的最大公因数(gcd, greatest common divisor)。最容易想到的方法就是，先对两数分别做质因数分解（factoring），然后将它们相同的质因子相乘得到结果。如：$1035=3^3\cdot 5\cdot 23,759=3\cdot 11\cdot 23,gcd(1035,759)=3\cdot 23=69$。不幸的是：

目前还没有足够高效的分解质因数的算法。

当数字非常大时，算法在有生之年内都不会运行完。有无更好的求解gcd的算法？2000年前其实就有了，叫做辗转相除法/欧几里得算法(Euclid’s algorithm)：

int euclid(int a, int b) {
	if (b == 0) return a;
	return euclid(b, a % b);
}

测试结果：

euclid_test
Input the first integer a: 18
Input the second integer b: 12
gcd of a and b is: 6
time = 2ms
======================================
euclid_test
Input the first integer a: 56
Input the second integer b: 128
gcd of a and b is: 8
time = 1ms
======================================
euclid_test
Input the first integer a: 37
Input the second integer b: 91
gcd of a and b is: 1
time = 1ms
======================================

证明上述代码的正确性，就是证明：

若 $a,b$都是正整数，且 $a\geq b$，则 $gcd(a,b)=gcd(b,a\mod b)$。（欧几里得法则）

$Proof$: 设$d=gcd(a,b)$，则$a=k_1d,b=k_2d$，其中$k_1,k_2$互质。$a-b=(k_1-k_2)d$，$k_2$必定与$k_1-k_2$互质，易得$gcd(b,a-b)=gcd(a,b)$。重复这一步骤，就得到了：$gcd(a,b)=gcd(b,a\mod b)$。假如$a<b$，那么进行一次递归后，$a,b$将会互换位置，所以不需要比较大小。▮

下面分析时间复杂度，就是要搞清楚每次$a\mod b$到底使$a$减小了多少。可以证明出：

若 $a\geq b$，则 $a\mod b < \frac{a}{2}$。

$Proof$: 若$b\leq\frac{a}{2}$，则$a\mod b < b\leq\frac{a}{2}$；若$b>\frac{a}{2}$，则$a\mod b=a-b<\frac{a}{2}$。▮

由此可以推断：任意两次递归调用，$a,b$的值就会变为原来的一半（减小一位），假设初始的$a,b$都是$n$ bits，那么该算法需要进行$2n$次递归调用。每次调用进行取模操作（即$n-bit$除法），故总复杂度为$O(n^3)$。

2. 扩展欧几里得算法

给定两个正整数$a,b$，我们可以运行Euclid’s algorithtm去计算它们的gcd。现在反过来，给定一个正整数$d$，如何验证$d$是否等于$gcd(a,b)$？下面的定理为我们提供了简便的方法：

若 $d$整除 $a$和 $b$，且存在整数 $x$和 $y$使 $d=ax+by$，则 $d=gcd(a,b)$。

$Proof$: $d$整除 $a$和 $b$，可以推出$d\leq gcd(a,b)$；$d=ax+by$，可以推出$gcd(a,b)$整除$d$，表明$d\geq gcd(a,b)$；综上，$d=gcd(a,b)$。▮

普通的Euclid algorithm只能求出$d$，有没有方法同时把$x$和$y$也给求出来？扩展欧几里得算法可以做到：

function extended_euclid(a, b){
    if b = 0: return (1, 0, a)
    (x', y', d) = extended_euclid(b, a mod b)
    return (y', x'-(a/b)y', d)
}

$Proof$:

乍一看有点乱，我们一步步分析。总体框架和普通欧几里得算法一致，把(a, b)变为(b, a mod b)；递归到最后一层时，返回a，也就是最终的d。到这里为止，与普通的欧几里得算法一模一样。

在最后一层递归时，易得$d=a\cdot 1 + b\cdot 0$，显然应该返回$x’=1,y’=0$。现在应该分析的就是，深层递归返回的$(x’,y’):gcd(b, a\mod b)=bx’+(a\mod b)y’$，与本层递归应该返回的$(x,y):gcd(a,b)=ax+by$，之间有何关系？下面是推导：

$$
d=gcd(a,b)=gcd(b,a\mod b)=bx’+(a\mod b)y’=bx’+(a-\lfloor a/b\rfloor b)y’=ay’+b(x’-\lfloor a/b\rfloor y’)
$$

$$
x=y’,y=x’-\lfloor a/b\rfloor y’
$$

综上，上述算法是正确的。▮

举例说明上述算法的步骤，计算$gcd(48,27)$，先计算出$d$：

$$
gcd(48,27):\underline{48}=1\cdot\underline{27}+21
$$

$$
gcd(27,21):\underline{27}=1\cdot\underline{21}+6
$$

$$
gcd(21,6):\underline{21}=3\cdot\underline{6}+3
$$

$$
gcd(6,3):\underline{6}=2\cdot\underline{3}+0
$$

$$
gcd(3,0):d=3
$$

下面开始自下而上计算$x,y$：

$$
gcd(3,0):d=3=1\cdot\underline{3}+0,(x,y)=(1,0)
$$

$$
gcd(6,3):0=\underline{6}-2\cdot\underline{3}=>d=3=1\cdot\underline{6}-1\cdot\underline{3},(x,y)=(1,-1)
$$

$$
gcd(21,6):3=\underline{21}-3\cdot\underline{6}=>d=3=-1\cdot\underline{21}+4\cdot\underline{6},(x,y)=(-1,4)
$$

$$
gcd(27,21):6=\underline{27}-1\cdot\underline{21}=>d=3=4\cdot\underline{27}-5\cdot\underline{21},(x,y)=(4,-5)
$$

$$
gcd(48,27):21=\underline{48}-1\cdot\underline{27}=>d=3=-5\cdot\underline{48}+9\cdot\underline{27},(x,y)=(-5,9)
$$

最终得出$(x,y,d)=(-5,9,3)$。

该算法的复杂度与普通欧几里得算法一致，都是$O(n^3)$。

具体的代码实现：

void extended_euclid(int& x, int& y, int& d, int a, int b) {
	if (b == 0) {
		x = 1; y = 0; d = a;
		return;
	}
	int xx; int yy;
	extended_euclid(xx, yy, d, b, a % b);
	x = yy;
	y = xx - (a / b) * yy;
}

测试结果：

extended_euclid_test
Input the first integer a: 18
Input the second integer b: 12
(x, y, d): (1, -1, 6)
time = 1ms
======================================
extended_euclid_test
Input the first integer a: 56
Input the second integer b: 128
(x, y, d): (7, -3, 8)
time = 1ms
======================================
extended_euclid_test
Input the first integer a: 37
Input the second integer b: 91
(x, y, d): (32, -13, 1)
time = 1ms
======================================

3. 除法的模运算

铺垫了这么久，终于讲到了除法的模运算。如何进行除法的模运算？最容易想到的方法，除法就是乘法的逆运算：

$$
5\cdot 6 \mod 7=2=>2/6 \mod 7=5
$$

但是，这样做有可能会产生多个结果，如：

$$
10/5\mod 15=2,5,8\cdots
$$

这种情况下，我们就认为没有结果，类似于整数除法的除零操作。为了保证答案只有一个，类比倒数，我们不如直接定义：

$$
x/y\equiv xy^{-1} (\mod N)
$$

找到满足要求的$y^{-1}$就可以保证解唯一。在等式两边同乘$y$（$y$和$N$互质）：

$$
xyy^{-1}\equiv x (\mod N)
$$

可以推出：

$$
yy^{-1}\equiv 1 (\mod N)
$$

观察上面的推论，求解除法的模运算关键是要找到那个$y^{-1}$。我们定义：

若 $ax\equiv 1(\mod N)$，则称 $x$为 $a\mod N$的乘法逆元（multiplicative inverse of $a$ modulo $N$）。

我们比较关注的第一个问题是：给定$a,N$，$a\mod N$的乘法逆元是否一定存在？举例：$a=2,N=6$时，就不存在。先给出结论：

$gcd(a,N)=1$时（互质），$a\mod N$的乘法逆元存在；$gcd(a,N)>1$时，则不存在。

$Proof$: $ax\mod N=ax+kN$，则$gcd(a,N)$整除$ax\mod N$。若$gcd(a,N)>1$，$ax\not\equiv 1 (\mod N)$。▮

因此，只要$gcd(a,N)=1$，我们就可以求出$a\mod N$的乘法逆元，如何求呢？

回想一下扩展欧几里得算法，若把$a,N$作为输入，会求出$ax+Ny=d=1$。两边同时取模，得到$ax\equiv 1 (\mod N)$。令人惊讶的是，$x$就是$a\mod N$的乘法逆元！（当然还有其他方法求解）

不过别急，这里求出的$x$是$a\mod N$唯一的乘法逆元吗？实际上：

若$a\mod N$的乘法逆元存在，则它是唯一的。

$Proof$: 反证法。假设$x_1\not\equiv x_2$，$ax_1\equiv 1(\mod N)$且$ax_2\equiv 1(\mod N)$，则$x_1\equiv x_1\cdot 1\equiv x_1\cdot ax_2 \equiv ax_1\cdot x_2\equiv 1\cdot x_2\equiv x_2 (\mod N)$，与$x_1\not\equiv x_2(\mod N)$矛盾。▮

综上所述，求解除法的模运算$a/b(\mod N)$的步骤就是运行扩展欧几里得算法$(x,y,d)=extended-euclid(b,N)$，若$d>1$，结果不存在；若$d=1$，输出$ax \mod N$。复杂度为$O(n^3)$。

具体实现如下：

// by extended_euclid
int get_multiplicative_inverse(int a, int N) {
	int x; int y; int d;
	extended_euclid(x, y, d, a, N);
	if (d != 1) return 0; // inverse not exist
	return x;
}

int mod_divide(int a, int b, int N) {
	int inverse_b = get_multiplicative_inverse(b, N);
    // 0 if inverse not exist, wrap around to 0~N-1
	return inverse_b < 0 ? (a * inverse_b) % N + N : (a * inverse_b) % N; 
}

测试结果：

mod_divide_test (N = 1337)
Input the first integer a: 5
Input the second integer b: 2
Output: 671
time = 0ms
======================================
mod_divide_test (N = 1337)
Input the first integer a: 29
Input the second integer b: 9
Output: 746
time = 1ms
======================================
mod_divide_test (N = 1337)
Input the first integer a: 3873587
Input the second integer b: 45
Output: 482
time = 0ms
======================================

Quick Link: 算法笔记整理