MAT-欧几里得及拓展欧几里得算法

欧几里得(Euclid)古希腊数学家，被称为"几何之父"，他最著名的著作《几何原本》是欧洲数学的基础，总结了平面几何五大定理，同时，在透视、圆锥曲线、球面几何学及数论方面也有突出贡献。在介绍算法之前，让我们首先对历史上为了人类的发展孜孜不倦，坚持真理的先辈们致以真诚的敬意和缅怀。

1.欧几里得算法

又称辗转相除法，用于计算两个正数a,b的最大公约数。算法的计算依赖于下面的定理：

gcd(a,b) = gcd(b,a mod b) (a>b 且 a mod b不为0)

证明：
a可以表示为a=kb+r,其中a,b,k,r都是整数，则r=a mod b 。假设d是a,b的一个公约数,则有d|a,d|b (即a,b能被d整除)， 而r=a-kb，所以d也能整除r,记 d|r.  因此d也是(b,a mod b)   的公约数， 因此(a,b)与(b,a mod b)的公约数是相等的，得证。gcd(a,b)=gcd(b,a mod b)

举个形象点的例子，假设a=16,b=12,按照算法所描述的，a=b+4 假设d为其公约数,那么有16=d*k1,12=d*k2,r = a-b = d(k1-k2) 即d|(a-b) = d|r,则  gcd(16,12) = gcd(12,4)；
同理可得 gcd(12,4)= gcd(4,0)  这里4就是a,b的最大公约数。

算法实现：

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

2.欧几里得拓展算法

定理：对于不完全为0的非负整数a,b,gcd(a,b)表示a,b的最大公约数，必然存在整数对x,y使得gcd(a,b)=ax+by

证明：
设a>b，当b=0,gcd(a,b)=a，此时x=1,y=0
当a,b!=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
=ay2+b(x2-(a/b)y2)
即x1=y2;y1=x2-(a/b)y2;
基于上面的事实，我们可以基于x2,y2的值使用递归来计算x1,y1,直到b=0时结束。

同样的，基于上面的事实，我们看个例子a=16,12,我们的目标是求解
gcd(16,12) = 16x+12y;
根据假设有：
16x1+12y1=gcd(16,12);
12x2+(16%12)y2=gcd(12,16%12)=gcd(12,4);
4x3+(12%4)y3=gcd(4,12%4)=gcd(4,0)
可以得到:
x1=y2; y1=x2-y2;
x2=y3; y2=x3-3y3;
x3=1; y3=0;
我们知道b=0时，递归停止，那么可以反推出
gcd(4,0)=4*1+0*0;  x3=1,y3=0;
gcd(12,4)=12x+4*y; x2=0,y2=1;
gcd(16,12)=16x+12y;x1=1,y1=-1;
算法实现:

int ext_gcd(int a,int b,int &x,int&y)
{
	if(b==0){
		x=1;y=0;
		return a;
	}
	int d=ext_gcd(b,a%b,x,y);
	int t=x;
	x=y;
	y=t-a/b*y;
	return d;
}

3.拓展欧几里得算法应用:

(1)求解不定方程

<1>.对于不定整数ax+by=c,如果c mod gcd(a,b)=0,则方程存在整数解，否则不存在整数解。
给方程两端同除以gcd(a,b)得到新的不定式a'x+b'y=c',gcd(a',b')=1(因为a,b已除过最大公约数)；
<2>.利用拓展的欧几里得算法求出a'x+b'y=1的一组整数解x0,y0，则n'x0,n'y0是方程a'x+b'y=n'的一组整数解;
这里的n'=c/gcd(a,b);
<3>.可以得到方程a'x+b'y=n'的所有其他整数解满足:
x=n'x0+b't;

y=n'y0-a't;(t为任意整数)

求解不定方程ax+by=c算法实现:

bool exgcd_equation(int a,int b,int c,int &x,int &y)
{
	int d = ext_gcd(a,b,x,y);
	if(c%d)
		return false;
	int t= c/d;
	x*=t;y*=t;
	return true;
}

(2)求解模线性方程

同余方程 ax≡b (mod n)对于未知数 x 有解，当且仅当 gcd(a,n) | b。且方程有解时，方程有 gcd(a,n) 个解。求解方程 ax≡b (mod n)相当于求解方程ax+ny=b, 设d=gcd(a,n),假如 整数x和y，满足d=ax+ny(拓展的欧几里得)，则存在解x0,y0使得方程满足 a*x0+n*y0=d,方程两边同乘以b/d,有 a*x0*b/d+n*y0*b/d=b，所以x=x0*b/d,y=y0*b/d为ax+ny=b的一个解，所以x=x0*b/d为ax≡b (mod n)的解，方程ax≡b (mod n)的一个解为 x0=x*(b/d)mod n,其他d个解分别为xi=(x0+i*(n/d))mod n{i=0,1,...d-1}

 我们看到解的之间的间隔为n/d，那么是如何得到呢？

 假设间隔为dx,那么有

 a*x=b(mod n);
 a*(x+dx)=b(mod n);
 两式相减，得到:
 a*dx(mod n) = 0;
也就是说a*dx是a的倍数，同时也是n的倍数，即a*dx是a和n的公倍数，为了求出dx，我们应该求出a和n的最小公倍数，此时对应的dx是最小的。

 设a和n的最大公约数为d，那么a和n的最小公倍数为(a*n)/d.

 即a*dx=a*n/d;所以dx =n/d ，这就是解的间隔。

 bool modular_equation(int a,int b,int n)
 {
	int x,y,x0,i;
	int d = exgcd(a,n,x,y);
	if(b%d)
		return false;
	x0 = x*(b/d)%n;
	for(i=1;i<d;i++)
		printf("%d\n"(x0+i*(n/d))%n);
	return true;
 }

 (3)求模的逆元
 同余方程ax≡b (mod n),如果gcd(a,n)=1,则方程只有唯一解。

  这时如果b=1,同余方程就是ax=1(mod n),gcd(a,n)=1.称求出的x为a的对模n乘法的逆元。对于同于方程ax=1(mod n),gcd(a,n)=1的求解就是求解方程ax+ny=1，x,y为整数。这个   可用拓展欧几里得算法求出，原同余方程的唯一就是用拓展欧几里得算法得出的x。