欧几里得算法和扩展欧几里得算法

欧几里得算法

欧几里德算法又称辗转相除法，用于计算两个整数a,b的最大公约数。基本算法：设a=qb+r，其中a，b，q，r都是整数，则gcd(a,b)=gcd(b,r)，即gcd(a,b)=gcd(b,a%b)。

基本代码实现：

1	int gcd( int a, int b)

2

{

3

if (b==0)

4

return a;

5

return

6	gcd(b,a%b);

7

}

扩展欧几里得算法

已知整数a、b，扩展欧几里得算法可以在求得a、b的最大公约数的同时，能找到整数x、y（其中一个很可能是负数），使它们满足贝祖等式

。有两个数a,b，对它们进行辗转相除法，可得它们的最大公约数——这是众所周知的。然后，收集辗转相除法中产生的式子，倒回去，可以得到ax+by=gcd(a,b)的整数解。

用类似辗转相除法，求二元一次不定方程47x+30y=1的整数解。

a，b为给定的值，通过gcd求解的方法求47与30的最大公约数（即a，b的最大公约数）

47=30*1+17
30=17*1+13
17=13*1+4
13=4*3+1
gcd（a，b）=gcd（b，a%b）

然后把它们改写成“余数等于”的形式

17=47*1+30*(-1) //式1
13=30*1+17*(-1) //式2
4=17*1+13*(-1) //式3
1=13*1+4*(-3)

然后把它们“倒回去”

1=13*1+4*(-3) //应用式3
1=13*1+[17*1+13*(-1)]*(-3)
1=13*4+17*(-3) //应用式2
1=[30*1+17*(-1)]*4+17*(-3)
1=30*4+17*(-7) //应用式1
1=30*4+[47*1+30*(-1)]*(-7)
1=30*11+47*(-7)

得解x=-7, y=11。

基本算法：对于不完全为 0 的非负整数 a，b，gcd（a，b）表示 a，b 的最大公约数，必然存在整数对 x，y ，使得 gcd（a，b）=ax+by。

其实大致意思就是a 和 b 的最大公约数是 gcd（a，b） ，那么，我们一定能够找到这样的 x 和 y ，使得: a*x + b*y = gcd 这是一个不定方程（其实是一种丢番图方程），有多解是一定的，但是只要我们找到一组特殊的解 x0 和 y0 那么，我们就可以用 x0 和 y0 表示出整个不定方程的通解：

x = x0 + (b/gcd)*t

y = y0 – (a/gcd)*t

为什么不是：

x = x0 + b*t

y = y0 – a*t

那是因为：

b/gcd 是 b 的因子， a/gcd 是 a 的因子是吧？那么，由于 t的取值范围是整数，你说 (b/gcd)*t 取到的值多还是 b*t 取到的值多？同理，(a/gcd)*t 取到的值多还是 a*gcd 取到的值多？那肯定又要问了，那为什么不是更小的数，非得是 b/gcd 和a/gcd ？

注意到：我们令 B = b/gcd ， A = a、gcd ， 那么，A 和 B 一定是互素的吧？这不就证明了 最小的系数就是 A 和 B 了吗？要是实在还有什么不明白的，看看《基础数论》（哈尔滨工业大学出版社），这本书把关于不定方程的通解讲的很清楚

现在，我们知道了一定存在 x 和 y 使得 ： a*x + b*y = gcd ， 那么，怎么求出这个特解 x 和 y 呢？只需要在欧几里德算法的基础上加点改动就行了。

我们观察到：欧几里德算法停止的状态是： a= gcd ， b = 0 ，那么，这是否能给我们求解 x y 提供一种思路呢？因为，这时候，只要 a = gcd 的系数是 1 ，那么只要 b 的系数是 0 或者其他值（无所谓是多少，反正任何数乘以 0 都等于 0 但是a 的系数一定要是 1），这时，我们就会有： a*1 + b*0 = gcd

当然这是最终状态，但是我们是否可以从最终状态反推到最初的状态呢？

假设当前我们要处理的是求出 a 和 b的最大公约数，并求出 x 和 y 使得 a*x + b*y= gcd ，而我们已经求出了下一个状态：b 和 a%b 的最大公约数，并且求出了一组x1 和y1 使得： b*x1 + (a%b)*y1 = gcd ， 那么这两个相邻的状态之间是否存在一种关系呢？

我们知道： a%b = a - (a/b)*b（这里的 “/” 指的是整除，例如 5/2=2 , 1/3=0），那么，我们可以进一步得到：

gcd = b*x1 + (a-(a/b)*b)*y1

= b*x1 + a*y1 – (a/b)*b*y1

= a*y1 + b*(x1 – a/b*y1)

