模逆——拓展欧几里得 - wuuconix's blog
模逆——拓展欧几里得 - 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的第一步啦！