淺談擴充套件歐幾里得演算法

什麼是擴充歐幾里得？簡單的說，就是求關於x,y的方程 ax + by = gcd(a,b) 的所有整數解

---

現在我們來解決四個問題

• 什麼是裴屬定理，如何證明裴屬定理？

• 怎麼用擴充套件歐幾里得來求ax + by = gcd(a,b) 的特解？

• 怎麼求由特解推出其他的所有解？

• 擴充套件歐幾里得的三大應用

---

一、裴蜀定理內容即證明

內容

設a, b是不全為零的整數，則存在整數x, y 使得 ax + by = gcd(a, b).

證明：

1. a, b中假如有一個為0，則上述定理完全吻合，如a = 0, 則gcd(a, b) = b , 顯然成立

2. 若a, b 都不等於0

   設gcd(a, b) = d,則d != 0

   等式兩邊同時除以d，得：

   a₁x + b₁y = 1

   再利用輾轉相除法：gcd(a, b) --> gcd(b, a % b) --> ……

   設模出來的餘數叫做x, 則有：

   gcd(a₁, b₁) = gcd(b₁, r₁) = gcd(r₁, r₂) = …… = gcd(r_n-1, r_n) = 1

   把輾轉相除法中的運算展開，做成帶餘數的除法，得

   • a₁ = x₁b + r₁

   • b₁ = x₂r₁ + r₂

   • r₁ = x₃r₂ + r₃

     ……

   • r_n-3 = x_n-1r_n-2 + r_n-1

   • r_n-2 = x_nr_n-1 + r_n

   • r_n-1 = x_n+1r_n

   不妨令輾轉相除法在除到互質的時候退出則 r_n = 1 , 由倒數第二行得

   • r_n-2 = x_nr_n-1 + 1 ---> 1 = r_n-2 - x_nr_n-1

   由倒數第三行得：r_n-1 = r_n-3 - x_n-1r_n-2帶入上式

   得：1 = (1 + x_nx_n-1)r_n-2 - x_nr_n-3

   採取同樣的方法，逐個消去r_n-2,……, r₁

   最終會得到1 = f(x)a₁ + g(x)b₁

   其中f(x)和g(x) 是關於x的函式

   設f(x) = x, g(x) = y，則可以得到1 = a₁x + b₁y, 得證!

   正常人都不喜歡看上面寫的奇奇怪怪的推導，那就看下面我這個具體的推導叭，這裡是假設到 r4 = 1，然後往前推導滴

   

二、求特解

我們都知道，歐幾里得公式可以由這個式子表達

gcd(a, b) = gcd(b, a % b)

通過這個式子，我們可以不斷遞推到b = 0， 此時a即為a和b的最大公約數

現在我們對這個式子進行展開，看看有什麼好玩的嘛

gcd(a, b) = a * x1 + b * x2

gcd(b, a % b) = b * x2 + (a % b) * y2

其中a % b = a - (a / b) * b

故由第一行的歐幾里得公式得：

a * x1 + b * y1 = b * x2 + {a - (a / b) * b}y2

化簡得a * x1 + b * y1 = a * y2 + b * {x2 - (a / b) * y2}

根據待定係數法得到

• x1 = y2

• y1 = x2 - (a / b) * y2

也就是說，如果有了gcd(b, a % b)的解x2, y2就可以推出gcd(a, b)的解x1, y1

那麼這就和求gcd的過程一樣，一直推到最後x = 1， y = 0，就可以回溯回去，得到gcd(a, b)的特解x1, y1

特解其實是最小的一個解

二、怎麼利用特解推出其他所有整數解

先說結論：

\[x = x0 + k * \frac{b}{gcd(a,b)} \]

\[y = y0 - k * \frac{a}{gcd(a,b)} \]

任意的x + y = x0 + y0

x0, y0就是方程的特解

對於a * x + b * y = g來說，讓x增加b/gcd,讓y增加a / gcd

這樣就相當於加a * b / gcd，減去a * b / gcd,一加一減，最終不變

三、擴充套件歐幾里得三大應用

• 用擴充套件歐幾里得演算法解不定方程ax+by=c

如果 c % gcd(a, b) == 0，即c 是 gcd(a, b)的整數倍，則方程有解，且解就在上面的基礎上乘以$$\frac{c}{gcd(a, b)}$$即可

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

int main(){
    int a, b, x, y, gcd;
  	cin>>a>>b;
    e_gcd(a, b, gcd, x, y);
    cout<<x<<' '<<y<<endl;
    cout<<gcd<<endl;
    return 0;
}

• 用擴充套件歐幾里得演算法求解模線性方程ax≡b (mod n)

  對於模線性方程ax≡b (mod n)可以化簡為ax + ny = b，設d = gcd(a, n) 當且僅當b % d == 0時有解，且有d個解

  當b % d == 0時，設ax + ny = d,的特解為x,y

  等式兩邊同時乘以\(\frac{b}{d}\)，則方程變成\(ax\frac{b}{d} + ny\frac{b}{d} = b\)

  則原方程ax + ny = b的解為\(x' = x*\frac{b}{d},y' = y*\frac{b}{d}\)

  故ax≡b (mod n)的特解為 x0 = x * (b / d) % n

  d個解為xi= x0 + i*(n/d) mod n，{i = 0……n-1}

  解的間隔為n/d

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

bool modular_linear_equation(int a, int b, int n){
    int x, y, x0, gcd;
    exgcd(a, n, gcd, x, y);
    int d = gcd;//d為gcd(a, n)
    if(b % d)return false;//無解
    x0 = x * (b/d) % n;//特解
    for(int i = 0; i < d; ++i)cout<<x0 + i * (n / d)<<endl;//所有的解
    return true;
}

• 用歐幾里德演算法求乘法逆元：

什麼叫乘法逆元？

ax ≡ 1 (mod p)

這裡我們稱x為a關於p的逆元

上述式子可以還原為：ax + py = 1,根據上面講的擴充套件歐幾里得演算法，我們知道gcd(a, p) = 1時才有解，通過擴充套件歐幾里得演算法得到的解x₀,利用x₀%p就可以得到最小解

為什麼呢？

因為通解為x = x₀ + p * k

那麼，也就是說， a 關於 p 的逆元是一個關於 p 同餘的，那麼根據最小整數原理，一定存在一個最小的正整數，它是 a 關於p 的逆元，而最小的肯定是在（0 , p）之間的，而且只有一個，這就好解釋了。

但是，由於問題的特殊性，有時候我們得到的特解 x₀ 是一個負數，還有的時候我們的 p 也是一個負數這怎麼辦？

當 p 是負數的時候，我們取 p 的絕對值就行了，當 x₀ 是負數的時候，x₀% p 的結果仍然是負數（在計算機計算的結果上是這樣的，雖然定義的時候不是這樣的），這時候，我們仍然讓 x₀ 對abs(p) 取模，然後結果再加上abs(p) 就行了，於是，我們不難寫出下面的程式碼求解一個數 a 對於另一個數 p 的乘法逆元：

int a, b, x, y, gcd, ans, p;

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

inline int cal(int a, int p){
    exgcd(a, p, gcd, x, y);
    if(1 % gcd != 0)return -1;
    x = x * 1 / gcd;
    p = abs(p);
    ans = x % p;
    if(ans <= 0)ans += p;
    return ans;
}