1.欧几里得算法
欧几里德算法又称辗转相除法,用于计算两个整数a,b的最大公约数 gcd(a,b)。基本算法:设 a = qb + r,其中a,b,q,r都是整数,则 gcd(a,b) = gcd(b,r),即 gcd(a,b) = gcd(b,a%b)。
证明:
a = qb + r
如果 r = 0,那么 a 是 b 的倍数,此时显然 b 是 a 和 b 的最大公约数。
如果 r ≠ 0,任何整除 a 和 b 的数必定整除 a - qb = r,而且任何同时整除 b 和 r 的数必定整除 qb + r = a,所以 a 和 b 的公约数集合与 b 和r 的公约数集合是相同的。特别的,a 和 b 的最大公约数是相同的。
递归实现:
int gcd(int a, int b)
{
return b == 0 ? a : gcd(b, a%b);
}
2.贝祖等式
在数论中,裴蜀等式(英语:Bézout’s identity)或贝祖定理(Bézout’s lemma)是一个关于最大公约数(或最大公约式)的定理。裴蜀定理得名于法国数学家艾蒂安·裴蜀,说明了对任何整数a、b和它们的最大公约数d,关于未知数x和y的线性丢番图方程(称为裴蜀等式):
ax + by = m 有整数解时当且仅当m是d的倍数。
裴蜀等式有解时必然有无穷多个整数解,每组解x、y都称为裴蜀数,可用扩展欧几里得算法(Extended Euclidean algorithm)求得。
例如,12和42的最大公因数是6,则方程12x+42y=6有解。事实上有(-3)×12 + 1×42 = 6及4×12 + (-1)×42 = 6。
特别来说,方程 ax+by=1 有整数解当且仅当整数a和b互素。
裴蜀等式也可以用来给最大公约数定义: d其实就是最小的可以写成ax+by形式的正整数。(这个定义的本质是整环中“理想”的概念。因此对于多项式整环也有相应的裴蜀定理。)
证明:
如果 a
和 b 中有一个是 0,比如 a=0,那么它们两个的最大公约数是 b。这时裴蜀等式变成 by=m,它有整数解 (x,y) 当且仅当 m 是 b 的倍数,而且有解时必然有无穷多个解,因为 x 可以是任何整数。定理成立。以下设 a和 b 都不为0。
设 A={xa+yb;(x;y)∈Z2} ,下面证明A中的最小正元素是 a 与 b 的最大公约数。
首先, A∩N∗ 不是空集(至少包含 |a| 和 |b|),因此由于自然数集合是良序的, A 中存在 最小正元素 d0=x0a+y0b。考虑 A中任意一个 正元素 p(=x1a+y1b)对 d0 的带余除法:设 p=qd0+r,其中 q 为正整数, 0≤r<d0。但是
r=p−qd0=x1a+y1b−q(x0a+y0b)∈A
而d0 已经是集合 A 中最小的正元素了,又 0≤r<d0,所以 r=0。
因此 d0 | p。也就是说,A中任意一个正元素p都是 d0 的倍数,特别地: d0 | a、d0 | b。因此 d0 是 a 和 b 的公约数。
另一方面,对 a 和 b 的任意正公约数 d,设 a=kd、b=ld,那么
d0=x0a+y0b=(x0k+y0l)d
x0k+y0l≥1 ,因此 d | d0。所以 d0 是 a 和 b 的最大公约数。
在方程 ax+by=m中,如果 m=m0d0,那么方程显然有无穷多个解:
相反的,如果 ax+by=m有整数解,那么 |m|∈A,于是由前可知 d0 | |m|(即 d0 | m)。
m=1时,方程有解当且仅当a、b互质。方程有解时,解的集合是
其中 (x0,y0) 是方程ax+by=d的一个解,可由辗转相除法得到。
所有解中,有且仅有一个解(x,y) 满足 −b≤x≤b,−a≤y≤a
。
3.扩展欧几里得算法
扩展欧几里得算法就是在求 a,b
的最大公约数 m=gcd(a,b) 的同时,求出贝祖等式 ax + by = m的一个解 (x,y)。
扩展欧几里得算法递归实现
有两个数 a,b
,对它们进行辗转相除法,可得它们的最大公约数——这是众所周知的。然后,收集辗转相除法中产生的式子,倒回去,可以得到 ax+by=gcd(a,b)的整数解。
先来看下这个几乎所有总结扩展欧几里得算法的帖子中都会用到的例子
(可能出自wikipedia,毕竟wikipedia上也是用的这个栗子):
用类似辗转相除法,求二元一次不定方程47x+30y=1的整数解。
47=30*1+17
30=17*1+13
17=13*1+4
13=4*3+1
然后把它们改写成“余数等于”的形式
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。
这些式子自己手写一遍才能理解的更清楚。
这个“倒回去”的过程对应的就是代码中递归向上返回的过程。现在的问题就是要找出这个可以倒回去的递推关系。
参照ACM之家中的证明:
设:a>b。
推理1,显然当 b=0,gcd(a,b)=a。此时 x=1,y=0;//推理1
推理2,ab!=0 时
设 ax1+by1=gcd(a,b);
bx2+(a mod b)y2=gcd(b,a mod b);
根据朴素的欧几里德原理有 gcd(a,b)=gcd(b,a mod b);
则:ax1+by1=bx2+(a mod b)y2;
即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+bx2-(a/b)*by2;
根据恒等定理得:x1=y2 ,y1=x2-(a/b)*y2;//推理2
这样我们就得到了求解 x1,y1 的方法:x1,y1 的值基于 x2,y2.
上面的思想是以递归定义的,因为 gcd 不断的递归求解一定会有个时候 b=0,所以递归可以结束。
这样我们就找到了递推关系:
x1=y2;
y1=x2−(a/b)∗y2;
递推的终止条件再上面也已经给出:
当b=0,gcd(a,b)=a。此时x=1,y=0;
由此我们可以得到递归的扩展欧几里得算法的代码:
int exgcd(int a, int b, int &x, int &y)
{
if(b == 0)
{//推理1,终止条件
x = 1;
y = 0;
return a;
}
int r = exgcd(b, a%b, x, y);
//先得到更底层的x2,y2,再根据计算好的x2,y2计算x1,y1。
//推理2,递推关系
int t = y;
y = x - (a/b) * y;
x = t;
return r;
}
现在我们知道了 a 和 b 的最大公约数是 gcd ,那么,我们一定能够找到这样的 x 和 y ,使得: a*x + b*y = gcd 这是一个不定方程(其实是一种丢番图方程),有多解是一定的,但是只要我们找到一组特殊的解 x0 和 y0 那么,我们就可以用 x0 和 y0 表示出整个不定方程的通解:
x = x0 + (b/gcd)*t
y = y0 – (a/gcd)*t
扩展欧几里得算法的非递归实现!
算法E
(推广的欧几里得算法). 给定两个正整数m和n, 计算他们的最大公因数d,并计算两个未必为正数的整数 a 和 b, 使得 am+bn=d.E1. [初始化.] 置 a′←b←1,a←b′←0,c←m,d←n.
E2. [除法.] 令 q 和 r 分别是用 d 除 c 所得的商和余数. (我们有 c=qd+r 和 0≤r<d.)
E3. [余数为0?] 如果 r=0 ,算法终止,此时有 am+bn=d , 正如所求.
E4. [循环.] 置 c←d,d←r,t←a′,a′←a,a←t−qa,t←b′,b′←b,b←t−qb , 然后返回 E2.
┃
假定 m = 1769 , n = 551, 我们相继得到:
1 | 0 | 0 | 1 | 1769 | 551 | 3 | 116 |
0 | 1 | 1 | -3 | 551 | 116 | 4 | 87 |
1 | -4 | -3 | 13 | 116 | 87 | 1 | 29 |
-4 | 5 | 13 | -16 | 87 | 29 | 3 | 0 |
答案是正确的: 5×1769 - 16×551 = 8845 - 8816 = 29, 这是1769和551的最大公因数。
《计算机程序设计艺术》用的数学归纳法证明。
首先我们发现 等式:
a′m+b′n=c,am+bn=d
在步骤 E2
执行后总是成立。 (1)
证明(1):
算法首次转到E2时,由上表,显然等式成立。
步骤 E2,E3 并不会改变等式的正确性。
所以只需要证明步骤 E4 不改变他们的正确性就可以了。
执行 E4 之前,
a'm + b'n = c (*)
am + bn = d (**)
成立。
a'm + b'n = c = qd + r;
a'm + b'n - q(am + bn) = r;
(a' - qa)m + (b' - qb)n = r (***)
执行 E4 之后,c←d, d←r, t←a', a'←a, a←t-qa, t←b', b'←b, b←t-qb;
所以
(**)变为 a' + b' = c
(***)变为 am + bn = d
所以步骤 E4 不改变等式的正确性,得证。
现在我们对n归纳,证明算法E
是正确的:如果 m 是 n 的倍数,算法显然能够正确执行,因为在第一次达到 E3 的时候就立即完成了。 当 n=1 时这种情况总是出现,剩下的情况仅是 n>1 且 m 不是 n 的倍数。在这种情况下,算法在第一次执行后置 c←n 和 d←r,并且由于 r←n, 按归纳法可以假设 , d 的最终值是 n 和 r 的最大公因数。根据上文欧几里得算法的论证,数偶 {m,n} 和 {n,r} 具有相同的公因数,特别的,它们具有相同的最大公因数。因此 d 是 m 和 n 的最大公因数, 而且由(1)式, am+bn=d。看完了优美的数学归纳法证明,我们来找找 a′m+b′n=c,am+bn=d
的递推公式:
首先我们发现a′m+b′n=c,am+bn=d
这两个式子是等价的,只是同一个式子的不同阶段。a′m+b′n=c 保存了 am+bn=d 的上一个阶段。所以我们只需要得到其中一个的递推公式就可以了。
实际上 a′ 只是保存了 a 之前的一个阶段。
下面我们求解 am+bn=d
的递推公式。
如果给上表标上序号:
余数等于 | |||||||
---|---|---|---|---|---|---|---|
1 | 0 | 1 | 1769 | 551 | 3 | 116 | 116 = m + (-3)n |
2 | 1 | -3 | 551 | 116 | 4 | 87 | 87 = -4m + 13n |
3 | -4 | 13 | 116 | 87 | 1 | 29 | 29 = 5m +(-16)n |
4 | 5 | -16 | 87 | 29 | 3 | 0 |
第 i
再 gcd() 的过程中,我们会得到一系列的余数,但我们不仅仅想得到这些余数的序列,我们还想得到 aim+bin=di 这样的等式。
因为 di=ri−1,所以 d 序列本身就是 r
的序列。从上表可以清楚的看到这二者之间的关系。
第一行:a1=0,b1=1
0×1769 + 1×551 = 551 初始条件是成立的,即 a1m+b1n=d1;
第二行:由第一行的 1769 = 3×551 + 116 ,把余数116表示为 116 = 1×1769 + (-3)×551 = m + (-3)n, 我们注意到第一行的余数 r1
等于第二行的 d2 所以 a2=1,b2=−3, a2m+b2n=r1=d2第三行:由第二行的 551 = 4×116 +87, 把余数87表示为 87 = 551 + (-4)×116, 其中551 = n, 而116已经求得为m+(-3)n, 所以 87 = n + (-4)(m+(-3)n) = -4m + 13n ,即 a3=−4,b3=13
, a3m+b3n=r2=d3;
第四行:由第三行的 116 = 1×87 + 29,把余数29表示为 29 = 116 + (-1)×87 = (m + (-3)n) + (-1)(-4m + 13n) = 5m +(-16)n,所以a4=5,b4=−16
, a4m+b4n=r3=d4;第五行:由第4行的 87 = 3×29 , r4=0,算法终止, d4即为最大公约数, a4,b4即为求得的满足 am+bn=gcd(m,n) 的 a,b
。
我们发现 :
di=ri−1=ci−1−qi−1di−1=di−2−qi−1di−1=(ai−2m+bi−2n)−qi−1(ai−1m+bi−1n)
=(ai−2−qi−1ai−1)m+(bi−2−qi−1bi−1)n
又 di=aim+bin
所以
ai=ai−2−qi−1ai−1;
bi=bi−2−qi−1bi−1;
实际上,这里的 ai−2,bi−2 就是 a′i,b′i。
于是,我们就找到了 am+bn=d 的系数 a 和 b
的递推公式。
最后,给出非递归的C++代码:
int exgcd(int m, int n, int &x, int &y) {
if (n == 0) {
x = 1; y = 0;
return m;
}
int a, a1, b, b1, c, d, q, r, t;
a1 = b = 1;
a = b1 = 0;
c = m; d = n;
q = c/d; r = c%d;
while (r) {
c = d;
d = r;
t = a1;
a1 = a;
a = t - q * a;
t = b1;
b1 = b;
b = t - q * b;
q = c/d;
r = c%d;
}
x = a; y = b;
return d;
}