1. 程式人生 > >求解最大公約數——歐幾里得演算法及其(解同餘方程)拓展歐幾里得

求解最大公約數——歐幾里得演算法及其(解同餘方程)拓展歐幾里得

/*
  擴充套件歐幾里得定理
  擴充套件歐幾里得定理:對於兩個不全為0的整數a、b,必存在一組解x,y,
  使得ax+by==gcd(a,b);
  拓展歐幾里得實現
  下面面的程式中,x和y用全域性變數儲存
  int gcd(int a,int b)
  {
    int t,d;
    if(b==0)
    {
        x=1;
        y=0; 
	//此時b==0,也就是說gcd(a,0)==a。原式變為ax+by==a=gcd(a,b)--> x==1,y==0
        return a; //返回a,b最大公約數的值 
    }
    d=gcd(b,a%b); //歐幾里得求最大公約數應用 
    t=x;
    x=y;
    y=t-(a/b)*y;  //不明處2
    return d;     //返回a,b最大公約數的值
  } 
  //不明處2 解釋 ax+by==gcd(a,b)(1) 
    說明規則,x,y表示第一次遞迴時的值,x1,y1表示第二次遞迴時的值。
	那麼gcd(a,b)==gcd(b,a%b),同時都代入式1,
	有ax+by==b*x1+(a%b)*y1。將右邊變形一下
    b*x1+(a%b)*y1==b*x1+(a-(a/b)*b)*y1==a*y1+b*(x1-(a/b)*y1),
	最終得到ax+by==a*y1+b*(x1-(a/b)*y1)
    也就是說: 
	上一深度的x等於下一深度的y1,
	上一深度的y等於下一深度的x1-(a/b)*y1。 
   *需要注意,上面推導時用的除法都是整型除法

    因此,得到了不定式ax+by==gcd(a,b)的一組解,x、y。
    那麼對於一般的不定式ax+by==c,它的解應該是什麼呢。
	很簡單,x1=x*(c/gcd(a,b)),y1=y*(c/gcd(a,b)) 
	
	再深入一點,就解出這麼一組解其實一般來說是解決不了什麼問題的,
	我們現在要得到所有的解,那麼這所有的解究竟是什麼呢?
	假設d=gcd(a,b). 那麼x=x0+b/d*t; y=y0-a/d*t;其中t為任意常整數。
  
*/ 


#include<stdio.h>
#include<stdlib.h>

/*
  求關於 x 的同餘方程 ax ≡ 1 (mod b)的最小正整數解。
  即求ax=mb+r 1=nb+r
  變形ax+(n-m)b=1,此方程即拓展歐幾里德的應用ax+by=gcd(a,b),(n-m相當於y)
  事實上ax ≡1(mod b) 有解的必要條件是gcd(a,b)|1,即gcd(a,b)=1; 
  使用拓展歐幾里得可知 ax+by=1(x,y是整數) 
*/

int  Extragcd(int a,int b,int *x,int *y)
{
	int d,t;
	if(b==0)  //遞迴呼叫終止條件,當根據歐幾里得輾轉相除法則,餘數為0停止 
	{
		*x=1;
		*y=0;
		return a;
	}
	else
	{
		d=Extragcd(b,a%b,x,y);
		t=*x;     //根據下一個x1,y1的值,倒推前一個x,y的值 
		*x=*y;
		*y=t-a/b*(*y);
		return d;
	}
}

int main()
{
	int a,b,x,y;
	int ans;
	freopen("mod.in","r",stdin);
	freopen("mod.out","w",stdout); 
	scanf("%d%d",&a,&b);
	ans=Extragcd(a,b,&x,&y);
	if(ans!=1) return 0;
	//根據若x是方程的一個解,則方程的所有解為x+k*b k為整數 
    x=x%b; //保證最小的正整數解x ,且x屬於{0,1,2,3...b-1} 
    while(x<=0)
        x+=b ;
    printf("%d\n",x);
    return 0;
} 



ax+by=c問題

問題:ax+by=c,已知a、b、c,求解使該等式成立的一組x,y。其中a、b、c、x、y均為整數

a,b的最大公約數為gcd(a,b)。如果c不是gcd(a,b)的倍數,則該等式無解,因為等式左邊除以gcd(a,b)是整數,

而等式右邊除以gcd(a,b)後為小數。(根據解方程的時候,在等式的左右兩邊同時除以非0的整數,等式依然成立)

因此,只有當c是gcd(a,b)的倍數的時候,該等式有解。這樣,可以通過計算使ax1+by1=gcd(a,b)成立的x1、y1,

然後有x=(c/gcd(a,b))*x1,y=(c/gcd(a,b))*y1,得到x,y。

問題就被轉換為求使ax+by=gcd(a,b)成立的一組x,y

。這可以用擴充套件歐幾里德演算法求解。如下:

根據歐幾里得的輾轉相除法,最終如果b為零,則gcd(a,b)=a,那麼x=1,y=0為一組解。(一組重要的特解)

如果b不為零,根據歐幾里德定理,可以設

ax1+by1=gcd(a,b)=gcd(b,a%b)=bx2+(a%b)y2=bx2+(a-(a/b)*b)y2

進一步化簡可以得到 -> a*y2+b*(x2-(a/b)*y2)

化簡後有x1=y2,y1=x2-(a/b)y2。因此x1,y1依賴於x2,y2,同理依次類推遞迴呼叫求出x3,y3,x4,y4……,

類似於輾轉相除,直到b=0時,求出xn,yn,便可以由後往前倒推出x1,y1的值。

擴充套件歐幾里德演算法:

// 擴充套件歐幾里德演算法,解gcd(a, b) = ax + by
// 結果儲存在x,y中,使用者呼叫時保證a、b、c都是整數
// 返回a,b的最大公約數
int extended_euclid(int a, int b, int &x, int &y)
{
    if(b == 0) // gcd(a, b) == a
    {
        x = 1;
        y = 0;
        return a;
    }
    int n = extended_euclid(b, a%b, x, y);
    int tmp = x;
    x = y;
    y = tmp - static_cast<int>(a/b)*y;
    return n;
}
使用者呼叫linear_equation求解線性方程:
// 等式ax+by=c,已知a、b、c,求x和y。
// 解該線性方程等同於解同餘式ax = c(mod b)
// 返回值表示是否有解,true有解,false無解
bool linear_equation(int a, int b, int c, int &x, int &y)
{
    int n = extended_euclid(a, b, x, y);
    if(c%n)
        return false;
    int k = c/n;
    x *= k;
    y *= k;
    return true;
}

linear_equation函式也可以用來解同餘式ax=c(mod b)。

由ax=c(mod b),可以得到ax = mb+r;c = nb+r。化簡可以得到ax+(n-m)b=c。呼叫linear_equation可以求出x。