扩展欧几里得算法------扩展欧几里德算法

扩展欧几里得算法及其应用

一、扩展欧几里得算法

扩展欧几里得算法：对于不完全为 0 的非负整数 a，b，若gcd(a,b)表示 a，b 的最大公约数，必然存在整数对x，y ，使得 ax+by = gcd(a,b)。
算法过程：
设 a>b，当 b=0时，gcd(a,b)=a。此时满足ax+by = gcd(a,b)的一组整数解为x=1，y=0；当a*b!=0 时，
设 a*x1+b*y1=gcd(a,b)；b*x2+(a mod b)*y2=gcd(b,a mod b)；
根据欧几里得原理知 gcd(a,b)=gcd(b,a mod b)；则：a*x1+b*y1=b*x2+(a mod b)*y2；
即:a*x1+b*y1=b*x2+(a-(a/b)*b)*y2=a*y2+b*x2-(a/b)*b*y2；（a/b：表示a除以b的商）
根据恒等定理得：x1=y2; y1=x2-(a/b)*y2；
按照上面的递归思想，我们不断对gcd进行递归，一定会出现b=0的情况，此时x=1，y=2，递归结束。这样我们就得到了求解一组x1,y1 的方法
将上面的过程写成代码如下：函数exgcd的代码：C/C++代码
①
可运行代码：

#include<stdio.h>
#include <iostream>
using namespace std;
/*x,y,r定义为全局变量，x，y的结果就是
方程ax+by=gcd(a,b)的一组解*/
int x,y,r;
//递归实现求一组x，y；函数exgcd无返还值
void exgcd(int a,int b)
{
if(b==0)
{
x=1;
y=0;
r=a;
}
else
{
exgcd(b,a%b);
int temp=x;
x=y;
y=temp-(a/b)*y;
}
}
int main()
{
int a,b;
//cin>>a>>b;
scanf("%d%d",&a,&b);
exgcd(a,b);
//cout<<q<<"=("<<x<<")*"<<a<<"+("<<y<<")*"<<b<<endl;
printf("%d=(%d)*%d+(%d)*%d\n",r,x,a,y,b);
return 0;
}

注意上面的代码是将x，y，q定义成了全局变量，其中x，y表示方程的一组整数解，q表示a，b的最大公约数。
核心代码：

int x,y,r;
void exgcd(int a,int b)
{
if(b==0)
{
x=1;
y=0;
r=a;
}
else
{
exgcd(b,a%b);
int temp=x;
x=y;
y=temp-(a/b)*y;
}
}

②

可运行代码一：递归算法C/C++

#include<stdio.h>
#include <iostream>
using namespace std;
/*递归，a,b代表方程ax+by=gcd(a,b)的a,b
x,y代表方程的一组解，其中x，y并未定义为
全局变量，而是将x，y的地址传过去，exgcd
的返回值时a，b的最大公约数*/
int exgcd(int a,int b,int &x,int &y)
{
int r,t;
if(b==0)
{
x=1;
y=0;
return a;
}
r=exgcd(b,a%b,x,y);
t=x;
x=y;
y=t-a/b*y;
return r;
}
int main()
{
int m,n,x,y,r;
scanf("%d%d",&m,&n);
r=exgcd(m,n,x,y);
printf("(%d)*%d+(%d)*%d=%d\n",x,m,y,n,r);
return 0;
}

核心代码：

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

可运行代码二：非递归C/C++

#include<stdio.h>
#include <iostream>
using namespace std;
/*非递归，a,b代表方程ax+by=gcd(a,b)的a,b
x,y代表方程的一组解*/
int exgcd(int a,int b,int &x,int &y)
{
int x1,y1,x0,y0,r,q;
x0=1; y0=0;
x1=0; y1=1;
x=0; y=1;
r=a%b;
q=(a-r)/b;
while(r)
{
x=x0-q*x1; y=y0-q*y1;
x0=x1; y0=y1;
x1=x; y1=y;
a=b; b=r; r=a%b;
q=(a-r)/b;
}
return b;
}
int main()
{
int m,n,x,y,r;
scanf("%d%d",&m,&n);
r=exgcd(m,n,x,y);
printf("(%d)*%d+(%d)*%d=%d\n",x,m,y,n,r);
return 0;
}

