对于不完全为 0 的非负整数 a,b,gcd(a,b)表示 a,b 的最大公约数,必然存在整
数对 x,y ,使得 gcd(a,b)=ax+by。
c++语言实现
#include<iostream>
#include<cstdio>
usingnamespacestd;
intx,y,q;
voidextend_Eulid(inta,intb){
if(b==0){
x=1;y=0;q=a;
}else{
extend_Eulid(b,a%b);
inttemp=x;x=y;y=temp-a/b*y;
}
}
intmain(){
inta,b;
cin>>a>>b;
extend_Eulid(a,b);
printf("%d=(%d)*%d+(%d)*%d\n",q,x,a,y,b);
return0;
}
求解 x,y的方法的理解
设 a>b。
1,显然当 b=0,gcd(a,b)=a。此时 x=1,y=0;
2,ab<>0 时
设 ax1+by1=gcd(a,b);
bx2+(a mod b)y2=gcd(b,a mod b);
根据朴素的欧几里德原理有 gcd(a,b)=gcd(b,a mod b);
则:ax1+by1=bx2+(a mod b)y2;
即:ax1+by1=bx2+(a-[a/b]*b)y2=ay2+bx2-[a/b]*by2;
也就是ax1+by1==ay2+b(x2-[a/b]*y2);
根据恒等定理得:x1=y2; y1=x2-[a/b]*y2;
这样我们就得到了求解 x1,y1 的方法:x1,y1 的值基于 x2,y2.
上面的思想是以递归定义的,因为 gcd 不断的递归求解一定会有个时候 b=0,所以递归可以结束。
扩展欧几里德算法
扩展欧几里德算法是用来在已知a, b求解一组x,y使得ax+by = Gcd(a, b) =d(解一定存在,根据数论中的相关定理)。扩展欧几里德常用在求解模线性方程及方程组中。下面是一个使用C++的实现:
int exGcd(int a, int b, int &x, int &y){if(b == 0){x = 1;y = 0;return a; //---很难找出一个这么实现的价值,因为扩展欧几里得还有更大的用途;个人认为定义全局数组更好,不用return r。}int r = exGcd(b, a % b, x, y);int t = x;x = y;y = t - a / b * y;return r;}
把这个实现和Gcd的递归实现相比,发现多了下面的x,y赋值过程,这就是扩展欧几里德算法的精髓。
可以这样思考:
对于a' = b, b' = a % b 而言,我们求得 x, y使得 a'x + b'y = Gcd(a', b')
由于b' = a % b = a - a / b * b (注:这里的/是程序设计语言中的除法)
那么可以得到:
a'x + b'y = Gcd(a', b') ===>
bx + (a - a / b * b)y = Gcd(a', b') = Gcd(a, b) ===>
ay +b(x - a / b*y) = Gcd(a, b)
因此对于a和b而言,他们的相对应的p,q分别是 y和(x-a/b*y)
使用扩展欧几里德算法解决不定方程的办法
对于不定整数方程pa+qb=c,若 c mod Gcd(a, b)=0,则该方程存在整数解,否则不存在整数解。
上面已经列出找一个整数解的方法,在找到p * a+q * b = Gcd(a, b)的一组解p0,q0后,/*p * a+q * b = Gcd(a, b)的其他整数解满足:
p = p0 + b/Gcd(a, b) * t
q = q0 - a/Gcd(a, b) * t(其中t为任意整数)
至于pa+qb=c的整数解,只需将p * a+q * b = Gcd(a, b)的每个解乘上 c/Gcd(a, b) 即可
在找到p * a+q * b = Gcd(a, b)的一组解p0,q0后,应该是
得到p * a+q * b = c的一组解p1 = p0*(c/Gcd(a,b)),q1 = q0*(c/Gcd(a,b)),p * a+q * b = c的其他整数解满足:
p = p1 + b/Gcd(a, b) * t
q = q1 - a/Gcd(a, b) * t(其中t为任意整数)
p 、q就是p * a+q * b = c的所有整数解。
那么在a*x+b*y=m这个线性方程成立的情况下,如何来求解x和y呢?
1.令a1=a/gcd(a,b),b1=b/gcd(a,b),m1=m/gcd(a,b)。如果我们能够首先求出满足a*x1+b*y1=gcd(a,b)这个方程的x1和y1,那么x=x1*m1,y=y1*m1就可以求出来了。由欧几里德算法gcd(a,b)=gcd(b,a%b),所以a*x1+b*y1=gcd(a,b)=gcd(b,a%b)=b*x2+(a%b)*y2,现在只要做一些变形就可以得到扩展欧几里德算法中的用到的式子了。令k=a/b(商),r=a%b(余数),那么a=k*b+r。所以r=a-k*b,带入上式,得到a*x1+b*y1=b*x2+(a-(a/b)*b)y2=a*y2+b*(x2-(a/b)*y2) => x1=y2,y1=x2-(a/b)*y2。有了这两个式子我们就知道了在用欧几里德求最大公约数的时候,相应的参数x,y的变化。现在再回过头来看一下扩展欧几里德算法的代码就很好理解了,实际上扩展欧几里德就是在求a和b的最大公约数的同时,也将满足方程a*x1+b*y1=gcd(a,b)的一组x1和y1的值求了出来。下面代码中突出的部分就是标准的欧几里德算法的代码。
__int64 exGcd(__int64 a,__int64 b,__int64 &x,__int64 &y){
if(b==0){
x=1;
y=0;
return a;
}
__int64 g=exGcd(b,a%b,x,y);
__int64 temp=x;
x=y;
y=temp-(a/b)*y;
return g;
}
2.那么x,y的一组解就是x1*m1,y1*m1,但是由于满足方程的解无穷多个,在实际的解题中一般都会去求解x或是y的最小正数的值。以求x为例,又该如何求解呢?还是从方程入手,现在的x,y已经满足a*x+b*y=m,那么a*(x+n*b)+b*(y-n*a)=m显然也是成立的。可以得出x+n*b(n=…,-2,-1,0,1,2,…)就是方程的所有x解的集合,由于每一个x都肯定有一个y和其对应,所以在求解x的时候可以不考虑y的取值。取k使得x+k*b>0,x的最小正数值就应该是(x+k*b)%b,但是这个值真的是最小的吗??如果我们将方程最有两边同时除以gcd(a,b),则方程变为a1*x+b1*y=m1,同上面的分析可知,此时的最小值应该为(x+k*b1)%b1,由于b1<=b,所以这个值一定会小于等于之前的值。在实际的求解过程中一般都是用while(x<0)x+=b1来使得为正的条件满足,为了更快的退出循环,可以将b1改为b(b是b1的倍数),并将b乘以一个倍数后再加到x上。
编程时 exgcd 更多用于求解“中国余数定理”相关知识 举个例子比如n除以5余2 除以13余3 那么n最小是多少,所有的n满足什么条件?
n(min)=42
n=42+k*65
欧几里德算法的扩展
扩展欧几里德算法不但能计算(a,b)的最大公约数,而且能计算a模b及b模a的乘法逆元,用C语言描述如下:
intgcd(inta,intb,int&ar,int&br){
intx1,x2,x3;
inty1,y2,y3;
intt1,t2,t3;
if(0==a){//有一个数为0,就不存在乘法逆元
ar=0;br=0;
returnb;
}
if(0==b){
ar=0;
br=0;
returna;
}
x1=1;x2=0;x3=a;y1=0;y2=1;y3=b;intk;
for(t3=x3%y3;t3!=0;t3=x3%y3){
k=x3/y3;
t2=x2-k*y2;t1=x1-k*y1;
x1=y1;
x2=y2;
x3=y3;
y1=t1;
y2=t2;
y3=t3;
}
if(y3==1){//有乘法逆元
ar=y2;br=x1;
return1;
}
else{//公约数不为1,无乘法逆元
ar=0;br=0;returny3;
}
}
扩展欧几里德算法对于最大公约数的计算和普通欧几里德算法是一致的。计算乘法逆元则显得很难明白。我想了半个小时才想出证明他的方法。
首先重复操作整除中的一个论断:
如果gcd(a,b)=d,则存在m,n,使得d = ma + nb,称呼这种关系为a、b组合整数d,m,n称为组合系数。当d=1时,有 ma + nb = 1 ,此时可以看出m是a模b的乘法逆元,n是b模a的乘法逆元。
为了证明上面的结论,我们把上述计算中xi、yi看成ti的迭代初始值,考察一组数(t1,t2,t3),用归纳法证明:当通过扩展欧几里德算法计算后,每一行都满足a×t1 + b×t2 = t3
第一行:1 × a + 0 × b = a成立
第二行:0 × a + 1 × b = b成立
假设前k行都成立,考察第k+1行
对于k-1行和k行有
t1(k-1) t2(k-1) t3(k-1)
t1(k) t2(k) t3(k)
分别满足:
t1(k-1) × a + t2(k-1) × b = t3(k-1)
t1(k) × a + t2(k) × b = t3(k)
根据扩展欧几里德算法,假设t3(k-1) = j t3(k) + r
则:
t3(k+1) = r
t2(k+1) = t2(k-1) - j × t2(k)
t1(k+1) = t1(k-1) - j × t1(k)
则
t1(k+1) × a + t2(k+1) × b
=t1(k-1) × a - j × t1(k) × a +
t2(k-1) × b - j × t2(k) × b
= t3(k-1) - j t3(k) = r
= t3(k+1)
得证
因此,当最终t3迭代计算到1时,有t1× a + t2 × b = 1,显然,t1是a模b的乘法逆元,t2是b模a的乘法逆元。