滑蒻稽的博客

扩展欧几里得算法(exgcd)

欧几里得算法

首先我们来回顾一下求解2个数的最大公约数(gcd,Greatest Common Divisor)的欧几里得算法:
这个算法的核心是 g c d ( a , b ) = g c d ( b , a   m o d   b ) gcd(a,b)=gcd(b,a \ mod \ b) gcd(a,b)=gcd(b,a mod b),即 a a a b b b的公约数等于 b b b a   m o d   b a \ mod \ b a mod b的公约数.
至于为什么,oi-wiki上有,我们就不证了,但为了理解这个公式,我们举个例子
g c d ( 3 , 5 ) = g c d ( 5 , 3 ) = g c d ( 2 , 3 ) = g c d ( 3 , 2 ) = g c d ( 1 , 2 ) = g c d ( 2 , 1 ) = g c d ( 1 , 1 ) = g c d ( 0 , 1 ) = g c d ( 1 , 0 ) gcd(3,5)=gcd(5,3)=gcd(2,3)=gcd(3,2)=\\gcd(1,2)=gcd(2,1)=gcd(1,1)=gcd(0,1)=gcd(1,0) gcd(3,5)=gcd(5,3)=gcd(2,3)=gcd(3,2)=gcd(1,2)=gcd(2,1)=gcd(1,1)=gcd(0,1)=gcd(1,0)
然后你会发现,我们将要用1去模0了,但这显然是没有意义的*(a模b的定义为a除以b的余数).并且,1和0是没有公约数的(0没有任何约数),那么我们的程序好像就无法继续进行了.我们又观察到,5和3因为互质,最大公约数是1,恰好是当b=0的时候a的值.那么,其他的数进行gcd操作也会有这样的性质吗?

再来个例子
g c d ( 4 , 2 ) = g c d ( 0 , 2 ) = g c d ( 2 , 0 ) gcd(4,2)=gcd(0,2)=gcd(2,0) gcd(4,2)=gcd(0,2)=gcd(2,0)
好像确实是这样的!我们只需输出当b等于0时a的值即可!

实际上,当b=0时,就代表着前一步的a%b==0,也就是 b ∣ a b|a ba,那么我们的最大公约数显然就是b(下一步的a)了.

*如果你编译运行1%0这段代码,编译器报错[Warning] division by zero [-Wdiv-by-zero],程序RE

这样就可以写出代码:

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

在进入下一部分前,我们先来分析一下这个算法的复杂度.细心的读者已经发现了,上面的例子中有些操作是无效的,他们仅仅交换了a和b的位置,但对于程序运行来说,这是必需的.但他们的数量显然小于等于有效操作的数量,对于时间复杂度的分析来说可以略去.

回顾到取模这个运算,在正数意义下它的代码可以写成这样(你是否知道负数的取模运算?可以参考这篇文章):

while(a>b) a-=b//=a%b

实际上,G++中实现的取模不是这么简单.具体来说,若 a = q b + r a=qb+r a=qb+r,则 a   m o d   b = r a \ mod \ b=r a mod b=r(例:7%(-3)=1 -7%(-3)=-1 -7%3=-1).并且,G++的原则是使商尽可能大
有点扯远了,最重要的是知道上面在正数意义下的代码!我们再来举几个例感性理解一下复杂度
不用感性理解了,请读者们自己思考吧,我直接说结论好了,gcd中的有效操作至少能让a减半.

那么,gcd的复杂度就为 l o g   n log \ n log n

补充知识:

  • 多个数的最大公约数/最小公倍数求法:每次取出2个数求得最大公约数/最小公倍数后将最大公约数/最小公倍数放回去
  • 算数基本定理(用于求lcm,Least Common Multiple,最小公倍数): g c d ( a , b ) × l c m ( a , b ) = a × b gcd(a,b) \times lcm(a,b)=a \times b gcd(a,b)×lcm(a,b)=a×b

扩展欧几里德算法(exgcd)

恶心gcd

