扩展欧几里得算法学习笔记

Posted 2021-05-19 llmmkk

扩展欧几里得算法学习笔记：

---

前言：

学了两周数据结构发现数论图论忘光了，所以回来补一下，顺便写下笔记。

---

前置需要：

欧几里得算法，裴蜀定理，脑子

欧几里得算法:即辗转相除法，\\(\\gcd(a,b)=\\gcd(b,a \\bmod b)\\)

裴蜀定理：若 \\(a,b\\) 是整数,且 \\(\\gcd(a,b)=d\\)，那么对于任意的整数\\(x\\)和\\(y\\),\\(ax+by\\) 都一定是\\(d\\)的倍数，特别地，一定存在整数\\(x,y\\),使\\(ax+by=d\\)成立。

脑子

---

用途：

扩展欧几里得算法\\((exgcd)\\)一般用于判断并求解形如 \\(ax+by=c\\) 的不定方程。（\\(a,b,c\\) 为已知数，\\(x,y\\) 为未知数）

---

思想：

首先由裴蜀定理可知，方程有解必须满足 \\(\\gcd(a,b)\\mid c\\),所以先讨 \\(ax+by=\\gcd(a,b)\\)，再推至所有情况。

可以由欧几里得算法往后推：
\\(ax+by=\\gcd(a,b)\\)
\\(=\\gcd(b,a\\bmod b)\\)
\\(=\\gcd(b,a-\\left \\lfloor \\frac{a}{b} \\right \\rfloor*b)\\)
\\(=bx_1+(a-\\left \\lfloor \\frac{a}{b} \\right \\rfloor*b)y_1\\)
\\(=bx_1+ay_1-\\left \\lfloor \\frac{a}{b} \\right \\rfloor *by_1\\)
\\(=ay_1+b(x_1-\\left \\lfloor \\frac{a}{b} \\right \\rfloor*y_1)\\)

发现只要有了 \\(\\gcd(b,a\\bmod b)\\)时的 \\(x_1,y_1\\) ，就可以推出 \\(\\gcd(a,b)\\) 时的 \\(x\\) 和 \\(y\\)。即：
\\(x=y_1,y=x_1-\\left \\lfloor \\frac{a}{b} \\right \\rfloor*y_1\\)

当欧几里得算法算到最后即 \\(\\gcd(a,0)\\) 时返回\\(a\\)为\\(\\gcd(a,0)\\)，显然此时有\\(a*1+0*0=\\gcd(a,0)\\)，即: \\(x=1,y=0\\)

综上可以发现，我们在递归 \\(\\gcd(a,b)\\) 回溯时同时返回当前的 \\(x_1,y_1\\)，就可以求出上一层的 \\(x,y\\)，这样就可以通过递归求解最开始的 \\(x,y\\) 了。

该算法依据于欧几里得算法 \\((gcd)\\)，因而得名扩展欧几里得算法 \\((exgcd)\\)。

---

实现:

#include<bits/stdc++.h>
using namespace std;
int exgcd(int a,int b,int &x,int &y){
	if(a<b)return exgcd(b,a,y,x);
	if(b==0){
		x=1;y=0;
		return a;
	}
	else{
		int tx,ty;
		int d=exgcd(b,a%b,tx,ty);
		x=ty;
		y=tx-a/b*ty;
		return d;
	}
}
int main(){
	int a,b,x,y;
	cin>>a>>b;
	int d=exgcd(a,b,x,y);
	cout<<x<<" "<<y<<endl;
	return 0;
}

---

另外这里

	int tx,ty;
	int d=exgcd(b,a%b,tx,ty);
	x=ty;
	y=tx-a/b*ty;

可以写成更美观的形式，一样的含义，即

	int x1;
	int d=exgcd(b,a%b,x1,x);
	y=x1-a/b*x;

不难理解。

---

例题：

1.P1082 [NOIP2012 提高组] 同余方程

推柿子：
\\(ax\\equiv 1\\pmod b\\)
\\(ax-1=bt\\)
\\(ax-bt=1\\)

发现是道模板题，但怎么保证 \\(x\\) 是最小整数解呢？
由于我们求出的 \\(x,y\\) ,是不定方程的一组特解，所以来考虑通解。
对于方程\\(ax+by=\\gcd(a,b)\\)，可以稍微变形得到

\\(ax+k*ab+by-k*ab=\\gcd(a,b)\\)

\\(a(x+kb)+b(y-ka)=\\gcd(a,b)\\)

有一组整数解为\\(x_0,y_0\\)

