扩展欧几里得算法求逆元_从辗转相除法到求逆元,数论算法初体验

今天是算法和数据结构专题的第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,毫无疑问这样的做法会超时。

辗转相除法

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

662fe7964dcf4c9fb4ab4514a4ef25dc.png

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

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

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

27a90241d420348deca533bcbf55bc52.png

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

569261f395bd3b173eb515ba52d5376c.png

这样我们就得到了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可以取任意整数。

我们代入算一下即可:

40d141955cfc92b5b9efc42b2d286966.png

所以我们求出了这样的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,这里的/是整除计算的意思,我们代入:

d19005e325d3abd4d084921a284d27ff.png

显然对于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,代入:

8c13056f641e53614470cfc6b049dae8.png

所以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为质数的时候,可以得到

e19045f617dccc53e7282fe6206d241e.png

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

ab0620130af82059e0108c74958487b1.png

也就是说我们求出a^{m-2}然后再对m取模就得到了a的逆元,我们使用快速幂可以很方便地求出来。但是这个只有m为质数的时候才可以使用。

总结

今天我们聊了欧几里得定理聊了辗转相除法还聊了拓展欧几里得和求解逆元,虽然这些内容单独来看并不难,合在一篇文章当中量还是不小的。这些算法底层的基础知识是数论,对于没有参加过竞赛的同学来说可能有些陌生,但是它也是算法领域一个很重要的分支。

如果喜欢本文,可以的话,请点个关注,给我一点鼓励,也方便获取更多文章。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/258740.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

[翻译] Fast Image Cache

https://github.com/path/FastImageCache Fast Image Cache is an efficient, persistent, and—above all—fast way to store and retrieve images in your iOS application. Part of any good iOS applications user experience is fast, smooth scrolling, and Fast Image …

centos 安装 MatConvNet (gpu)

1. 安装准备 matlab2017a &#xff0c;参考&#xff1a;《centos 安装matlab2017a(无root权限)》 GCC 4.8(支持c11) 键入&#xff1a;sudo yum install gcc gcc-c &#xff08;建议sudo装&#xff09; 至少CUDA 7.5&#xff0c;&#xff08;本人选择cuda8.0&#xff…

php练习 租房子

