exgcd

扩展欧几里得(exgcd):
此算法是用来求方程ax+by=gcd(a,b)的一组可行解的算法(不一定最优最小)。

以下内容摘自洛谷p1082题解第一篇:

我们拿到了一组 a,b。目标是求出满足 ax + by = gcd(a,b)(①) 的整数 x 与 y。

根据普通欧几里得算法,gcd(a,b) = gcd(b, a\mod b)(②)。

如果我们先前已经求出了另一组数 x_2, y_2,它们满足 bx_2 + (a \mod b)y_2 = gcd(b, a \mod b)(③),(提示一下,这个等式与①的结构一致,系数则是根据②替换的。此处的递归形态已经有所显露,别着急)

结合②③,一定有

bx_2 + (a \mod b)y_2 =gcd(b, a \mod b)=gcd(a,b)
在这个“如果”实现的时候,我们的目标变成了求出满足 ax + by = bx_2 + (a \mod b)y_2(④)的 x 与 y。

其中a,b,x_2,y_2都已知,x,y待求。未知数比方程更多,没有唯一解。我们先求出一组必然存在的解,最后将在答案处理时转变为使正整数 x 最小的解。

取模运算是 a \mod b = a - b×(a/b),能理解吧?

等式④实际上是:

ax + by = bx_2 + (a-b×(a/b))y_2

ax + by = bx_2 + ay_2 - b × (a/b)y_2

ax + by = ay_2 + b(x_2-(a/b)y_2)

看上面这个方程,一组必然存在的解出现了!

x = y_2, y = x_2 - (a/b)y_2 ⑤

可见我们只要求出 x_2,y_2,就能得出正确的x,y。问题是 x_2,y_2怎么求。

脑海里抛掉前面的x,y,现在我们手上是 b,a \mod b这两个系数,而目标是求出满足 bx_2 + (a \mod b)y_2 = gcd(b,a \mod b)(③)的 x_与 y_2,以便于待会回馈给⑤。

等一下。

比较于③,我再去看看原方程①:

ax + by = gcd(a,b)(①)

原方程中的 a 变成了③中的 b,原方程中的 b 变成了③中的 a \mod b而已。

按照上面的一模一样下来(其实都只是推导过程),我们发现,最好有 x_3,y_3来支撑 x_2, y_2。

再一模一样下来,我们又需要 x_4,y_4来支撑 x_3, y_3。

……

这个递归中 a,b 不断变成 b, a \mod b,逼近最后:

b 将等于最早的 a,b 的最大公因数,就像普通 gcd 的结果。

但我们再往下一层。此时由于 a \mod b = 0,下一层的 b = 0。

是时候直接回馈了,我们需要一组 x_n,y_n满足 a_nx_n + b_ny_n = gcd(a_n,b_n)。

因为gcd(a_n,b_n)=gcd(a,b)所以我们需要一组x_n,y_n满足 a_nx_n + b_ny_n = gcd(a,b)。而此时a_n=gcd(a,b)所以只要取x=1即可成为答案。最后一层回馈的 y_n建议用 0,但由于 b = 0,回馈其它数值等式一定也成立。意料之外的是,这样的程序有时候会求错解。如果回馈 3,那么你将收获 70 分的好成绩。如果回馈 -20002,那么你将收获 40 分的好成绩。为什么呢?

原因仅仅是,被回馈(却与 0 等效)的 y_n在回溯时滚雪球式增长,容易数值越界。

最后一层结束后,就开始返回,直到最上层。每一层可以轻松地根据下层的 x_{k+1},y_{k+1}求出当前层的 x_k, y_k。

当我们求出了第一层的x时,我们就可以根据x算出y的值,这就是方程的一组可行解。

如果要求x的最小正整数解,则将x批量地减去或加上b。
证明:

ax+by=1

ax + by + k×ba - k×ba = 1

a(x+kb) + (y-ka)b = 1

1.显然这并不会把方程中任何整数变成非整数。

2.“加上或减去 bb”不会使 xx 错过任何解。

小提示:

上面的算法中,以辗转相除的方式向下递进的好处(目的)是缩小系数,直到出现一个解可以确定的情况。

以下是 [p1082同余方程](https://www.luogu.org/problemnew/show/P1082)的程序。

#include<bits/stdc++.h>
using namespace std;
long long x,y;
int aa,bb;
void exgcd(long long a,long long b)
{
	if(b==0)
	{
		x=1;
		y=0;
		return;
	}
	exgcd(b,a%b);
	long long cunx=x,cuny=y;//cunx yiqiuchudex
	x=cuny;
	y=cunx-a/b*cuny;
}
int main()
{
	cin>>aa>>bb;
	exgcd(aa,bb);
	x=((x%bb)+bb)%bb;
	cout<<x;
	return 0;
}

  

posted @ 2019-07-23 10:31  zhenyan2003  阅读(320)  评论(0编辑  收藏  举报