核心代码：

int exgcd(int a,int b,int &x,int &y)
{
int x1,y1,x0,y0,r,q;
x0=1; y0=0;
x1=0; y1=1;
x=0; y=1;
r=a%b;
q=(a-r)/b;
while(r)
{
x=x0-q*x1; y=y0-q*y1;
x0=x1; y0=y1;
x1=x; y1=y;
a=b; b=r; r=a%b;
q=(a-r)/b;
}
return b;
}

对比①②的两种算法程序，我感觉②较好一些。
方程的多组解：

当求出一组解x，y之后，我们可以求出所有的解，其它的整数解满足
X=x+(b/gcd(a,b))*t （t为任意整数） Y=y-(a/gcd(a,b))*t （t为任意整数）
二、扩展欧几里得应用[b]（解不定方程、解线性同余方程、解模的逆元）[/b]

1.扩展欧几里得解不定方程
不定方程：形如px+qy=c的方程称为不定方程。

使用扩展欧几里得算法解不定方程的办法：

对于不定整数方程px+qy=c，若 c mod gcd(p,q)=0,则该方程存在整数解，否则不存在整数解。

上面的欧几里得算法已经详细给出求方程ax+by= gcd(a,b)一个整数解的方法，在找到方程ax+by = gcd(a,b)的一组解x1,y1后，方程ax+by = gcd(a,b)的其他整数解满足：