对比之前我们的状态：求一组 x 和 y 使得：a*x + b*y = gcd ，是否发现了什么？

这里：

x = y1

y = x1 – a/b*y1

01

证明：设 a>b。

02

03	推理1，显然当 b=0，gcd（a，b）=a。此时 x=1，y=0； //推理1

04

05	推理2，ab!=0 时

06

07	设 ax1+by1=gcd(a,b);

08

09	bx2+(a mod b)y2=gcd(b,a mod b);

10

11	根据朴素的欧几里德原理有 gcd(a,b)=gcd(b,a mod b);

12

13	则:ax1+by1=bx2+(a mod b)y2;

14

15	即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+bx2-(a/b)*by2;

16

17	根据恒等定理得：x1=y2; y1=x2-(a/b)*y2; //推理2

18

19	这样我们就得到了求解 x1,y1 的方法：x1，y1 的值基于 x2，y2.

20

21	上面的思想是以递归定义的，因为 gcd 不断的递归求解一定会有个时候 b=0，所以递归可以结束。

扩展欧几里德的递归代码：

01	#include <iostream>

02	using namespace std;

03

04	int exgcd( int a, int b, int & x, int & y){

05	if (b == 0){

06	//根据上面的推理1，基本情况

07

x = 1;

08

y = 0;

09

return a;

10

}

11	int r = exgcd(b, a%b, x, y);

12

//根据推理2

13	int t = y;

14	y = x - (a/b)*y;

15

x = t;

16

return r;

17

}

18

19	int main() {

20

int x,y;

21	exgcd(47,30,x,y);

22	cout << "47x+30y=1 的一个整数解为: x=" << x << ", y=" << y << endl;

23

return 0;

24

}

非递归实现，比上面的看上去要复杂了不少，不熟悉的话直接用上面的就可以：

01	int exgcd( int m, int n, int &x, int &y)

02

{

03	int x1,y1,x0,y0;

04	x0=1; y0=0;

05	x1=0; y1=1;

06

x=0; y=1;

07	int r=m%n;

08	int q=(m-r)/n;

09

while (r)

10

{

11	x=x0-q*x1; y=y0-q*y1;

12	x0=x1; y0=y1;

13	x1=x; y1=y;

14	m=n; n=r; r=m%n;

15	q=(m-r)/n;

16

}

17

return n;

18

}

扩展欧几里德算法的应用

（1）求解不定方程

用扩展欧几里得算法解不定方程ax+by=c;

这个应该比较好理解了，两个可以同乘以k

1	bool linear_equation( int a, int b, int c, int &x, int &y)

2

{

3	int d=exgcd(a,b,x,y);

4

if (c%d)

5	return false ;

6	int k=c/d;

7	x*=k; y*=k; //求得的只是其中一组解

8	return true ;

9

}

（2）求解模线性方程（线性同余方程）

同余方程 ax≡b (mod n) (也就是 ax % n = b) 对于未知数 x 有解，当且仅当 gcd(a,n) | b (也就是 b % (gcd(a,n))==0 )。且方程有解时，方程有 gcd(a,n) 个解。

求解方程 ax≡b (mod n) 相当于求解方程 ax+ ny= b, (x, y为整数)

1	在方程 3x ≡ 2 (mod 6) 中， d = gcd(3,6) = 3 ，3 不整除 2，因此方程无解。

2

3	在方程 5x ≡ 2 (mod 6) 中， d = gcd(5,6) = 1，1 整除 2，因此方程在{0,1,2,3,4,5} 中恰有一个解: x=4。

证明略去，直接说算法：

首先看一个简单的例子：

5x=4(mod3)

解得x = 2,5,8,11,14…….

由此可以发现一个规律，就是解的间隔是3.

那么这个解的间隔是怎么决定的呢？

如果可以设法找到第一个解，并且求出解之间的间隔，那么就可以求出模的线性方程的解集了.

我们设解之间的间隔为dx.

那么有

a*x = b(mod n);

a*(x+dx) = b(mod n);

两式相减，得到:

a*dx(mod n)= 0;

也就是说a*dx就是a的倍数，同时也是n的倍数，即a*dx是a 和 n的公倍数.为了求出dx,我们应该求出a 和 n的最小公倍数,此时对应的dx是最小的.

设a 和 n的最大公约数为d,那么a 和 n 的最小公倍数为(a*n)/d.

即a*dx = a*n/d;

所以dx = n/d. (d = gcd(a,n) )

因此解之间的间隔就求出来了.

01	bool modular_linear_equation( int a, int b, int n)

02

{

03	int x,y,x0,i;

04	int d=exgcd(a,n,x,y);

05

if (b%d)

06	return false ;

07	x0=x*(b/d)%n; //特解

08	for (i=1;i<d;i++)

09	printf ( "%d\n" ,(x0+i*(n/d))%n);

10	return true ;

11

}

