扩展欧几里得算法及其应用
一、扩展欧几里得算法
扩展欧几里得算法:对于不完全为 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; }
#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为任意整数)
二、扩展欧几里得应用(解不定方程、解线性同余方程、解模的逆元)
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; }
乘法逆元:若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; } }
扩展欧几里得算法------扩展欧几里德算法,布布扣,bubuko.com
原文:http://blog.csdn.net/yanghuaqings/article/details/38440789