登录  
 加关注
查看详情
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

银河里的星星

落在人间

 
 
 

日志

 
 

轻松理解扩展欧几里德算法  

2009-09-02 22:01:28|  分类: 算法与acm |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

转载请注明作者:phylips@bmy 出处:http://duanple.blog.163.com/blog/static/7097176720098210128972/

欧几里德算法

用来求两个数的最大公约数的算法。具体如下
int gcd(int a,int b){
   if(b == 0) return a;
   return gcd(b,a%b); 
}

首先看正确性证明,实际上需要证明gcd(a,b)=gcd(b,a%b),我们只要证明gcd(a,b)=gcd(a-b,b)即可,因为可以由此逐步扩展为gcd(a,b) = gcd(a-k*b,b),而 gcd(a-k*b,b)=gcd(a%b,b)。
因为a,b的公约数必然是a-b,b的公约数故 gcd(a,b) <= gcd(a-b,b);另a-b b的公约数也必然是a b的公约数,gcd(a,b) >= gcd(a-b,b).所以gcd(a,b) = gcd(a-b,b)。证明完毕。

再看复杂度,我们可以看到,如果a >= b 那么 a%b < a/2(可以从b>a/2 <= a/2两种情况考虑),这样我们可以看到经过两次调用gcd(a,b)=gcd(b,a%b)=gcd(a%b,b%(a%b)),原来的a就变成了a%b,也就是说a的规模降低了1/2,通过这点可以观察出是log级的复杂度。

还有一点值得思考。我们可以看到递归时采用了这样一个手段gcd(a,b)=gcd(b,a%b),递归中把a b交换了。为什么这样呢,首先考虑gcd(a,b) = gcd(a%b,b),这个递归一次以后就终止了无法保证a b可以继续减小。当然实际上这两个等式都是成立的,但是gcd(b,a%b)才能保证求解的顺利进行。当我们碰到类似的情况时,就可以考虑能否互换一下顺序就可以继续递归下去呢。

另有一点,我们可以看到在上面的过程中,实际上存在一种假设或者说是不变性:a >= b。实际上gcd(b,a%b)就可以完成这种保证,因为b > a%b 。看到这种性质,我们就可以直接判断b==0就可以结束程序,因为肯定是b首先=0。这个不变性只有在第一次调用时可能不遵守,比如我们调用gcd(3,5),但是我们可以发现它并不影响这个程序的正确性,因为经过一次调用后gcd(3,5)=gcd(5,3)自动完成了调整。

 

扩展欧几里德算法


欧几里德算法比较简单直接,但是扩展形式就不是那么一目了然了。以前也存在一个错误的理解,一直认为扩展欧几里德是用来求解同余方程的,而ax+by = c是先转换为同余方程之后再求解的,实际上这个认识是刚好相反的。正是这样的一个理解,导致以前看到ax+by=c都是想转换为同余方程。最近在水木上也碰到了一个类似的问题ax-by=c,实际上这个只需要看成ax+b*(-y)就可以转换成ax+by=c的形式了,求解出来把y反号既可以了。

而扩展欧几里德算法实际上正是为了求解ax+by=gcd(a,b)这样的方程创造的。

思路是这样的:利用递归。

设d= gcd(a,b) = gcd(b,a%b),结合欧几里德算法,我们首先可以把gcd(a,b)=ax+by转换为另一个方程gcd(b,a%b) = bx'+a%by'.为什么进行这种转换呢?因为通过欧几里德算法我们知道这种转换一定会停止,并且最终变成gcd(a,b)=ax+by会变成gcd(a,0)=ax + 0*y而这个方程很容易解,只要x =1 y = 0就可以满足了(一个猜想:y可以是任何数,不过选择不同的y对应最后解集合中不同的一组解)。

但是进行这样的转换后,很明显x' y'就不是原来的解了。根据递归的思想,假设现在我们已经求出了x' y',剩下的关键就是如何用x' y'求出x y.我们观察gcd(b,a%b) = bx'+a%by',因为实际上gcd(b,a%b)=gcd(a,b)只要把右边重新写成 =ax+by就可以了所以需要对bx'+a%by'进行变形,因为a%b = a-ceil(a/b)b,故bx'+a%by' = bx'+(a-ceil(a/b)b)y' = ay' + b(x'-ceil(a/b)y') ,这样便可得出 x = y' y = x'-ceil(a/b)y'。

于是递归表达式可以写出了。

int(x,y,d) ext_euclid(a,b)
if b == 0 return (1,0,a)
(x',y',d) = ext_euclid(b,a%b)
return (y',x'-ceil(a/b)y',d)

这样我们就得到了ax+by=gcd(a,b)的一组解。ax0 + by0 = gcd(a,b),实际上很容易得出一系列解
x = x0+b*i
y = y0-a*i
带入可以观察到 a(x0+b*i) + b(y0-a*i) = ax0+by0 = gcd(a,b)。

相关应用

从上面的过程可以看到,ax+by=gcd(a,b)一定有解,再看一般形式ax+by=c,只有当且仅当c是gcd(a,b)的整数倍时才有解,否则无解。因为如果我们两边同除gcd(a,b),会出现这样一种情况a'x+b'y=c/gcd(a,b),如果c不是gcd(a,b)的整数倍,那么左边是整数,右边是分数,显然无解。

ax+by=c的求解可以先求出ax+by=gcd(a,b),然后将x y扩大c/gcd(a,b)倍就可以了。

再看ax+by=c,实际上它等价于同余方程,ax =c mod b,这样同余方程就可以用扩展欧几里德算法求解。

另外如果a b互质那么ax+by=gcd(a,b) = 1,转换为同余方程ax = 1 mod b这样x就变成了a的关于mod b的逆元。也就是说求逆元也不过是扩展欧几里德算法的一个特殊情况。而这个操作将是RSA公钥系统的一个基本操作。下面再具体介绍整个RSA系统的理论基础。

  评论这张
 
阅读(2359)| 评论(1)

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018