扩展欧几里得定理(Extended Euclidean algorithm, EXGCD),常用于求 a x + b y = g c d ( a , b ) ax+by=gcd(a,b) ax+by=gcd(a,b)的一组可行解。

P.S 原来extra指额外,extend才是扩展…
定理: x 1 = y 2 , y 1 = x 2 − ⌊ a b ⌋ x_1=y_2,y_1=x_2-\lfloor \frac a b \rfloor x1=y2,y1=x2ba
证明(摘自oi-wiki,但手写一遍才能更好的弄懂):


a x 1 + b y 1 = gcd ⁡ ( a , b ) ax_1+by_1=\gcd(a,b) ax1+by1=gcd(a,b)
b x 2 + ( a   m o d   b ) y 2 = gcd ⁡ ( b , a   m o d   b ) bx_2+(a\bmod b)y_2=\gcd(b,a\bmod b) bx2+(amodb)y2=gcd(b,amodb)
由欧几里得定理可知: gcd ⁡ ( a , b ) = gcd ⁡ ( b , a   m o d   b ) \gcd(a,b)=\gcd(b,a\bmod b) gcd(a,b)=gcd(b,amodb)
所以 a x 1 + b y 1 = b x 2 + ( a   m o d   b ) y 2 ax_1+by_1=bx_2+(a\bmod b)y_2 ax1+by1=bx2+(amodb)y2
又因为 a   m o d   b = a − ( ⌊ a b ⌋ × b ) a\bmod b=a-(\lfloor\frac{a}{b}\rfloor\times b) amodb=a(ba×b)
所以 a x 1 + b y 1 = b x 2 + ( a − ( ⌊ a b ⌋ × b ) ) y 2 ax_1+by_1=bx_2+(a-(\lfloor\frac{a}{b}\rfloor\times b))y_2 ax1+by1=bx2+(a(ba×b))y2
a x 1 + b y 1 = a y 2 + b x 2 − ⌊ a b ⌋ × b y 2 = a y 2 + b ( x 2 − ⌊ a b ⌋ y 2 ) ax_1+by_1=ay_2+bx_2-\lfloor\frac{a}{b}\rfloor\times by_2=ay_2+b(x_2-\lfloor\frac{a}{b}\rfloor y_2) ax1+by1=ay2+bx2ba×by2=ay2+b(x2bay2)
因为 a = a , b = b a=a,b=b a=a,b=b ,所以 x 1 = y 2 , y 1 = x 2 − ⌊ a b ⌋ y 2 x_1=y_2,y_1=x_2-\lfloor\frac{a}{b}\rfloor y_2 x1=y2,y1=x2bay2

证明中比较重要的就是这个式子(其实就是取模运算的数学表达式):
a   m o d   b = a − ⌊ a b ⌋ × b a \ mod \ b=a-\lfloor \frac a b \rfloor \times b a mod b=aba×b
知道了 x x x y y y的递归表达式,我们就可以进行程序设计了.

既然是递归,就一定有最底层的 x x x y y y能够让我们轻易求出值,回到基础gcd算法的递归结束位置, a n o w = g c d ( a , b ) , b n o w = 0 a_{now}=gcd(a,b),b_{now}=0 anow=gcd(a,b),bnow=0,带入 a x + b y = g c d ( a , b ) ax+by=gcd(a,b) ax+by=gcd(a,b),得 g c d ( a , b ) × x + 0 × y = g c d ( a , b ) {gcd(a,b)} \times x +0 \times y=gcd(a,b) gcd(a,b)×x+0×y=gcd(a,b)

那么显然 x x x取1, y y y取任意值即可(我们选择取0).

考虑到在源代码上修改,我们仍然返回gcd的值,在过程中计算x和y即可.

代码:

int exgcd(int a,int b,int &x,int &y)
{
	if(b==0)
	{
		x=1,y=0;
		return a;
	}
	int ret=exgcd(b,a%b,x,y);
	int temp=x;
	x=y;
	y=temp-(a/b)*y;
	return ret; 
}
posted @ 2021-01-31 20:33  huaruoji  阅读(94)  评论(0编辑  收藏  举报