詳解歐幾里得輾轉相除法求逆元及代碼實現

雖然網上有很多博客詳細的介紹了歐幾里得輾轉相除法求逆元，但是在當初我沒有這些數學基礎的時候，着實看着很蛋疼，所以特意寫一篇文章來彌補當時菜雞的心靈。

0x00 廣義歐幾里得除法計算最大公因數

這個算法是用來計算兩個整數\(a,b\)的最大公因數。
具體有以下的過程:

設\(a,b\)是任意兩個正整數，記\(r_{-2}=a,r_{-1}=b\)。反復運用歐幾里得除法有:

\[r_{-2}=q_0 · r_{-1}+ r_{0}, 0<r_0<r_{-1}, \]

\[r_{-1}=q_1 · r_{0}+ r_{1}, 0<r_1<r_{0}, \]

\[r_{0}=q_2 · r_{1}+ r_{2}, 0<r_0<r_{-1}, \]

\[. \]

\[r_{n-3}=q_{n-1} · r_{n-2}+ r_{n-1}, 0<r_{n-1}<r_{n-2}, \]

\[r_{n-2}=q_{n} · r_{n-1}+ r_{n}, 0<r_{n-1}<r_{n-2}, \]

\[r_{n-1}=q_{n+1} · r_{n}+ r_{n+1}, r_{n+1}=0 \]

經過有限步驟，比然存在\(n\)使得\(r_{n+1}=0\)，這是因為

\[0=r_{n+1}<r_{n}<r_{n-1}<···<r_{1}<r_{0}<r_{-1}=b, \]

且\(b\)是有限正整數.

其中可以證明的是\(n\le 5log b\) 具體證明不做贅述，詳情參見《信息安全數學基礎 第二版》第23面。

對於上面的算法，我們舉例進行說明:

例1 設\(a=1859,b=1573\)，計算(\(a,b\))

解:
運用廣義歐幾里得除法，有

\[1859=1·1573+286 \]

\[1573=5·286+143 \]

\[286=2·143+0 \]

由此可得，1859和1573的最大公約數為143

以下是基於此算法的C++代碼

#include <iostream>
using namespace std;
int GetGreatestCommonDivisor(int a, int b)
{
    a = abs(a);
    b = abs(b);
    if (a < b)
    {
        return GetGreatestCommonDivisor(b, a);
    }
    int temp;
    int remainder;
    int i = 0;

    for (; i * b <= a; i++)
    {
    }
    remainder = a - b * (i - 1);
    while (remainder != 0)
    {
        i = 0;
        for (; i * remainder <= b; i++)
        {
        }
        temp = b;
        b = remainder;
        remainder = temp - remainder * (i - 1);
    }
    ;
    return b;
}
int main()
{
    int a, b;
    while (true)
    {
        cin >> a >> b;
        int gcd = GetGreatestCommonDivisor(a, b);
        cout << gcd << endl;
    }
    return 0;
}

0x01 貝祖等式

從廣義歐幾里得除法的表示中，逐次消去\(r_{n-1},r_{n-2},...,r_2,r_1,r_0\)，可找到整數\(s,t\)使得

\[s·a+t·b=(a,b) \]

以上的等式就叫做Bézout(貝祖)等式.

這里再舉個例子，首先引用上述對於例1的解答。

例1 設\(a=1859,b=1573\)，計算(\(a,b\))

解:
運用廣義歐幾里得除法，有

\[1859=1·1573+286 \]

\[1573=5·286+143 \]

\[286=2·143+0 \]

由此可得，1859和1573的最大公約數為143

根據此例,我們有

\(143\)	\(=\)	\((-5)·286\)	\(+\)	\(1573\)
	\(=\)	\((-5)·((-1)·1573+1859)\)	\(+\)	\(1573\)
	\(=\)	\(6·1573\)	\(+\)	\((-5)·1859\)

綜上，有\(s=-5,t=6\)，滿足\(s·a+t·b=(a·b)\)

將以上的計算方法，通過代碼的形式給出來(注:這段代碼只是基於以上所述的基礎來寫的)

#include <iostream>
#include <vector>
using namespace std;

vector<int> GetGreatestCommonDivisor(int a, int b)
{

    a = abs(a);
    b = abs(b);
    if (a < b)
    {
        return GetGreatestCommonDivisor(b, a);
    }
    vector<int> result;
    result.push_back(a);

    int temp;
    int remainder;
    int i = 0;

    for (; i * b <= a; i++)
    {
    }
    remainder = a - b * (i - 1);
    while (remainder != 0)
    {
        i = 0;
        for (; i * remainder <= b; i++)
        {
        }
        temp = b;
        result.push_back(b);
        b = remainder;
        remainder = temp - remainder * (i - 1);
    };
    result.push_back(b);
    return result;
}

