模逆——拓展欧几里得 - wuuconix's blog

背景

在准备用python实现AES的时候，遇到了求伽罗华域下一个多项式的逆的问题。我发现，我不光把域的知识忘光了，别说多项式的逆了，我连如何用python实现求一个整数的逆都不知道。

我只了解手动手动推的过程，主要利用到了贝组等式和欧几里得的逆推。

手动推的过程

基于贝组等式ax+by=gcd(a,b)ax+by = gcd(a,b)ax+by=gcd(a,b)有解。

比如我们这里求a在模b意义下的逆c。

因为逆如果存在，则a和b的最大公因数一定是1。

贝组等式变成了这样：ax+by=1ax+by=1ax+by=1

这时候，我们如果等式两边同时模b。

就变成了ax=1 mod bax = 1\ mod\ bax=1 mod b。

那x不就是我们要求的a在模b意义下的逆了嘛！

所以转化一下思路，我们只要求出了贝组等式中的x，那么我们就求出了a的逆！

所以这个问题实际上是一个二元一次方程式求解的问题。

在大二下半学期的数学基础课程中，我做为530知名做题家之一，已经牢牢记住了这个x是怎么求的。

下面拿求 7 在 模 26 下的逆做例子。

首先我们写出辗转相除法的所有推导过程。

26 = 3 * 7 + 5
7  = 1 * 5 + 2
5  = 2 * 2 + 1

然后从最后一个式子5 = 2 * 2 + 1出发，把1移到一边，另外的移到另一边。

1 = 5 - 2 * 2

之后把2 用移项后的第二个式子2 = 7 - 1 * 5进行代替

1 = 5 - 2 * (7 - 1 * 5)

整理一下。

1 = -2 * 7 + 3 * 5

接下来，类似的，再把5换掉。

1 = -2 * 7 + 3 * (26 - 3 * 7)

整理

1 = 3 * 26 - 11 * 7

是不是很熟悉了？两边再换一下位置！

3 * 26 - 11 * 7 = 1

这不就是贝组等式嘛！而且 7 前面的 系数我们也知道了，是 -11，也就是15

所以 15 就是 7 的逆了。

思考

上面手动推导的过程实际上是利用已有的一些式子（求gcd写的一些式子）来进行变量替换，从1 * 5 - 2 * 2 = 1最终变为了3 * 26 - 11 * 7 = 1。

如果用算法完全模拟这个过程，我们首先要生成并保存这些式子，然后进行一些替换，感觉实现起来充满了繁琐。

所以我们需要简化这个过程，最好能找出每一组的关系，得到某个简单的关系来进行迭代。就像gcd(a,b)=gcd(b,a%b)gcd(a,b)=gcd(b,a\% b)gcd(a,b)=gcd(b,a%b)这种漂亮的式子一样。

这种关系，拓展欧几里得算法Extended Euclidean algorithm给出来了。这里 主要参考了知乎大佬不抱怨的世界的回答。

这里写一下我理解的推导出两个式子之间关系的过程。

我们仔细观察手动推出的中间式子，去掉那些因为替换变量产生的中间式子，可以得到以下一些重要的式子。

1 =  1 * 1 +  0  * 0
1 =  0 * 2 +  1  * 1
1 =  1 * 5  - 2  * 2
1 = -2 * 7  + 3  * 5
1 =  3 * 26 - 11 * 7

这些式子分别是1与0、2与1、5与2、7与5的贝组等式和我们要求的26和7的贝组等式。

所以我们实际上是从一个贝组等式出发，经过某种规律得到了上一个贝组等式，逐步推导得到了我们需要的贝组等式，从而求出逆。

我们的注意点应该放在前面的系数上，也就是x和y上，因为我们最后求得逆就是这个系数。

这里列两个贝组等式。

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

bx+(a%b)y=gcd(b,a%b)bx + (a \% b)y = gcd(b, a\%b)bx+(a%b)y=gcd(b,a%b)

假设一个式子得解为(x1,y1)(x_1, y_1)(x1​,y1​)。第二个式子得解为(x2,y2)(x_2,y_2)(x2​,y2​)，代入等式。

ax1+by1=gcd(a,b)ax_1+by_1=gcd(a,b)ax1​+by1​=gcd(a,b)

bx2+(a%b)y2=gcd(b,a%b)bx_2 + (a \% b)y_2 = gcd(b, a\%b)bx2​+(a%b)y2​=gcd(b,a%b)

而两个等式右边是相等得，就是辗转相除法嘛2333。

所以 ax1+by1=bx2+(a%b)y2ax_1 + by_1 = bx_2 + (a\%b)y_2ax1​+by1​=bx2​+(a%b)y2​

