欧几里德算法

淩亂°似流年 2022-05-08 07:00 227阅读 0赞

原文转载：博主——AdemJensen[https://blog.csdn.net/rentenglong2012/article/details/68944518][https_blog.csdn.net_rentenglong2012_article_details_68944518]

写得很好

>  
> 
>  *  欧几里德算法
>     
>      *  定义
>      *  求法推演
>      *  编程实现
>         
>          *  CPP代码递归
>          *  CPP代码迭代
>  *  扩展欧几里德算法
>     
>      *  定义
>      *  算法推演
>         
>          *  推演方式1
>             
>              *  考虑情况1
>              *  考虑情况2
>              *  考虑情况3
>          *  推演方式2
>      *  编程实现
>         
>          *  CPP定义
>          *  CPP精简
>      *  关于最小正整数解的问题
>  *  应用
>     
>      *  求解不定方程二元一次方程的最小正整数解
>         
>          *  内容
>          *  算法
>          *  证明正确性
>          *  编程实现
>             
>              *  CPP代码
>          *  关于最小正整数解的问题
>      *  求乘法逆元
>         
>          *  定义
>          *  求法
>          *  编程实现
>             
>              *  CPP代码
>      *  解线性同余方程
>         
>          *  定义
>          *  解法
>          *  编程实现
>             
>              *  CPP代码
> 
>

欧几里德算法。。。用来求最大公约数嘛。。。

咳，至于教练说欧几里德算法“很简单的东西”。。。嘛。。。理解起来确实有困难呢Orz。。。

不过好像最近（明明是一直）考的不少。。。写个博客整理一下。。。

# **欧几里德算法** #

首先是欧几里德算法。。。

## **定义** ##

我们设（射）函数 gcd(a,b)

为求解两个**自然数**的最大公因数的函数。例如：gcd(25,35)=5

那么。。。咋求（球）呢。。。

## **求法推演** ##

首先考虑一下：

对于任意两个正整数 a,b

，都有：

a=kb+r(k,r∈N)

所以有：r=a%b

（在这里，%指的是**取余运算**）

然后我们假设 c

是 a 和 b

的最大公约数，即c=gcd(a,b)

然后，我们就能得到：c|a      c|b

（在这里当然包括除了在这里的任何地方， a|b表示 b 能够整除 a）

然后又因为上面那个式子，有：

> r=a−kb

所以有：

> c|r

整合一下上面的式子，我们可以得到：

> c=gcd(b,r)

即

> gcd(a,b)=gcd(b,a%b)

## **编程实现** ##

然后我们就可以编程咯。。。递归的方法是再好不过了。。。

但是还要考虑一个问题：什么时候停止递归。。。

我们由乘除法的性质，可以得到：0除以任何一个非0的数都得0，即：

> 0÷x=0
> 
> (x∈(−∞,0)∪(0,+∞))

所以0与一个数的最大公约数是这个数本身，即：

> gcd(a,0)=gcd(0,a)=a

那就太好办了，反正 a%b

一定是在不断变小，并且总会有它等于0的时候，那么我们不妨将它设置为递归终点。

### **CPP代码（递归）** ###

int gcd(int a,int b)
    {
        if(b==0) return a;
        return gcd(b,a%b);
    }

这种写法就是**辗转相除法**。

当然为了防止某些情况爆栈（比如说高精度运算什么的），还可以不用递归来写。。。

### **CPP代码（迭代）** ###

int gcd(int a,int b)
    {
        while(b)
        {
            a=a%b;
            swap(a,b);
        }
        return a;
    }

当然本质上这两种计算方式都是一模一样的23333

# **扩展欧几里德算法** #

然后就是邪恶的扩展欧几里得算法。。。

## **定义** ##

扩展欧几里德算法旨在解决这样一个问题：对于以下二元一次方程：

> ax+by=gcd(a,b)
> 
> (a,b∈N)

求解出其中的一组关于 x和 y的解。啊啊啊注意啊这里的一组并不是最小的那一组！！！

然后我们来看看。。。

## **算法推演** ##

还是从最基本的谈起。。。

在这里提供两种推演方式，一个是自己弄得，另外一个好像是网上比较常见的。。。

### **推演方式1：** ###

**考虑情况1：**

a=b

的情况

我们可以得到：

> ax+ay=a

所以：

> x+y=1

所以：

> x=1
> 
> y=0

或者

> x=0
> 
> y=1

**考虑情况2：**

b=0

的情况

当b=0

时，必有：

> ax+0y=gcd(0,a)

即：

> ax=a

所以：

> x=1
> 
> y=−1,0,1,2,3,4,......,+∞

这样，无论y

取何值，都是一组有效的整数解。所以为了方便，我们取y=0，即：

> x=1

y=0

**考虑情况3：**

a,b≠0

的情况

由那个二元一次方程，我们可以推导出：

> ax+by=bx′+(a%b)y′

（假设 x′,y′是另外一组解）

把中间那个“a%b

”拆开（注意那个除号是整除）：

> ax+by=bx′+(a−a/b⋅b)y′

然后稍微整理一下：

> ax+by=bx′+ay′+(a/b⋅b)y′

然后再稍微整理一下，就得到：

> ax+by=ay′+b(x′−a/b⋅y′)

然后就可以得到：

> x=y′
> 
> y=x′−a/b⋅y

然后继续往下算，递归计算，直到遇到前面说的两种情况就好了。

### **推演方式2：** ###

或者说你看不惯我的证明放肆方式，那么下面还有一个，是我同学当时写的一个手记，我用马克当写了一下：

以下转载自同学的手记

已知原方程 ax+by=gcd(a,b)

>   
> 当a=b时 :  
> x=1,y=0 或者 x=0,y=1  
> 当a>b时 :  
> 当b=0时 : gcd(a,b)=a x=1,y=0  
> a⋅b≠0时：  
> 设以下方程组：  
> a1x1+b1y1=gcd(a1,b1)  
> a2x2+b2y2=gcd(a2,b2)  
> a3x3+b3y3=gcd(a3,b3)  
> ······  
> aixi+biyi=gcd(ai,bi)
> 
>  
> 
>  
> 
> 那么a1=a
> 
> b1=b (x1,y1)即一组解  
> 因为gcd(ai,bi)=gcd(bi,ai%bi)  
> 所以ai=bi−1 ，bi=ai−1%bi−1  
> 所以gcd(ai,bi)=gcd(bi,ai%bi)=gcd(ai+1,bi+1)
> 
>  
> 
>  
> 
> 又因为gcd(ai,bi)=gcd(a2,b2)=⋅⋅⋅⋅⋅⋅=gcd(ai,bi)
> 
>   
> 所以a1x1=b1y1=a2x2+b2y2  
> 所以a2=b1 b2=a1%b2=a1−(a1/b1)b1  
> 所以a1x1+b1y1  
> =b1x2+(a1−(a1/b1)b1)y2  
> =a1y2+b1(x2−(a1/b1)y2)  
> 所以x1=x2 ，y1=x2−(a/b)y2  
> 依次类推：  
> 求得 (a/gcd)x+(b/gcd)y=1  
> 所以(a/gcd)(x+(b/gcd))+(b/gcd)(y−(a/gcd))=1  
> x的周期是b/gcd

## **编程实现** ##

按照上面所述的a,b≠0

的情况进行，直到遇到前面说的两种递归终点就好了。

### **CPP（定义）** ###

int exgcd(int a,int b,int &x,int &y)
    {
        if(b==0)
        {
            x=1;
            y=0;
            return a;
        }
        int nx,ny;
        int q=exgcd(b,a%b,nx,ny);
        x=ny;
        y=nx-a/b*ny;
        return q;
    }

#include <iostream>
     
    using namespace std;
     
    int exgcd(int a, int b, int &x, int &y)     ///“&x”，对x的引用，在函数里面也可更改x值，而不仅仅是更改原先值的副本
    {
      if(b == 0)
      {
        x = 1;
        y = 0;
        return a;
      }
      int r=exgcd(b, a % b, x, y);
      int t = x;
      x = y;
      y = t - a / b * y;
      return r;                   ///最终，r依旧是a, b的最大公约数
    }
     
    int main()
    {
      int x,y; 
      cout<<exgcd(40, 30, x, y)<<endl;
      cout<<x<<endl<<y;
    }

在这里，函数的返回值即是a

和b的gcd

。

使用地址引用的方式传入x

和y两个变量，执行完函数后，原先的x和y

变量就已经被存入了一组解。

### **CPP（精简）** ###

int exgcd(int a,int b,int &x,int &y)
    {
        if(b==0)
        {
            x=1;
            y=0;
            return a;
        }
        int q=exgcd(b,a%b,y,x);
        y-=a/b*x;
        return q;
    }

关于这个程序对于非递归终点的情况处理为什么如此精神焕发活力四射姬情满满奇怪，在这里做一下解释哈

为什么再下一层递归中要**将本层递归的y作为下一层递归的x，将本层递归的x作为下一层递归的y**呢。。。

别着急，先接着往下看。

你会注意到，将本层的x作为下一层的y，正好满足了上面的x=y′这一递推式。

然后，y被减去一个a/b∗x，即是**下面的递归结束后，返回的x′****值被减去一个a/b∗y′**。

即是y=x′−a/b⋅y′，也满足了递推关系式。

所以总结一下：我们直接将本层递归的x和y变量交换传入下一层，达到了以下目的：

*  省略掉了中间变量占用的时间和空间复杂度
 *  最大化利用了函数的地址引用

然而更推荐第一种写法。。。并不怎么容易写错。。。第二种看似精简，其实很容易出错的。。。（其实我本人很喜欢第二种的说。。。但是没办法，在考场上确实很容易写错）

## **关于最小正整数解的问题** ##

先普及一下，最小正整数解指的是**其中一个数是正整数**就好啦，不用两个都是。。。

我们来想一下如何将x

值变为最小正整数。

在原有方程左右两侧同除一个gcd(a,b)

，得（以下将gcd(a,b)简化为gcd

）：

> (a/gcd)x+(b/gcd)y=1

然后在左右两侧同加一个ab/gcd2，并移项，合并同类项，得到：

> (a/gcd)(x+(b/gcd))+(b/gcd)(y−(a/gcd))=1

然后我们得到：每当x的值增加b/gcd时，都需要y的值对应减小a/gcd来维持等式。

所以，满足的解集：

> x,y  
> x+(b/gcd),y−(a/gcd)  
> x+2(b/gcd),y−2(a/gcd)  
> x+3(b/gcd),y−3(a/gcd)

那么我们称，这个**方程的解**具有**周期性**，x的周期为b/gcd，y的周期为a/gcd。

也就是说，满足( x+n(b/gcd),y−n(a/gcd))就是有效的解。

然后我们对得出来的x值进行一下这样的操作：

> 1.对b/gcd取余  
> 2.加上b/gcd  
> 3.对b/gcd取余

因为x求出的是任意一组可能。

这样得出来的x值一定是最小整数解。然后我们再推导一下怎么求y就行了。

为啥？因为如果我们对b/gcd取余，如果是负数，就会得出负的结果，不满足要求，所以要再加上b/gcd，使它变成正数。然后如果一开始是正数，那么一加就不是最小的整数解了，所以再.对b/gcd取余。这样无论是正数还是负数，都能得到最小正整数解。

如果想让y值是最小正整数解，只要对y进行相同的操作，但是注意是对a/gcd取余和相加就行了。即：对y

进行如下操作：

> 1.对a/gcd取余  
> 2.加上a/gcd  
> 3.对a/gcd取余

然后再用相同的方法把x

算出来就行了。

# **应用** #

## **求解不定方程（二元一次方程）的最小正整数解** ##

这玩意绝逼是个坑。。。

### **内容** ###

求解形如：

> ax+by=c

形式的**有整数解的**方程。

这个东西。。。已经没有什么好说的了。。。一开始以为多么高端的一个东西。。。结果发现它只能解出**有整数解的二元一次方程**。。。

### **算法** ###

算法。。。很简单。。。

一开始的时候，扩展欧几里德算法算出来的不是形如：

> ax+by=gcd(a,b)

的方程嘛。。。

解方程的方法也超级简单，你只需要先算出c是gcd(a,b)的几倍，然后将扩展欧几里德求出来的x和y分别乘以这个求出来的倍数。。。就是答案了。。。（喵了个咪的就是这个东西害我以为欧几里得算法是螺旋上天法力无边结果发现原来只是这样一个东西）

### **证明正确性** ###

然后我们要注意了。。。由于要解出的是**整数解**。。。所以就**必须**要满足：

> gcd(a,b)|c

否则方程是**木有整数解**滴。。。

道洗呆？

先简化这个问题。首先来讲，从上面的解法中可以看出，右边的c

变化，在有整数解的情况下，x和y

只是对应的放大。

所以问题就成了为什么只有满足：

> n⋅ax+n⋅by=n⋅gcd(a,b)

才有正整数解。。。

分情况讨论一下：

*  a和b
 *  相等的情况：

这种情况下，方程的形式就变成了：

> n⋅a(x+y)=n⋅a

的形式。所以如果右边一旦不是n⋅a

的形式，就无法保证方程有整数解。

*  a和b不相等的情况：
 *  a和b
 *  的最大公因数等于1的情况：这种情况不管怎么样都有整数解，想想为啥。
 *  a和b这种情况下，我们设a<b。反正如果不满足的话互换他俩就是咯。然后我们想：如果a和b的最大公约数不为1，那么一定满足：

> a=n⋅gcd(a,b)
> 
> b=m⋅gcd(a,b）

所以在有整数解的方程ax+by=c

中，一定满足：

> n⋅x⋅gcd(a,b)+m⋅y⋅gcd(a,b)=q⋅gcd(a,b)

即：

> q=nx+my

即：c的值一定是最大公约数的倍数，证毕。

### **编程实现** ###

原理很简单，遵循算法按照定义来就好了。。。

这里就不再重复写exgcd的代码了，直接引用上面给出的就好啦。

**CPP代码**

bool eqa_solve(int &x,int &y,int a,int b,int c)
    {
        int q=exgcd(a,b,x,y);
        if(c%q!=0) return false;
        int tms=c/q;
        x*=tms;
        y*=tms;
        return true;
    }

解释用法：将要存储结果的变量填写在x和y的位置，然后按照方程ax+by=c的形式填写a、b、c

三个位置，返回值为方程是否有整数解。

### **关于最小正整数解的问题** ###

跟扩展欧几里德算法还是这玩意：

> (a/gcd)(x+(b/gcd))+(b/gcd)(y−(a/gcd))=1

这个**方程的解**同样具有**周期性**，x的周期为b/gcd，y的周期为a/gcd。也就是说，满足( x+n(b/gcd),y−n(a/gcd))就是有效的解。

然后我们对得出来的x值进行一下这样的操作：

> 1.对b/gcd取余  
> 2.加上b/gcd  
> 3.对b/gcd取余

这样得出来的x值一定是最小整数解。然后我们再推导一下怎么求y就行了。

如果想让y值是最小正整数解，只要对y进行相同的操作，但是注意是对a/gcd取余和相加就行了。即：对y进行如下操作：

> 1.对a/gcd取余  
> 2.加上a/gcd  
> 3.对a/gcd取余

然后再用相同的方法把x算出来就行了。