//vector內的數據格式
//第一位是a,b的最大公約數
//剩余四位分別對應
//sa+tb


vector<int> bezout(int a, int b)
{
    vector<int> gcdProcess = GetGreatestCommonDivisor(a, b);
    vector<int> result;
    //此處對應輸入了兩個相等的數的情況
    if (gcdProcess.size() == 1)
    {
        result.push_back(gcdProcess[0]);
        result.push_back(gcdProcess[0]);
        result.push_back(1);
        result.push_back(0);
        result.push_back(0);
        return result;
    }

    int a_b = gcdProcess[gcdProcess.size() - 1];
    result.push_back(a_b);
    result.push_back(0);
    result.push_back(0);
    result.push_back(0);
    result.push_back(0);
    for (int i = gcdProcess.size() - 2; i >= 1; i--)
    {
        if (result[2] == 0)
        {
            //適用於第一次
            result[2] = gcdProcess[i];
            result[3] = 1;
            result[4] = gcdProcess[i - 1];
            int j = 0;
            for (; result[4] + j * result[2] != a_b; j--)
            {
            }
            result[1] = j;
        }
        else
        {
            int outside = result[1];
            int inside = 0;
            for (; inside * gcdProcess[i] + gcdProcess[i - 1] != result[2]; inside--)
            {
            }
            result[1] = result[3] + outside * inside;
            result[2] = gcdProcess[i];
            result[3] = outside * 1;
            result[4] = gcdProcess[i - 1];
        }
    }

    return result;
}

int main()
{
    vector<int> result = bezout(169, 121);
    cout << result[0] << " = (" << result[1] << ") * ( " << result[2] << " ) + (" << result[3] << ") * (" << result[4] <<" )"<< endl;
    return 0;
}

其中，貝祖等式是可以被證明的，限於水平和篇幅，此處不作證明。根據其證明的過程，可以通過以下的方式,具體計算整數\(s,t\)使得

\[sa+ta=(a,b) \]

首先，令

\[r_{-2}=a,r_{-1}=b \]

\[s_{-2}=1,s_{-1}=0 \]

\[t_{-2}=0,t_{-1}=1 \]

(1) 如果\(r_{-1}=0\),則令

\[s=s_{-2},t=t_{-2} \]

否則，計算

\[q_0=[\frac{r_{-2}}{r_{-1}}],r_0=(-q_0)r_{-1}+r_{-2} \]

(2) 如果\(r_0=0\),則令

\[s=s_{-1},t=t_{-1} \]

否則，計算

\[s_0=(-q_0)s_{-1}+s_{-2},t_0=(-q_0)t_{-1}+t_{-2} \]

以及

\[q_1=[\frac{r_{-1}}{r_{0}}],r_1=(-q_1)r_{0}+r_{-1} \]

(3) 如果\(r_1=0\),則令

\[s=s_{0},t=t_{0} \]

否則，計算

\[s_1=(-q_1)s_{0}+s_{-1},t_0=(-q_1)t_{0}+t_{-1} \]

以及

\[q_2=[\frac{r_{0}}{r_{1}}],r_1=(-q_2)r_{1}+r_{0} \]

\[. \]

\[. \]

\[. \]

(\(j+2\)) 如果\(r_j=0\),則令

\[s=s_{j-1},t=t_{j-1} \]

否則，計算

\[s_j=(-q_j)s_{j-1}+s_{j-2},t_j=(-q_j)t_{j-1}+j_{j-2} \]

以及

\[q_{j+1}=[\frac{r_{j-1}}{r_{j}}],r_{j+1}=(-q_{j+1})r_{j}+r_{j-1} \]

最后一定有\(r_{n+1}=0\),這時，令

\[s=s_n,t=t_n \]

總之，可以找到整數\(s,t\)，使得

\[sa+tb=r_n=(a,b) \]

上述的表述過程可以用以下的表來表示

\(j\)	\(s_j\)	\(t_j\)	\(q_{j+1}\)	\(r_{j+1}\)
-3				\(a\)
-2	1	0		\(b\)
-1	0	1	\(q_0\)	\(r_0\)
0	\(s_0\)	\(t_0\)	\(q_1\)	\(r_1\)
1	\(s_1\)	\(t_1\)	\(q_2\)	\(r_2\)
...	...	....	...	...
n-2	\(s_{n-2}\)	\(t_{n-2}\)	\(q_{n-1}\)	\(r_{n-1}\)
n-2	\(s_{n-1}\)	\(t_{n-1}\)	\(q_{n}\)	\(r_{n}\)
n	\(s_{n}\)	\(t_{n}\)	\(q_{n+1}\)	\(r_{n+1}=0\)

