扩展欧几里得算法

Extended Euclidean Algorithm

一、欧几里得算法

在整数域中，存在如下事实：

$(a,b)=(a-kb,b),k\in \mathbb{Z}$

因此在求解两整数$a,b(a>b)$的最大公约数时，令上式中$k=\lfloor \frac{a}{b}\rfloor$，有：

$(a,b)=(a-\lfloor \frac{a}{b}\rfloor b,b)=(amodb,b)$

借此便缩小了所求问题的规模.分别用$b,amodb$替代$a,b$，将这一操作重复地进行下去，直到$b=0$，此时$(a,0)=a$，其中的$a$便是最初两整数的最大公约数.

二、扩展欧几里得算法

扩展欧几里得算法用于求解下述问题：

对于给定的$a,b\in \mathbb{Z}$，找到一组对应的$x,y\in \mathbb{Z}$，使得$ax+by=(a,b)$

1.递归实现

(1)递归关系式

根据裴蜀定理 (Bézout’s lemma) ，存在整数$x,y,x^{\prime },y^{\prime }$，使得：

$\{\begin{array}{l}ax+by=(a,b)\\ b{x}^{\prime }+(a-\lfloor \frac{a}{b}\rfloor b)y^{\prime }=(b,a-\lfloor \frac{a}{b}\rfloor b)=(b,amodb)\end{array}$

又由于$(b,amodb)=(a,b)$，即两式右侧相等.得到：

$ax+by=b{x}^{\prime }+(a-\lfloor \frac{a}{b}\rfloor b)y^{\prime }$

整理得到：

$a(x-y^{\prime })+b(y-(x^{\prime }-\lfloor \frac{a}{b}\rfloor y^{\prime }))=0$

上式对于任意$a,b$均成立，因此：

$\{\begin{array}{l}x=y^{\prime }\\ y=x^{\prime }-\lfloor \frac{a}{b}\rfloor y^{\prime }\end{array}$

只要知道了$x^{\prime },y^{\prime }$，就能得到$x,y$，借此便缩小了求解问题的规模.

(2)终止条件

每次递归时，分别用$b,amodb$替代$a,b$，直到$b=0$.此时$x=1,y=0$.

代码详见Extended Euclidean Algorithm.py.

2.迭代实现

递归关系用矩阵表示：

$\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 1& -d_1\end{array}\right)\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)$

因此：

$\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 1& -d_1\end{array}\right)\left(\begin{array}{cc}0& 1\\ 1& -d_2\end{array}\right)...\left(\begin{array}{cc}0& 1\\ 1& -d_n\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)$

$d_i$是第$i$轮迭代时的$\lfloor \frac{a}{b}\rfloor$.

代码详见Extended Euclidean Algorithm.py.

3.手工计算

记$ax+by=d$，分别令$x=1,y=0;x=0y=1$：

$$\begin{gathered} 1\cdot a+0\cdot b=a(1) \\ 0\cdot a+1\cdot b=b(2) \end{gathered}$$

利用$(1)(2)$两式，计算$(1)-\lfloor \frac{a}{b}\rfloor (2)$：

$a-\lfloor \frac{a}{b}\rfloor b=amodb(3)$

类似地，利用$(2)(3)$两式，计算$(2)-\lfloor \frac{b}{amodb}\rfloor (3)$

反复执行这一操作，直到等式右侧为$(a,b)$

三、在GF( p )中求乘法逆元

在$GF(p)$中求$b$的乘法逆元$y=b^{-1}modp$，即是求解满足$px+by=1modp$的$y$.

利用上述扩展欧几里得算法求解即得.

def recursive_ext_euclid(a: int, b: int):
    # recursive
    if a < b:
        return recursive_ext_euclid(b, a)
    if b == 0:
        return 1, 0
    else:
        x, y = recursive_ext_euclid(b, a % b)
        q = a // b
        return y, (x - q * y)

def iterative_ext_euclid(a: int, b: int):
    # iterative
    if a < b:
        return iterative_ext_euclid(b, a)
    # diag(2)
    m11, m12, m21, m22 = 1, 0, 0, 1
    while True:
        if b == 0:
            return m11, m21
        q = a // b
        a, b = b, a % b
        m11, m12, m21, m22 = m12, m11 - q*m12, m22, m21 - q*m22