将a%ba\%ba%b等价转化为a−(a//b)∗ba - (a //b) * ba−(a//b)∗b，再代入。

⇒ax1+by1=bx2+(a−(a//b)∗b)y2\Rightarrow ax_1 + by_1 = bx_2 + (a - (a//b) * b)y_2⇒ax1​+by1​=bx2​+(a−(a//b)∗b)y2​

⇒ax1+by1=bx2+ay2−(a//b)∗by2\Rightarrow ax_1 + by_1 = bx_2 + ay_2 - (a//b)*by_2⇒ax1​+by1​=bx2​+ay2​−(a//b)∗by2​

⇒ax1+by1=ay2+b(x2−(a//b)y2)\Rightarrow ax_1 + by_1 = ay_2 + b(x_2 - (a//b)y_2)⇒ax1​+by1​=ay2​+b(x2​−(a//b)y2​)

⇒x1=y2 and y1=(x2−(a//b)y2)\Rightarrow x_1 = y_2\ and\ y_1 = (x_2 - (a//b)y_2)⇒x1​=y2​ and y1​=(x2​−(a//b)y2​)

ok，至此我们已经得到了上下两个贝组等式之间系数的关系了。

在这里(x1,y1)(x_1,y_1)(x1​,y1​)对应的是一开始(a,b)(a,b)(a,b)最大的贝组等式的解。

我们可以根据(b,a%b)(b, a \% b)(b,a%b)贝组等式解(x2,y2)(x_2,y_2)(x2​,y2​)来求出(x1,y1)(x_1, y_1)(x1​,y1​)。

而要知道(x2,y2)(x_2,y_2)(x2​,y2​)，我们就得知道(x3,y3)(x_3,y_3)(x3​,y3​)，就这样深入下去，我们需要一个终点，即一组确定的解。

然后从这组确定的解再原路返回，在回去的过程中将算出每一个解，知道开始的(x1,y1)(x_1,y_1)(x1​,y1​)。

这个过程实际上就是递归的过程。（现在想来递归是不是就是一个栈呢？一开始依次入栈，到达某种条件，像俄罗斯方块拼成一行了，就依次出栈消除

而那个终点实际上也很容易取。

我们用辗转相除法的时候的等式 gcd(a,b)=gcd(b,a%b)gcd(a, b) = gcd(b, a\%b)gcd(a,b)=gcd(b,a%b)。最后结束的情况就是b为零，那时候a就是它们的最大公因数。

所以我们考虑gcd(a,0)gcd(a, 0)gcd(a,0)的贝组等式。

很显然，我们能够构造出这样的等式 1∗a+0∗0=a1*a + 0 * 0 = a1∗a+0∗0=a。

所以这里的解就是(1,0)(1, 0)(1,0)。这便是终点。

代码实现

这里给出百度百科的exgcd的python代码实现。

def ext_euclid(a, b):     
    if b == 0:         
        return 1, 0, a     
    else:         
        x, y, q = ext_euclid(b, a % b)      
        x, y = y, (x - (a // b) * y)         
        return x, y, q

它的返回值有三个值，分别为x, y和 gcd。

也不难理解哈哈，它叫拓展欧几里得，普通的欧几里得的功能（得到最大公因子）当然也不能落下啦！

然后如何求模逆呢？

很简单，首先得判断它得gcd是不是1，只有a和b互素情况下，a才有逆。

然后x就是a的逆啦！

def exgcd(a, b):
	if (b == 0):
		return 1, 0, a
	else:
		x, y, gcd= exgcd(b, a % b)
		x, y = y, (x - (a // b) * y)
		return x, y, gcd

print(exgcd(7, 26))

还可以使用libnum这个强大的库来直接求模逆。这里把它的相关函数copy过来。

def xgcd(a, b):
    """
    Extented Euclid GCD algorithm.
    Return (x, y, g) : a * x + b * y = gcd(a, b) = g.
    """
    if a == 0:
        return 0, 1, b
    if b == 0:
        return 1, 0, a

    px, ppx = 0, 1
    py, ppy = 1, 0

    while b:
        q = a // b
        a, b = b, a % b
        x = ppx - q * px
        y = ppy - q * py
        ppx, px = px, x
        ppy, py = py, y

    return ppx, ppy, a
    
def invmod(a, n):
    """
    Return 1 / a (mod n).
    @a and @n must be co-primes.
    """
    if n < 2:
        raise ValueError("modulus must be greater than 1")

    x, y, g = xgcd(a, n)

    if g != 1:
        raise ValueError("no invmod for given @a and @n")
    else:
        return x % n

实际上就是拓展欧几里得，它用的拓展欧几里得是迭代版的，下次再去学习一下。

然后它判断了一下a和b的最大公因子是不是1，不是的话报错，提示无法求逆。

​

战术总结

今天总算是把拓展欧几里得的原理弄懂了，还经过7att1ce的推荐知道了libnum这个强大库的存在。总算是迈出了代码实现AES的第一步啦！