题目要求 1.封装类 <?php class DBDA {public $fuwuqi"localhost"; //服务器地址public $yonghuming"root";//用户名public $mima"";//密码 public $dbconnect;//连接对象//操作数据库的方法//$sql代表需要执行的SQL语句//$type代表SQL语…

【SHARE】WEB前端学习资料

参考资料&#xff1a;https://github.com/karlhorky/learn-to-program学习网站&#xff1a;http://www.codecademy.com/learn https://www.codeschool.com/ 制作网站&#xff1a;https://webmaker.org/zh-CN/explore JavaScript2015&#xff1a;https://esdiscuss.org/topic/ja…

python软件安装和使用方法_aws cli的安装及使用(内含python的安装方法)

liunx环境(使用bundled installer)&#xff1a;1.wget https://s3.amazonaws.com/aws-cli/awscli-bundle.zip //下载bundled installer2.unzip awscli-bundle.zip3.sudo ./awscli-bundle/install -i /usr/local/aws -b /usr/local/bin/aws如果你没有sudo权限或者是你想在当…

centos 安装boost(caffe需要)

安装 由于安装caffe&#xff0c;要求boost的版本在1.55以上&#xff0c;而服务器上的刚好是1.54,所以进行了重装。 参考&#xff1a;《CentOS 7下编译安装Boost_1_57_0 》 不过由于pycaffe需要boost.python,因此需要在./b2时修改为./b2 –stage debug 才可以。而不能去掉py…

JAVA正则表达式介绍和使用

本文引用自 http://www.cnblogs.com/android-html5/archive/2012/06/02/2533924.html 技术博客 1.Java中在某个字符串中查询某个字符或者某个子字串 Java代码 String s "Shang Hai Hong Qiao Fei Ji Chang";    String regEx "a|F"; //表示a或F Pat…

集合框架中的接口及其实现类

Collection&#xff1a;集合层次中的根接口&#xff0c;JDK没有提供这个接口直接地实现类。Set&#xff1a;不能包含重复的元素。SortedSet是一个按照升序排列元素的Set。List&#xff1a;是一个有序的集合&#xff0c;可以包含重复的元素。提供了按索引访问的方式。Map&#x…

C# 多线程 Parallel.For 和 For 谁的效率高?那么 Parallel.ForEach 和 ForEach 呢?

还是那句话&#xff1a;十年河东&#xff0c;十年河西&#xff0c;莫欺少年穷。 今天和大家探讨一个问题&#xff1a;Parallel.For 和 For 谁的效率高呢&#xff1f; 从CPU使用方面而言&#xff0c;Parallel.For 属于多线程范畴&#xff0c;可以开辟多个线程使用CPU内核&#x…

cuda、cudnn相关问题链接

1. cuda&#xff0c;cudnn安装 <caffe安装系列——安装cuda和cudnn> 2. 查看已有的cuda等版本 cuda 版本 cat /usr/local/cuda/version.txtcudnn 版本 cat /usr/local/cuda/include/cudnn.h | grep CUDNN_MAJOR -A 23. cudnn的安装&#xff0c;路径和版本问题 http://…

bigdecimal 小于等于0_图解小于 K 的两数之和

点击蓝色“五分钟学算法”关注我哟加个“星标”&#xff0c;天天中午 12:15&#xff0c;一起学算法作者 | P.yh来源 | 五分钟学算法题目描述 题目来源于 LeetCode 上第 1099 号问题&#xff1a;小于 K 的两数之和。给你一个整数数组 A 和一个整数 K&#xff0c;请在该数组中找出…

用STS创建Maven的Web项目转

右键New——>other——》Maven——》Maven Project 弹出框中点击Next&#xff0c;在Filter中写上&#xff1a;webapp. 然后在下面的框中选择org.apache.maven.archetypes&#xff0c;点击Next 在新弹出的窗口中写上Group Id和Artifact Id&#xff0c;Finish即可成功。 创建完…

img超出div width时, jQuery动态改变图片显示大小

参考&#xff1a; 1. http://blog.csdn.net/roman_yu/article/details/6641911 2. http://www.cnblogs.com/zyzlywq/archive/2012/02/23/2364292.html转载于:https://www.cnblogs.com/carlo/p/4584008.html

《TOGAF 9.1IT企业架构》什么是企业IT架构

2. 什么是企业IT架构 现在有越来越多的企业IT架构定义。在这一章&#xff0c;你会学习到一些企业IT架构的方法&#xff0c;我们会给你深入解释一种实用的方法&#xff0c;这种方法视企业架构师为CIO(译注&#xff1a;CIO首席信息官&#xff0c;是负责一个公司信息技术和系统所有…

pdf 深入理解kotlin协程_Kotlin协程实现原理:挂起与恢复

今天我们来聊聊Kotlin的协程Coroutine。如果你还没有接触过协程&#xff0c;推荐你先阅读这篇入门级文章What? 你还不知道Kotlin Coroutine?如果你已经接触过协程&#xff0c;但对协程的原理存在疑惑&#xff0c;那么在阅读本篇文章之前推荐你先阅读下面的文章&#xff0c;这…

编译py-faster-rcnn的问题汇总及解决方法

按照官网 的提示&#xff0c;我开始安装faster rcnn&#xff0c;但是出现了很多问题&#xff0c;我将其汇总了起来&#xff0c;并提出了解决办法。 先说明一下我的配置&#xff1a; python : anaconda2linux: centos 6.9 安装faster rcnn请先参考&#xff1a;《cuda8cudnn4 F…

openWRT自学---针对backfire版本的主要目录和文件的作用的分析整理

特别说明&#xff1a;要编译backfire版本&#xff0c;一定要通过svn下载:svn co svn://svn.openwrt.org/openwrt/branches/backfire&#xff0c;而不能使用http://downloads.openwrt.org/backfire/10.03/中的源码包&#xff1a;backfire_10.03_source.tar.bz2 结合文档《OpenWr…

自然语言交流系统 phxnet团队 创新实训 项目博客 (五)

3DMax方面所涉及的专业知识&#xff1a; &#xff08;1&#xff09;一下的关于3DMax中对于人物的设计和操作均需要在对3DMax基础知识熟练掌握的情况下进行的。 &#xff08;2&#xff09;骨骼架设&#xff1a;首先对导入到3DMax中的人物模型进行架设骨骼…

linux 安装python-opencv

三种方法&#xff1a; 1. pip 安装 &#xff1a; pip install opencv-python &#xff0c;最新版为opencv3安装后>>> import cv2 >>> print cv2.__version__参考&#xff1a;http://www.cnblogs.com/lclblack/p/6377710.html 2. anaconda的conda安装 ,可以指…

《你的灯亮着吗》读书笔记Ⅲ

转载于:https://www.cnblogs.com/yue3475975/p/4586220.html