（3）求解模的逆元；

同余方程ax≡b (mod n)，如果 gcd(a,n)== 1，则方程只有唯一解。

在这种情况下，如果 b== 1，同余方程就是 ax=1 (mod n ),gcd(a,n)= 1。

这时称求出的 x 为 a 的对模 n 乘法的逆元。

对于同余方程 ax= 1(mod n )， gcd(a,n)= 1 的求解就是求解方程

ax+ ny= 1，x, y 为整数。这个可用扩展欧几里德算法求出，原同余方程的唯一解就是用扩展欧几里德算法得出的 x 。

扩展欧几里德算法

谁是欧几里德？自己百度去

先介绍什么叫做欧几里德算法

有两个数 a b，现在，我们要求 a b 的最大公约数，怎么求？枚举他们的因子？不现实，当 a b 很大的时候，枚举显得那么的naïve ，那怎么做？

欧几里德有个十分又用的定理： gcd(a, b) = gcd(b , a%b) ，这样，我们就可以在几乎是 log 的时间复杂度里求解出来 a 和 b 的最大公约数了，这就是欧几里德算法，用 C++ 语言描述如下：



由于是用递归写的，所以看起来很简洁，也很好记忆。那么什么是扩展欧几里德呢？

现在我们知道了 a 和 b 的最大公约数是 gcd ，那么，我们一定能够找到这样的 x 和 y ，使得: a*x + b*y = gcd 这是一个不定方程（其实是一种丢番图方程），有多解是一定的，但是只要我们找到一组特殊的解 x0 和 y0 那么，我们就可以用 x0 和 y0 表示出整个不定方程的通解：

x = x0 + (b/gcd)*t

y = y0 – (a/gcd)*t

为什么不是：

x = x0 + b*t

y = y0 – a*t

这个问题也是在今天早上想通的，想通之后忍不住喷了自己一句弱逼。那是因为：

b/gcd 是 b 的因子， a/gcd 是 a 的因子是吧？那么，由于 t的取值范围是整数，你说 (b/gcd)*t 取到的值多还是 b*t 取到的值多？同理，(a/gcd)*t 取到的值多还是 a*gcd 取到的值多？那肯定又要问了，那为什么不是更小的数，非得是 b/gcd 和a/gcd ？

注意到：我们令 B = b/gcd ， A = a、gcd ， 那么，A 和 B 一定是互素的吧？这不就证明了 最小的系数就是 A 和 B 了吗？要是实在还有什么不明白的，看看《基础数论》（哈尔滨工业大学出版社），这本书把关于不定方程的通解讲的很清楚

现在，我们知道了一定存在 x 和 y 使得 ： a*x + b*y = gcd ， 那么，怎么求出这个特解 x 和 y 呢？只需要在欧几里德算法的基础上加点改动就行了。

我们观察到：欧几里德算法停止的状态是： a= gcd ， b = 0 ，那么，这是否能给我们求解 x y 提供一种思路呢？因为，这时候，只要 a = gcd 的系数是 1 ，那么只要 b 的系数是 0 或者其他值（无所谓是多少，反正任何数乘以 0 都等于 0 但是a 的系数一定要是 1），这时，我们就会有： a*1 + b*0 = gcd

当然这是最终状态，但是我们是否可以从最终状态反推到最初的状态呢？

假设当前我们要处理的是求出 a 和 b的最大公约数，并求出 x 和 y 使得 a*x + b*y= gcd ，而我们已经求出了下一个状态：b 和 a%b 的最大公约数，并且求出了一组x1 和y1 使得： b*x1 + (a%b)*y1 = gcd ， 那么这两个相邻的状态之间是否存在一种关系呢？

我们知道： a%b = a - (a/b)*b（这里的 “/” 指的是整除，例如 5/2=2 , 1/3=0），那么，我们可以进一步得到：

gcd = b*x1 + (a-(a/b)*b)*y1

= b*x1 + a*y1 – (a/b)*b*y1

= a*y1 + b*(x1 – a/b*y1)

对比之前我们的状态：求一组 x 和 y 使得：a*x + b*y = gcd ，是否发现了什么？

这里：

x = y1

y = x1 – a/b*y1

以上就是扩展欧几里德算法的全部过程，依然用递归写：



依然很简短，相比欧几里德算法，只是多加了几个语句而已。

这就是理论部分，欧几里德算法部分我们好像只能用来求解最大公约数，但是扩展欧几里德算法就不同了，我们既可以求出最大公约数，还可以顺带求解出使得： a*x + b*y = gcd 的通解 x 和 y