则方程通解为
\\(\\begin{cases} x=x_0+kb\\\\ y=y_0-ka \\end{cases}\\) \\((k\\in \\mathbb{Z})\\)
也就是说 \\(x\\) 是每 \\(b\\) 个数中存在一个，\\(y\\) 是每 \\(a\\) 个数中存在一个。
看原式显然 \\(x\\) 是不等于0或 \\(b\\) 的，所以只需取到 \\(1\\) 到 \\(b-1\\) 之间的 \\(x\\)。
具体操作只需要 \\(x\\) 对 \\(b\\) 取模一次，当然 \\(x\\) 可能为负，需要再加一个 \\(b\\) 变成正，而正的 \\(x\\) 加 \\(b\\) 又大于了 \\(b\\)，再取一次模即可，可见代码。

另外，方程有解的情况为 \\(\\gcd(a,b)\\mid 1\\)，题目又保证有解，则 \\(\\gcd(a,b)=1\\)，也就是说 \\(a,b\\) 互质。（虽然并没有什么用。

---

代码：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
ll exgcd(ll a,ll b,ll &x,ll &y){
	if(a<b)return exgcd(b,a,y,x);
	if(b==0){
		x=1;y=0;
		return a;
	}
	else{
		int x;
		int d=exgcd(b,a%b,x1,x);
		y=x1-a/b*x;
		return d;
	}
}
int main(){
	ll a,b,x,y;
	cin>>a>>b;
	ll d=exgcd(a,b,x,y);
	cout<<(x%b+b)%b;
	return 0;
}

---

2.P1516 青蛙的约会

设跳 \\(s\\) 次后两青蛙相遇，可以列出同余方程:

\\(x+sm\\equiv y+sn \\pmod L\\)

\\(x+sm-y-sn=tL\\)

由于 \\(s\\) 和 \\(t\\) 是未知数，所以化为

\\(s(n-m)+tL=x-y\\)

发现这就是不定方程，设 \\(a=n-m,b=L,c=x-y\\)

我们就得到了\\(sa+tb=c\\)

但此题 \\(c\\) 是 \\(\\gcd(a,b)\\) 的倍数，只需先求出
\\(sa+tb=\\gcd(a,b)\\) 的解，再方程两边同乘\\(\\frac{c}{\\gcd(a,b)}\\),即原方程答案为 \\(s\\frac{c}{\\gcd(a,b)}\\) 和 \\(t\\frac{c}{\\gcd(a,b)}\\)。

这就做完了吗？并没有。注意到 \\(gcd\\) 对负数是没有意义的，所以我们要保证 \\(\\gcd(a,b)\\) 中的 \\(a,b\\) 都不为负，\\(b=L\\),而 \\(L\\) 是保证大于0的，但 \\(a=n-m\\) 是可能为负的，所以方程两边要同乘 \\(-1\\) ，但为了保证 \\(b\\) 为负，所以我们只用把 \\(a\\) 和 \\(c\\) 都取负，对于 \\(b\\)，我们不操作。
这样方程不就不是原来的方程了吗，是的，但他们是有联系的。原本的方程 \\(sa+tb=c\\) 中\\(s\\) 和 \\(t\\) 相当于一般形式中的 \\(x\\) 和 \\(y\\) ，而现在的方程变为了 \\(-sa+tb=-c\\) 即 \\(sa-tb=c\\) ，
其中 \\(s\\) 和 \\(-t\\) 相当于一般形式中的 \\(x\\) 和 \\(y\\) 。
所以只需对最后得出的 \\(t\\) 取负就能得出原方程的解，另外也可以将 \\(s\\) 带入反推 \\(t\\)。当然本题对 \\(t\\) 没有要求，但不代表以后不会求类似的解。

理解了这一点代码就是模板了。

---

代码：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
ll exgcd(ll a,ll b,ll &x,ll &y){
	ll d,t;
	if(b==0){
		x=1;y=0;
		return a;
	}
	d=exgcd(b,a%b,x,y);
	t=x-a/b*y;
	x=y;
	y=t;
	return d;
}
ll x,y,m,n,l,a,b,c,tx,ty,d,r;
int main(){
	cin>>x>>y>>m>>n>>l;
	a=n-m,b=l,c=x-y;
	if(a<0){
		a=-a;
		c=-c;
	}
	d=exgcd(a,b,tx,ty);
	r=b/d;
	if(c%d==0) cout<<((c/d*tx)%r+r)%r;
	else cout<<"Impossible"<<endl;
	return 0;
}

---

题外话：

写了这么多，自己又分析了一通理解++了，希望以后不会忘。