1. 程式人生 > >擴充套件歐幾里德演算法解二元一次不定方程

擴充套件歐幾里德演算法解二元一次不定方程

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

已知兩個不完全為 0 的非負整數 a,b,必然存在整數對 x,y ,使它們滿足貝祖等式:

ax+by = gcd(a, b) =d

解一定存在,根據數論中的相關定理。下面給出程式碼:

int extgcd(int a, int b, int& x, int& y) {
    int gcd = a;
    if (b != 0) {
        gcd = extgcd(b, a % b, y, x);
        y -= (a / b) * x;
    }
    else {
        x = 1; y = 0;
    }
    return gcd;
}

解二元一次不定方程(ax + by = c):

擴充套件歐幾里德演算法很大的一個用處就是在解不定方程。在計算時一般先是求出一對特解,再根據x與y的變化比例,構造通解。

注意我們求的所有解,都是整數解。

我們先計算 ax+by = gcd(a, b)  的解,等式兩邊同時乘以

 c / gcd(a,b) 

 原式子就被構造成

a x(c / gcd(a,b) )+by(c / gcd(a,b) ) = c

所以就求出了不定方程的一對特解

X=x*c / gcd(a,b)\, ,\, Y=y*c / gcd(a,b)

這裡不存在 c / gcd(a,b) 不是整數的情況,因為當 c\, \, mod\, \, gcd(a,b) \neq 0 時,方程不存在整數解,不在考慮範圍內。

證明:

a=k1*gcd(a,b)\, ,\, b=k2*gcd(a,b)

ax+by=gcd(a,b)(x*k1+y*k2)=c

若 c\, \, mod\, \, gcd(a,b) \neq 0,則 x*k1+y*k2 非整數

求一組特解程式碼:

bool Indefinite_equation(int a, int b, int n, int& x, int& y) {    //求一組特解
    int gcd = extgcd(a, b, x, y);       
    if (n % gcd) return false;
    int k = n / gcd;
    x *= k;
    y *= k;
    return true;
}

現在我們來考慮構造通解( 特解是X、Y ):

由於已經滿足 aX+bY=c , 兩個未知數的關係類似於加權和固定,若一個增大,則另一個一定按比例減小,顯然滿足關係:

a(X+nb)+b(Y-na)=c

X 每變化 b 的整數倍,Y 就反向變化 a 的同等倍數。似乎構造出了通解表示式 X+nb 、Y-na 。

但是我們會發現一個問題,就是 a 、b 不一定是最小的步長,可能會 “跳過” 許多解。

為了求得最小步長,我們應該對 a 、b 同時除以某個整數,使得商也是整數,就求出了每次最小的變化量。很明顯,應除以

gcd(a, b)

所以,通解應該是

 X+n*b/gcd(a,b)\, \, ,\, \, Y-n*a/gcd(a,b)

求特殊的解:

最大/最小解

很多情況下,我們需要求的是 x 或 y 的最小正整數解。這裡以 x 的最小解為例:

有一種方案是在求出特解的基礎上,進行一個迴圈

  • 當特解大於零時:  while (x + b/gcd(a,b) > 0) ,x每次減 b/gcd(a,b)
  • 當特解不大於零時: while (x <= 0) ,x每次加 b/gcd(a,b)

但是這個方法存在效率問題,因為 b/gcd(a,b) 可能遠小於 x , 達到退出迴圈條件需要花很長時間。這裡我採用的是取模的方法

  • 當特解大於零時: x\, \, mod\, \, b/gcd(a,b)  便是小正整數解
  • 當特解不大於零時:x\, \, mod\, \, b/gcd(a,b) +b/gcd(a,b)  便是小正整數解

x 與 y 最接近的解

以求 x < y時 , x 與 y 最接近的解為例。通過 y - x 的值來考慮。設 x , y 的步長分別是 dx , dy ,則 y - x 的最短步長為 dx + dy。

  • 當特解 X < Y 時,\left \lfloor (y-x)/(dx+dy) \right \rfloor  便是通解中 n 的值
  • 當特解 X > Y 時,   \left \lfloor (y-x)/(dx+dy) \right \rfloor-1 便是通解中 n 的值。但若 (y - x) \, mod\, (dx + dy)=0 時,不需要減 1 

附上程式碼:

//求 x < y 時,x, y 最接近的解
bool Indefinite_equation(int a, int b, int n, ll& x, ll& y) {      
    int gcd = extgcd(a, b, x, y);
    if (n % gcd) return false;                         //無整數解
    int k = n / gcd, dx = b / gcd, dy = a / gcd;       //dx,dy分別為對應的最小變化值
    x *= k;
    y *= k;                                            //此為某特解
    ll cnt = (y - x) / (dx + dy);
    if (cnt < 0 && (y - x) % (dx + dy)) cnt --;        //處理x>y的情況
    x += dx * cnt;
    y -= dy * cnt;
    return true;
}