其中，對於\(j=0,1,2,...,n\)，有

• \(s_j=(-q_j)s_{j-1}+s_{j-2}\)
• \(t_j=(-q_j)t_{j-1}+t_{j-2}\)
• \(q_{j+1}=[\frac{r_{j-1}}{r_j}]\)
• \(r_{j+1}=(-q_{j+1})r_j+r_{j-1}\)

上述的計算過程為
\(s_j\) -> \(t_j\) -> \(q_{j+1}\) -> \(r_{j+1}\)
當\(r_{n+1}=0\)時,\(s=s_n,t=t_n\)，使得

\[sa+tb=r_n=(a,b) \]

這里接着以例1的\(a=1859,b=1573\)為例，進行一次計算的演示

\(j\)	\(s_j\)	\(t_j\)	\(q_{j+1}\)	\(r_{j+1}\)
-3				1859
-2	1	0		1573
-1	0	1	1	286
0	1	-1	5	143
1	-5	6	2	0

所以，存在\(s=-5,t=6\)使得

\[(-5)·1859+6·1573=143 \]

用代碼實現，有以下的過程

#include <iostream>

using namespace std;

void function(int a, int b)
{
    a = abs(a);
    b = abs(b);
    if (a < b)
    {
        function(b, a);
    }
    // 對應於r_{-2},r_{-1}
    int r1 = a;
    int r2 = b;
    // 對應於s_{-2},s_{-1}
    int s11 = 1;
    int s12 = 0;
    // 對應於t_{-2},t_{-1}
    int t21 = 0;
    int t22 = 1;

    // 對應於q_j
    int q = r1 / r2;

    // 幾個中間變量
    int tempS = s12;
    int tempT = t22;
    int tempR = r1;
    // 求余數
    r1 = r2;
    r2 = -q * r2 + tempR;

    while (r2 != 0)
    {
        tempS = s12;
        tempT = t22;
        s12 = (-q) * s12 + s11;
        t22 = (-q) * t22 + t21;
        s11 = tempS;
        t21 = tempT;
        q = r1 / r2;
        tempR = r1;
        r1 = r2;
        r2 = -q * r2 + tempR;
    }

    cout<<s12<<" * "<<a<<" + "<<t22<<" * "<<b<<" = "<<r1<<endl;
}

int main()
{
    function(737,635);
    return 0;
}

0x03 逆元的求解

首先明確一下逆元的概念

設\(m\)是一個正整數,a是一個整數，如果存在整數\(a'\)使得

\[a·a'≡a'·a≡1(\bmod m) \]

成立，則\(a\)叫做模\(m\)可逆元.

逆元存在的充分必要條件是\((a,m)=1\)，且存在時，其解是唯一的。
此處對於充分必要性,以及解的唯一性，在書上完整的證明，不做贅述。

因為\((a,m)=1\),所以存在\(sa+tm=1\),這里計算出來的\(s\)就是\(a\)的逆元\(a'\)。
將代碼，稍加改動就成了如下的形式:

//歐幾里得輾轉相除法求逆元
#include <iostream>
using namespace std;
int function(int a, int m)
{
    a = abs(a);
    m = abs(m);
    // 對應於r_{-2},r_{-1}
    int r1 = a;
    int r2 = m;
    // 對應於s_{-2},s_{-1}
    int s11 = 1;
    int s12 = 0;
    // 對應於t_{-2},t_{-1}
    int t21 = 0;
    int t22 = 1;

    // 對應於q_j
    int q = r1 / r2;

    // 幾個中間變量
    int tempS = s12;
    int tempT = t22;
    int tempR = r1;
    // 求余數
    r1 = r2;
    r2 = -q * r2 + tempR;

    while (r2 != 0)
    {
        tempS = s12;
        tempT = t22;
        s12 = (-q) * s12 + s11;
        t22 = (-q) * t22 + t21;
        s11 = tempS;
        t21 = tempT;

        q = r1 / r2;
        tempR = r1;
        r1 = r2;
        r2 = -q * r2 + tempR;
    }

    if (r1 == 1)
    {
        return s12 > 0 ? s12 : s12 + m;
    }
    else
    {
        return -1;
    }
}

int main()
{
    cout << function(3, 7) << endl;
    cout << function(3, 26) << endl;
    return 0;
}

0x04 總結

其實實際上，歐幾里得輾轉相除求逆元並不是很難，只是少了前面的那些基礎知識的鋪墊，所以硬看網上的文章確實頭皮發麻，當學習了基礎知識在縷一遍，並用代碼實現，確實會有很深的理解(可能吧)