欧几里德与扩展欧几里德算法详解

欧几里得算法：

定理1：设a,b,c,q都为整数，且b>0。如果 a = q*b+c，那么 gcd(a, b) = gcd(b, c)

证明方法用了集合的方法，就是说明一个的约数必定是另一个数的约数，从而两个数相等。

证明过程：由于可以写成a = q*b+c（a>b），那么设a和b的最大公约数为g；

证明g为c的公约数：因为a%g==0,所以(q*b+c)%g==0,又应为b%g==0,所以c%g==0,所以g是c的公约数;

证明g为b和c的最大公约数：如果g不是b和c的最大公约数，那么存在g'是b和c的最大公约数，那么g'也就是a和b的最大公约数，所以与已知矛盾;

因此利用欧几里德定理求a和b最大公约数代码如下：

template<typename T>
T gcd(T a,T b)
{
if(a < b) swap(a,b);
if(b == 0) return a;
return gcd(b,a%b);
}

拓展的欧几里得算法（Extended Euclidean Algorithm）

可以用来求二元一次方程组的整数解问题:a*x+b*y+c=0

步骤1直接使用公式a*x+b*y=gcd(a,b)-->步骤2求出gcd(a,b)-->步骤3求gcd(a,b)的同时可以获得a*x+b*y1=gcd(a,b)的一组原始解(x1,y1)-->步骤4使用x=x1+b/gcd(a,b)*t,y=y2-a/gcd(a,b)*t(其中t可以为任意整数)｛利用参数方程变形而来｝-->（x,y）即为所求

证明过程：令a1=a/gcd(a,b)，b1=b/gcd(a,b)，c=c/gcd(a,b)。如果我们能够首先求出满足a*x1+b*y1=gcd(a,b)这个方程的x1和y1，那么x=x1*(-c),y=y1*(-c){由于可以a*x+b*y+c=0写成a*x+b*y=-c}就可以求出来了。由欧几里德算法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的值求了出来。(x1,y1)就是原始解，通过参数方程可以得到x=x1+t1(其中t1为任意整数)，y=y1+t2(其中t2为任意整数)，代入方程a*x+b*y=gcd(a,b)-->a*(x1+t1)+b*(x1+t2)=gcd(a,b)-->a*x1+b*x1+a*t1+b*t2=gcd(a,b),其中a*x1+b*x2=gcd(a,b)->a*t1+b*t2=0则得到t1与t2的关系，由于t1和t2必须是整数，所以t1必须是b/gcd(a,b)的整数倍，所以就有x=x1+b/gcd(a,b)*t，将x代入方程则能算出y=y2-a/gcd(a,b)*t；

扩展欧几里德算法的代码：

template<typename T>
T extendGcd(T a,T b,T &x,T &y)
{
    if(b==0)
    {
       x=1,y=0;
       return a;
    }
    else
    {
        T g=gcd(b,a%b,y,x);
        y-=a/b*x;
        return g;
    }
}

例题：来自SUG的106题

There is an equation ax + by + c = 0. Given a,b,c,x1,x2,y1,y2 you must determine, 
how many integer roots of this equation are satisfy to the following conditions: 
x1<=x<=x2,y1<=y<=y2. 
Integer root of this equation is a pair of integer numbers (x,y).
Input
Input contains integer numbers a,b,c,x1,x2,y1,y2 delimited by spaces and line breaks. All numbers are not greater than 108 by absolute value.
Output
Write answer to the output.
Sample Input
1 1 -3
0 4
0 4
Sample Output
4

代码(未写main函数，可以直接在main函数中调用Equation返回的即是输出结果)：

#define Long_int long long int
#define LinearEquation(x,y,a,b,c) ((a)*(x)+(b)*(y)+(c))
template<typename T>
T gcd(T a,T b,T &x,T &y)
{
	if(b==0)
    {
       x=1,y=0;
       return a;
    }
    else
    {
        T g=gcd(b,a%b,y,x);
        y-=a/b*x;
        return g;
    }
}
//数据交换
template<typename T>
void sawp(T &t1,T &t2)
{
	T temp;
	temp = t1;t1 = t2;t2 = temp;
}
//主程序
Long_int Equation(void)
{
	Long_int a,b,c,temp;
	cin>>a>>b>>c;
	Long_int arr_x[2],arr_y[2];
	cin>>arr_x[0]>>arr_x[1];
	cin>>arr_y[0]>>arr_y[1];
	
	//找到坐标的上下界
	if(arr_x[0] > arr_x[1]) swap(arr_x[0],arr_x[1]);
	if(arr_y[0] > arr_y[1]) swap(arr_y[0],arr_y[1]);
	//判断条件1
	if((0==a)&&(0==b))
	{
		if(0 == c)
			return (arr_x[1]-arr_x[0]+1)*(arr_y[1]-arr_y[0]+1);
		else
			return 0;
	}
	//判断条件2
	if(0 == a)
	{
		if(0 == abs(c)%b)
		{
			temp = -1*c/b;
			if((temp>=arr_y[0])&&(temp<=arr_y[1]))
			{
				return (arr_x[1]-arr_x[0]+1);
			}
		}
		else
			return 0;
	}
	if(0 == b)
	{
		if(0 == abs(c)%a)
		{
			temp = -1*c/a;
			if((temp>=arr_x[0])&&(temp<=arr_x[1]))
			{
				return (arr_y[1]-arr_y[0]+1);
			}
		}
		else
			return 0;
	}
	//判断条件3
	Long_int x0,y0;
	Long_int d = gcd(a,b,x0,y0);
	if(c%d == 0)
	{
		c /= d;
		x0 *= -c;y0 *= -c;
		Long_int kx1,kx2,ky1,ky2;
		Long_int x1 = arr_x[0],x2 = arr_x[1],y1 = arr_y[0],y2 = arr_y[1];
		kx1 = (x1<=x0||(x1-x0)*d%b==0)?(x1-x0)*d/b:(x1-x0)*d/b+1;
		kx2 = (x2>=x0||(x0-x2)*d%b==0)?(x2-x0)*d/b:(x2-x0)*d/b-1;
		ky1 = (y1<=y0||(y1-y0)*d%a==0)?(y0-y1)*d/a:(y0-y1)*d/a-1;
		ky2 = (y2>=y0||(y0-y2)*d%a==0)?(y0-y2)*d/a:(y0-y2)*d/a+1;
		if(kx1 > kx2) swap(kx1,kx2);
		if(ky1 > ky2) swap(ky1,ky2);
		
		Long_int k_low,k_up;
		//cout<<" "<<kx1<<" "<<kx2<<" "<<ky1<<" "<<ky2<<endl;
		k_low = kx1>ky1?kx1:ky1;//cout<<"the k_low : "<<k_low<<endl;
		k_up = kx2>ky2?ky2:kx2;//cout<<"the k_up : "<<k_up<<endl;
		return (k_up-k_low)<0?0:(k_up-k_low+1);
	}
	else
		return 0;
}