x = x1+(b/gcd(a,b)）*t (t为任意整数) y=y1-(a/gcd(a,b）)*t (t为任意整数)

那么px+qy=c的整数解，只需将ax+by = gcd(a,b)的每个解乘上 c/ gcd(a,b) 即可。

代码如下：C/C++代码

可运行代码：

#include<stdio.h>
#include <iostream>
using namespace std;
/*递归，a,b代表方程ax+by=gcd(a,b)的a,b;
x,y代表方程的一组整数解*/
int exgcd(int a,int b,int &x,int &y)
{
int r,t;
if(b==0)
{
x=1;
y=0;
return a;
}
r=exgcd(b,a%b,x,y);
t=x;
x=y;
y=t-a/b*y;
return r;
}
//扩展欧几里得解不定方程px+qy=c的一组整数解x,y
bool linear_equation(int a,int b,int c,int &x,int &y)
{
int d=exgcd(a,b,x,y);//a,b的最大公约数
if(c%d)
return false;
else
{
int k=c/d;//k为两方程ax+by=gcd(a,b)与px+qy=c相差的系数
x*=k; y*=k;//求得的只是其中一组解
return true;
}
}
int main()
{
int i,p,q,c,x=0,y=0,r;
scanf("%d%d%d",&p,&q,&c);
if(linear_equation(p,q,c,x,y))
printf("x=%d y=%d\n",x,y);
return 0;
}

核心代码：

int exgcd(int a,int b,int &x,int &y)
{
int r,t;
if(b==0)
{
x=1;
y=0;
return a;
}
r=exgcd(b,a%b,x,y);
t=x;
x=y;
y=t-a/b*y;
return r;
}
bool linear_equation(int a,int b,int c,int &x,int &y)
{
int d=exgcd(a,b,x,y);//a,b的最大公约数
if(c%d)
return false;
else
{
int k=c/d;//k为两方程ax+by=gcd(a,b)与px+qy=c相差的系数
x*=k; y*=k;//求得的只是其中一组解
return true;
}
}

当求出一组解x，y之后，我们可以求出所有的解，其它的整数解满足

X=x+(b/gcd(a,b))*t （t为任意整数） Y=y-(a/gcd(a,b))*t （t为任意整数）

2.扩展欧几里得算法解线性同余方程

线性同余方程：形如ax≡b (mod n)的方程称为线性同余方程。

性质：当且仅当 b 能够被 a 与 n 的最大公约数d整除（记作 gcd(a,n) | b）时，方程有解，并且方程在[0,n-1]之间恰有gcd(a,n)个解（注意不是方程只有gcd(a,n)个解）。

例如：

①在方程6x ≡ 7 (mod 8)中，d=gcd(6,8)=2，7不能整除2，因此方程无解。

②在方程6x ≡ 10 (mod 8)中，d=gcd(6,8)=2，10能整除2，因此方程有解，并且方程在[0,7]上有2个解，

分别为x=3，x=7。

③在方程21x ≡ 7 (mod 3)中，d=gcd(21,3)=3，7不能整除3，因此方程无解。

④在方程10x ≡ 5 (mod 5)中，d=gcd(10,5)=5，5能整除5，因此方程有解，并且方程在[0,4]上有5个解，

分别为x=0，x=1，x=2，x=3，x=4。

⑤在方程30x ≡ 10 (mod 15)中，d=gcd(30,15)=15，10不能整除15，因此方程无解。

⑥在方程15x ≡ 40 (mod 10)中，d=gcd(15,10)=5，40能整除5，因此方程有解，并且方程在[0,14]上有5个解，

分别为x=0，x=2，x=4，x=6，x=8。

求解：求解线性同余方程ax≡b (mod n)相当于求解方程ax+ny= b, (x, y均为整数)。根据线性同余方程的性质可知，当gcd(a,n) | b时，方程有解，在这里我们讨论有解的情况。

设d= gcd(a,n)，则方程ax+ny=d的解可以用扩展欧几里得算法求出一组x0,y0，则有a*x0+n*y0=d，方程两边同时乘以b/d（因为d能整除b）就得到a*x0*b/d+n*y0*b/d=d*b/d=b，所以x=x0*b/d,y=y0*b/d是方程ax+ny=b的一组整数解，所以x=x0*b/d是线性同余方程ax≡b (mod n)的一个解。（设ans=x*(b/d),s=n/d，则方程ax≡b
(mod n)的最小整数解为(ans%s+s)%s）

方程ax≡b (mod n)的一个解为 x1= x* (b/ d ) mod n，且方程的 d 个解分别为 Xi= (x1+ i* (n/ d ))mod n {i= 0... d-1}

求解d个解的代码如下：

输入：a b n

输出：No answer!或d个解

#include<stdio.h>
#include<iostream>
using namespace std;
int exgcd(int a,int b,int &x,int &y)
{
int r,t;
if(b==0)
{
x=1;
y=0;
return a;
}
r=exgcd(b,a%b,x,y);
t=x;
x=y;
y=t-a/b*y;
return r;
}
bool modular_linear_equation(int a,int b,int n)
{
int x,y,x0,i;
int d=exgcd(a,n,x,y);
if(b%d)
{
printf("No answer!\n");
return false;
}
x0=(x*(b/d)%(n/d)+n/d)%(n/d);//最小整数解
for(i=0;i<d;i++)
printf("x=%d\n",(x0+i*(n/d))%n);
return true;
}
int main()
{
int a,b,n,x,y;
scanf("%d%d%d",&a,&b,&n);
modular_linear_equation(a,b,n);
return 0;
}

线性同余方程若有解，则有多个解，那么解与解之间有没有关系呢？有什么关系呢？

对上面性质中的举例，我们可以些许发现，当方程有解时，它们的解的间隔是相同的。那么这个间隔与方程的a,b,n有怎样的关系呢？

证明：我们设解之间的间隔为dx，因为ax≡b(mod n)所以a(x+dx)≡b (mod n)，两式相减得a*dx (mod n)=0，那么a*dx即是a的倍数也是n的倍数，则a*dx是a,n的公倍数，要求出dx的最小值，我们只需求出a,n的最小公倍数，则此时对应的dx最小。

设d是a和n的最大公约数，则a和n的最小公倍数为(a*n)/d。所以a*dx=a*n/d，则dx=n/d。

求线性同余方程的解的代码如下：

#include<stdio.h>
#include<iostream>
using namespace std;
int exgcd(int a,int b,int &x,int &y)
{
int r,t;
if(b==0)
{
x=1;
y=0;
return a;
}
r=exgcd(b,a%b,x,y);
t=x;
x=y;
y=t-a/b*y;
return r;
}
bool modular_linear_equation(int a,int b,int n)
{
int x,y,x0,i;
int d=exgcd(a,n,x,y);
if(b%d)
return false;
x0=x*(b/d)%n;//特解
for(i=0;i<100;i++)//控制输出多少解
printf("%d\n",x0+i*(n/d));
return true;
}int main()
{
int a,b,n,x,y;
scanf("%d%d%d",&a,&b,&n);
modular_linear_equation(a,b,n);
return 0;
}

核心代码：

bool modular_linear_equation(int a,int b,int n)
{
int x,y,x0,i;
int d=exgcd(a,n,x,y);
if(b%d)
return false;
x0=x*(b/d)%n;//特解
for(i=0;i<100;i++)//控制输出多少解
printf("%d\n",x0+i*(n/d));
return true;
}

3.求模的乘法逆元

乘法逆元：若ax≡1 mod f, 则称x为a关于模f的乘法逆元。也可表示为ax≡1(mod f)。

当a与f互素时，a关于模f的乘法逆元有唯一解（因为方程ax≡1 mod f在 内只有一个解）。如果不互素，则无解。如果f为素数，则从1到f-1的任意数都与f互素，即在1到f-1之间都恰好有一个关于模f的乘法逆元。

代码如下：

#include <stdio.h>
#include<iostream>
using namespace std;
/* 扩展欧几里得算法求乘法逆元*/
int ExtendedEuclid(int f,int d,int *result)
{
int x1,x2,x3,y1,y2,y3,t1,t2,t3,q;
x1=y2=1;
x2=y1=0;
x3=(f>=d)?f:d;
y3=(f>=d)?d:f;
while(1)
{
if(y3==0)
{
*result=x3;//两个数不互素则result为两个数的最大公约数，此时返回值为零
return 0;
}
if(y3==1)
{
*result = y2;//两个数互素则resutl为其乘法逆元，此时返回值为1
return 1;
}
q=x3/y3;
t1=x1-q*y1;
t2=x2-q*y2;
t3=x3-q*y3;
x1=y1;
x2=y2;
x3=y3;
y1=t1;
y2=t2;
y3=t3;
}
}
int main()
{
int x,y,z;
z = 0;
printf("输入两个数：\n");
scanf("%d%d",&x,&y);
if(ExtendedEuclid(x,y,&z))
printf("%d和%d互素，乘法的逆元是：%d\n",x,y,z);
else
printf("%d和%d不互素，最大公约数为：%d\n",x,y,z);
return 0;
}

核心代码：

int ExtendedEuclid(int f,int d,int *result)
{
int x1,x2,x3,y1,y2,y3,t1,t2,t3,q;
x1=y2=1;
x2=y1=0;
x3=(f>=d)?f:d;
y3=(f>=d)?d:f;
while(1)
{
if(y3==0)
{
*result=x3;
return 0;
}
if(y3==1)
{
*result = y2;
return 1;
}
q=x3/y3;
t1=x1-q*y1;
t2=x2-q*y2;
t3=x3-q*y3;
x1=y1;
x2=y2;
x3=y3;
y1=t1;
y2=t2;
y3=t3;
}
}