求解最大公約數——歐幾里得演算法及其(解同餘方程)拓展歐幾里得
阿新 • • 發佈:2019-02-08
/* 擴充套件歐幾里得定理 擴充套件歐幾里得定理:對於兩個不全為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。