欧几里德算法（最大公约数算法）完整分析

Posted 2020-09-26

　　欧几里德算法又称为辗转相除法，用于计算两个非负整数的最大公因数。其伪代码如下： 

gcd(a, b) //要求保证传入的a>=b

　　if(b == 0)

　　　　return a

　　return gcd(b, a % b)

　　首先说明这个函数能返回a与b的最大公因数。但是我们不从代码到原理，我们要从原理返回代码。（下面的出现的所有符号均为非负整数）

　　在a与b均非0且a>=b的情况下，若c是a和b的最大公因数(c>0)，那么就有c|a和c|b的同时成立。显然a=i*c，b=j*c，此处应满足1<=j<=i。记i=k*j+r，其中0<=r<j，那么i%j=r，floor(i / j) = k 。而对应的a=i*c=(k*j+r)*c=k*j*c+r*c=k*b+r*c，此处r*c<j*c=b，因此有floor(a / b) = k且a % b = r * c。到此我们发现了一个有趣的现象，a%b与b拥有公因数c。

　　再说明c必定是a%b与b的最大公因数。假设存在d>c，使得d是a%b与b的公因数，那么有a%b=q*d，b=w*d。那么a=e*b+a%b=(e*w+q)*d，显然d也同时是a和b的公因数，这与c是a和b的最大公因数的前提相违背，因此假设不成立，故所有a%b与b的公因数都不会超过c，因此c是a%b与b的最大公因数。

　　到了这里我们可以总结为当a与b均非0且a>=b的情况下，a与b的最大公因数与a%b与b的最大公因数完全一致。

　　是否感觉逐渐明朗起来了。若a与b均非0且a>=b，那么可以通过计算a%b与b的最大公因数得到a与b的最大公因数c。而同样若a%b不为0，则可以计算b%(a%b)与a%b的公因数来变相得到a与b的最大公因数c。一直到某个参数为0为止，注意到任意正整数x与0的最大公因数均为x，此时简单返回x即可。这与我们gcd的代码完全一致。

　　分析完了算法的正确性后，还需要进一步分析这个程序的时间复杂度和空间复杂度。空间复杂度为O(1)应该毋庸质疑，因为在递归的过程中并没有分配额外的空间。下面说明时间复杂度：

　　注意到每一次gcd的调用（不计算递归部分）所花费的时间都是常数，因此调用一次gcd的时间复杂度为O(1)。而在一次最大公因数的计算中gcd又总共被调用了多少次呢？注意到在递归调用的过程中，其传入的参数为b,a%b。其中b<=a，而a%b<=min(a-b,b)<=a/2。因此我们发现每次gcd的调用，一个传入参数值不被改变，而另外一个传入参数值至少减少了一半。由gcd的定义知道当任意一个参数减少到0时递归就将结束，因此在对a，b求最大公因数的过程中，最多调用log2(a)+log2(b)次的gcd，其中log2(a)次令a减半，log2(b)次令b减半。因此总共的时间复杂度为O(1)*(log2(a)+log2(b))=O(log2(max(a,b))。

---

　　欧几里德的算法不仅可以用于得出a与b的最大公因数，还能找到这样一个值i，使得

即i是可以满足条件i*a=j*b的最小自然数，其中j是也是自然数。

　　这该如何做到呢，实际上记c为a和b的最大公因数，则x*a=y*b可以等价地写作x*(a/c)=y*(b/c)，由于c是最大公因数，因此a/c与b/c的最大公因数为1，即二者互质。由算术基本定理可以得出(a/c)|y且(b/c)|x整除。而由于x是正整数，因此x的可能取值为(b/c),2(b/c),3(b/3),...。有趣的是(b/c)*(a/c)=(a/c)*(b/c)，即分别令x=(b/c),y=(a/c)就可以得到公式xa=yb满足，且次数x的取值又正好是可能取值中的最小值，因此i=x=b/c。而正巧欧几里德算法提供了高效计算c的方案，因此i=b/gcd(a,b)。

　　在已知c的情况下，计算i的时空复杂度均为O(1)。

　　事实上这里由所有x的可能取值也可以发现另外一个有趣的性质，凡是满足x*a=y*b的函数，若x为正整数，则x必然是i=(b/c)的倍数。

---

　　再说明欧几里德的一个用途，就是计算xa+yb=c中系数x和y的一个可能取值，其中c是a和b的最大公因数（a与b不能同时为0）。

 　　不妨假设a>=b。若b=0时，则显然c=a，而此时1*a+0*b=c是非常完美的解决方案。

　　现在假设已经找到了q*b+w*(a%b)=c中q和w的取值，那么可以推出c=q*b+w*(a%b)=q*b+w*(a-floor(a/b)*b)=w*a+(q+w*floor(a/b))*b，令x=w，y=(q+w*floor(a/b))，公式xa+yb=c得以满足。

　　而这又与欧几里德算法有什么关系呢？当然有！我们称欧几里德算法中最后一次递归（即出现一个参数为0的情况的那次递归操作）称为情况1，而对应的调用情况1的那次递归称为情况2，调用情况2的称为情况3，这样如是继续下去，则首次调用（即参数为a和b的那次）则必然可以被认为是情况n。由与情况1下已经保证能找到x和y使得等式成立，且对于任意情况k>=1，若情况k下能找到对应的系数使得参数ak,bk满足xk*ak+yk*bk=c，则能找到x_k+1*a_k+1+y_k+1*b_k+1=c，即情况k+1下也能找到使公式成立的对应参数。利用数学归纳法，我们得知对于任意情况i（i为自然数），都能找到对应的系数使得公式成立。因此自然我们也能找到令xa+yb=c的系数x和y。

coe(a, b) //a要保证大于b

　　if(b == 0)

　　　　return (1, 0)

 　　r = coe(b, b % a)

　　return (r.y, r.x + r.y * floor(a / b))

　　coe的时间复杂度与gcd一样，使用同样的分析方法就可以得出。

　　如果希望找到所有这类系数中x值非负且最小的情况，可以结合上一部分的方案。首先找到

其次令x加上或减去（取决于当前x的符号）i，同时另外减去或加上j。此时c=c+0=(x+i)a+(y-j)b=(x-i)a+(y+j)b。这样循环使得x不断靠近0，直到i>x>=0的情况发生。这里要说明这时候所能取得的x确实是最小非负整数。若存在更小的非负整数m和对应的n使得ma+nb=c，那么有(xa+yb)-(ma+nb)=(x-m)*a+(y-n)b=c-c=0，而这里定义的i>x>=x-m，这与i是符合这类公式(ia+jb=0)最小正整数的定义相悖，因此假设不成立，故此时得到的x才是公式(xa+yb=c)的最小非负整数。这里利用乘法和除法可以一步计算出来，因此在已知任意一组x和y的情况下，这个操作的时空复杂度为O(1)。

minX(a, b, x, y)

　　if(x < 0)

　　　　return (x % i + i, y - (floor(x / i) + 1) * j)

　　return ( x % i, y - floor(x / i) * j)