最简单的数论法（详解数论算法之欧几里得与求逆元）

今天是算法和数据结构专题的第22篇文章，我们一起来聊聊辗转相除法，下面我们就来聊聊关于最简单的数论法?接下来我们就一起去了解一下吧!

最简单的数论法

今天是算法和数据结构专题的第22篇文章，我们一起来聊聊辗转相除法。

辗转相除法又名欧几里得算法，是求最大公约数的一种算法，英文缩写是gcd。所以如果你在大牛的代码或者是书上看到gcd，要注意，这不是某某党，而是指的辗转相除法。

在介绍这个算法之前，我们先来看下最大公约数问题。

暴力解法

这个问题应该很明确了，我们之前数学课上都有讲过。给我们纸笔让我们求都没有问题，分解因数找下共同的部分，很快就算出来了。但是用代码实现怎么做呢？

用代码实现的话，首先排除分解因数的方法。因为分解因数复杂度太高了，也很容易想明白，既然要分解因数，那么首先需要获得一定量的质数吧。有了质数之后还要遍历质数，将整数一点一点分解，显然很麻烦，还不如直接暴力了。暴力解法并不复杂，我们直接从1开始遍历，记录下来同时能够整除这两个数的最大数即可。我们暴力的范围也不大，从1到n。

很容易写出代码：

def gcd(a, b): ret = 0 for i in range(min(a, b)): if a % i == 0 and b % i == 0: ret = i return ret

这个很简单，也许你可能还会想出一些优化，比如说首先判断一下a和b之间是否有倍数关系，如果有的话直接就可以得到结果了。再比如说我们i的遍历范围其实可以不用到min(a, b)，如果a和b没有倍数关系的话min(a, b) / 2就可以了。这些都是没有问题的，但是即使加上了这些优化依然改变不了这是一个O(n)算法的本质。

比如说a是1e9，b是1e9-1，毫无疑问这样的做法会超时。

辗转相除法

接下来就轮到正主——辗转相除法出场了，这个算法在《九章算术》当中曾经出现过，叫做更相减损术。不管叫什么，原理都是一样的，它的最核心本质是下面这个式子：

这个式子就是著名的欧几里得定理，这里的r可以看成是a对b取余之后的结果，也就是说a和b的最大公约数等于b和r的最大公约数。这样我们就把a和b的gcd转移成了b和r，然后我们可以继续转移，直到这两个数之间存在倍数关系的时候就找到了答案。

在我们写代码之前，我们先来看一下这个定理的证明。

我们假设u同时整除a和b，显然这样的u一定存在，因为u至少可以是1，所以：

所以可以得到u也整除r，同样我们可以证明能够整除b和r的整数也可以整除a。我们假设v可以同时整除b和r：

这样我们就得到了v也可以整除a。也就是说a和b的每一个因子都是b和r的因子，同样b和r的每一个因子也是a和b的因子，那么可以得出a和b的最大公约数就是b和r的最大公约数。

以上就是欧几里得定理的简单证明，如果看不懂也没有关系，我们记住这个定理的内容就可以了。

接下来就是用代码实现了，我们把这个公式套进递归当中非常容易：

def gcd(a, b): if a < b: a, b = b, a if a % b == 0: return b return gcd(b, a % b)

我们首先判断了a和b的大小关系，如果a小于b的话，我们就交换它们的值，保证a大于b。如果a和b取模的结果为0，那么说明a已经是b的倍数了，显然它们之间的最大公约数就是b。

但其实我们没有必要判断a和b的大小，我们假设a小于b，那么显然a % b = a，于是会递归调用b和a % b，也就是b和a，也就是说算法会自动调整两者的顺序。这么一来，这个代码还可以进一步简化，只需要一行代码。

def gcd(a, b): return a if b == 0 else gcd(b, a % b)

所以听到有人说自己用一行代码实现了一个算法，不要觉得它在装逼，有可能他真的写了一个gcd。

拓展欧几里得

拓展欧几里得本质上就是gcd，只是在此基础上做了一定的拓展，从而来解决不定方程。不定方程就是ax by = c的方程，方程要有解充要条件是(a, b) | c，也就是说a和b的最大公约数可以整除c。

也就是说求解ax by = gcd(a, b)的解。假如说我们找到了这样一组解x0和y0，那么x0 (b / gcd) * t和y0 - (a / gcd) * t也是方程的解，这里的t可以取任意整数。

我们代入算一下即可：

所以我们求出了这样的x0和y0之后就相当于求出了无数组解，那么这个x0和y0怎么求呢，这就需要用到gcd算法了。

我们观察一下gcd算法的递归代码，可以发现算法的终止条件是a=gcd，b=0。对于这样的a和b来说，我们已经找到了一组解使得ax by=gcd，比如很明显，x=1，y=0。实际上y可以为任何值，因为b=0。

我们回到递归的上一层的a和b，假设我们已经求出了b和a%b的最大公约数，并且求出了一组解x0和y0。使得b*x0 (a%b)* y0 = gcd。那么我们能不能倒推得到a和b时候的解呢？

因为a % b = a - (a/b)*b，这里的/是整除计算的意思，我们代入：

显然对于a和b来说，它的一组解就是y0和x0 - (a/b)*b*y0，我们把这几行计算加在代码当中即可，非常简单：

def exgcd(a, b, x=1, y=0): # 当b=0的时候return if b == 0: return a, x, y # 递归调用，获取b, a%b时的gcd与通项解 gcd, x, y = exgcd(b, a%b, x, y) # 代入，得到新的通项解 x, y = y, x - a//b*y return gcd, x, y

这里我建议大家不要死记代码，都去推导一下递归的这个推导公式。这个公式搞明白了，即使代码记不住也没有关系，后面临时用到的时候再推导也可以。不然的话，即使背下来了代码也不记得什么意思，如果碰到的场景稍微变动一下，可能还是做不出来。

逆元与解逆元

拓展欧几里得算法我们理解了，但是好像看不出来它到底有什么用。一般情况下我们也碰不到让我们计算通解的情况，但其实是有用的，用的最多的一个功能就是计算逆元。

在解释逆元之前先来看一个问题，我们有两个数a和b，和一个模底数p。我们可以得到(a b) % p = (a%p b%p)%p，也可以得到 (a - b)%p = (a%p - b%p)%p。甚至还可以得到 (a*b)% p =(a%p * b%p) %p，这些都是比较明确的，但是(a / b) % p = (a % p / b % p) % p，这个式子成立吗？

最后的式子是不成立的，因为模数没有除法的传递性，我们可以很方便举出反例。比如a是20， b是10，p是4，(a/b)%p=2，而(a %p / b%p) % p = 0。

这就导致了一个问题，假如说我们在一连串计算当中，由于最终的结果特别大，我们无法存储精确的值，希望存储它关于一个模底数取模之后的结果。但是我们的计算当中又涉及除法，这个时候应该怎么办？

这个时候就需要用到逆元了，逆元也叫做数论倒数。它其实起到一个和倒数类似的效果，假设a关于模底数p的逆元是x，那么可以得到：ax = 1 (mod p)

所以我们想要算 (a / b) % p，可以先求出b的逆元假设是inv(b)，然后转化成(a%p * inv(b)%p)%p。

这个逆元显然不会从天上掉下来，需要我们设计算法去求出来，这个用来求的算法就用到拓展欧几里得，我们下面来看一下推导过程。

假设a和b互质，那么gcd(a, b) = 1，代入：

所以x是a关于b的逆元，反之可以证明y是b关于a的逆元。

这么计算是有前提的，就是a和b互质，也就是说a和b的最大公约数为1。否则的话这个计算是不成立的，也就是说a没有逆元。那么整个求解逆元的过程其实就是调用拓展欧几里得的过程，把问题说清楚花了很多笔墨，但是写成代码只有两三行：

def cal_inv(a, m): gcd, x, y = exgcd(a, m) # 如果gcd不为1，那么说明没有逆元，返回-1 return (x % m m) % m if gcd == 1 else -1

在return的时候我们对x的值进行了缩放，这是因为x有可能得到的是负数，我们把它缩放到0到m的范围当中。

逆元的求解方法除了拓展欧几里得之外，还有一种算法，就是利用费马小定理。根据费马小定理，在m为质数的时候，可以得到

等式两边同时除以a，也就是乘上a的逆元，可以得到：