先介绍什么叫做欧几里德算法(辗转相除法)
有两个数 a b,现在,我们要求 a b 的最大公约数,怎么求?枚举他们的因子?不现实,当 a b 很大的时候,枚举显得那么的na?ve ,那怎么做?
欧几里德有个十分又用的定理: gcd(a, b) = gcd(b , a%b) ,这样,我们就可以在几乎是 log 的时间复杂度里求解出来 a 和 b 的最大公约数了,这就是欧几里德算法,用 C++ 语言描述如下:
#include
#include
#include
using namespace std;
int gcd(int a,int b)
{
if(b==0)return a;
else gcd(b,(a%b));
}
int main()
{
int n,m;
scanf("%d%d",&n,&m);
printf("gcd=%d\n",gcd(n,m));
printf("lcm=%d\n",n*m/gcd(n,m));
return 0;
}
由于是用递归写的,所以看起来很简洁,也很好记忆。那么什么是扩展欧几里德呢?
现在我们知道了 a 和 b 的最大公约数是 gcd ,那么,我们一定能够找到这样的 x 和 y ,使得: a*x + b*y = gcd 这是一个不定方程(其实是一种丢番图方程),有多解是一定的,但是只要我们找到一组特殊的解 x0 和 y0 那么,我们就可以用 x0 和 y0 表示出整个不定方程的通解:
x = x0 + (b/gcd)*t
y = y0 – (a/gcd)*t
通解推导过程:ax + by = gcd(a,b) ①
ax0 + by0 = gcd(a,b) ②
①﹣②:a(x-x0)+b(y-y0)=0
a(x-x0)=b(y0-y)
将两边同除以gcd(a,b),得 a/gcd(a,b)与b/gcd(a,b)互质,
所以 a/gcd(a,b)(x-x0)=b/gcd(a,b)(y0-y)
所以 x-x0=b/gcd(a,b) *t;y0-y=a/gcd(a,b)*t
所以 x=x0+b/gcd(a,b)*t
y=y0-a/gcd(a,b)*t
扩展欧几里德算法的全部过程,依然用递归写
int gcd(int a,int b)
{
int t,d;
if(b==0)
{
x=1;
y=0;
//不明处1
return a;
}
d=gcd(b,a%b);
t=x;
x=y;
y=t-(a/b)*y;
//不明处2
return d;
}
上面的程序中,x和y我是用全局变量保存的
我个人觉得第一次看到这个程序你会有以上两个不明白的地方(见注释),下面我分别解释
不明处1:由扩展欧几里得定理:ax+by==gcd(a,b)---式1,而此时b==0,也就是说gcd(a,0)==a。原式变为ax+by==a --> x==1,y==0。应该够清楚了吧
【扩展欧几里德算法详解(通解推导过程)】 不明处2:这里先说明一下我的一些规则,x,y表示第一次递归时的值,x1,y1表示第二次递归时的值。那么
gcd(a,b)==gcd(b,a%b),同时都代入式1,有ax+by==b*x1+(a%b)*y1。将右边变形一下
b*x1+(a%b)*y1==b*x1+(a-(a/b)*b)*y1==a*y1+b*(x1-(a/b)*y1),最终得到ax+by==a*y1+b*(x1-(a/b)*y1)
也就是说,上一深度的x等于下一深度的y1,上一深度的y等于下一深度的x1-(a/b)*y1。需要注意,上面推导时用的除法都是整型除法
到这里为止,我们便得到了不定式ax+by==gcd(a,b)的一组解,x、y。
那么对于一般的不定式ax+by==c,它的解应该是什么呢。很简单,x1=x*(c/gcd(a,b)),y1=y*(c/gcd(a,b))。
依然很简短,相比欧几里德算法,只是多加了几个语句而已。
这就是理论部分,欧几里德算法部分我们好像只能用来求解最大公约数,但是扩展欧几里德算法就不同了,我们既可以求出最大公约数,还可以顺带求解出使得: a*x + b*y = gcd 的通解 x 和 y
扩展欧几里德有什么用处呢?
扩展欧几里德算法的应用主要有以下三方面:
(1)求解不定方程;
(2)求解模线性方程(线性同余方程);
(3)求解模的逆元;
(1)使用扩展欧几里德算法解决不定方程的办法:
对于不定整数方程pa+qb=c,若 c mod Gcd(p, q)=0,则该方程存在整数解,否则不存在整数解。
上面已经列出找一个整数解的方法,在找到p * a+q * b = Gcd(p, q)的一组解p0,q0后,p * a+q * b = Gcd(p, q)的其他整数解满足:
p = p0 + b/Gcd(p, q) * t
q = q0 - a/Gcd(p, q) * t(其中t为任意整数)
至于pa+qb=c的整数解,只需将p * a+q * b = Gcd(p, q)的每个解乘上 c/Gcd(p, q) 即可。
在找到p * a+q * b = Gcd(a, b)的一组解p0,q0后,应该是得到p * a+q * b = c的一组解p1 = p0*(c/Gcd(a,b)),q1 = q0*(c/Gcd(a,b)),
p * a+q * b = c的其他整数解满足:
p = p1 + b/Gcd(a, b) * t
q = q1 - a/Gcd(a, b) * t(其中t为任意整数)
p 、q就是p * a+q * b = c的所有整数解。
相关证明可参考:http://www.cnblogs.com/void/archive/2011/04/18/2020357.html
用扩展欧几里得算法解不定方程ax+by=c;
代码如下:
bool linear_equation(int a,int b,int c,int &x,int &y)
{
int d=exgcd(a,b,x,y);
if(c%d)
return false;
int k=c/d;
x*=k;
y*=k;
//求得的只是其中一组解
return true;
}
(2)用扩展欧几里德算法求解模线性方程的方法:
同余方程 ax≡b (mod n)对于未知数 x 有解,当且仅当 gcd(a,n) | b。且方程有解时,方程有 gcd(a,n) 个解。
求解方程 ax≡b (mod n) 相当于求解方程 ax+ ny= b, (x, y为整数)
设 d= gcd(a,n),假如整数 x 和 y,满足 d= ax+ ny(用扩展欧几里德得出)。如果 d| b,则方程
a* x0+ n* y0= d, 方程两边乘以 b/ d,(因为 d|b,所以能够整除),得到 a* x0* b/ d+ n* y0* b/ d= b。
所以 x= x0* b/ d,y= y0* b/ d 为 ax+ ny= b 的一个解,所以 x= x0* b/ d 为 ax= b (mod n ) 的解。
ax≡b (mod n)的一个解为 x0= x* (b/ d ) mod n,且方程的 d 个解分别为 xi= (x0+ i* (n/ d ))mod n {i= 0... d-1}。
设ans=x*(b/d),s=n/d;
方程ax≡b (mod n)的最小整数解为:(ans%s+s)%s;
相关证明:
证明方程有一解是: x0 = x'(b/d) mod n;
由 a*x0 = a*x'(b/d) (mod n)
a*x0 = d (b/d) (mod n)(由于 ax' = d (mod n))
= b (mod n)
证明方程有d个解: xi = x0 + i*(n/d)(mod n);
由 a*xi (mod n) = a * (x0 + i*(n/d)) (mod n)
= (a*x0+a*i*(n/d)) (mod n)
= a * x0 (mod n)(由于 d | a)
= b
首先看一个简单的例子:
5x=4(mod3)
解得x = 2,5,8,11,14.......
由此可以发现一个规律,就是解的间隔是3.
那么这个解的间隔是怎么决定的呢?
如果可以设法找到第一个解,并且求出解之间的间隔,那么就可以求出模的线性方程的解集了.
我们设解之间的间隔为dx.
那么有
a*x = 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是最小的.
设a 和 n的最大公约数为d,那么a 和 n 的最小公倍数为(a*n)/d.
即a*dx = a*n/d;
所以dx = n/d.
因此解之间的间隔就求出来了.
代码如下:
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=1;
i
(3)用欧几里德算法求模的逆元:
同余方程ax≡b (mod n),如果 gcd(a,n)== 1,则方程只有唯一解。
在这种情况下,如果 b== 1,同余方程就是 ax=1 (mod n ),gcd(a,n)= 1。
这时称求出的 x 为 a 的对模 n 乘法的逆元。
对于同余方程 ax= 1(mod n ), gcd(a,n)= 1 的求解就是求解方程
ax+ ny= 1,x, y 为整数。这个可用扩展欧几里德算法求出,原同余方程的唯一解就是用扩展欧几里德算法得出的 x 。
求解形如 a*x +b*y = c 的通解,但是一般没有谁会无聊到让你写出一串通解出来,都是让你在通解中选出一些特殊的解,比如一个数对于另一个数的乘法逆元
什么叫乘法逆元?
文章图片
这里,我们称 x 是 a 关于 m 的乘法逆元
这怎么求?可以等价于这样的表达式: a*x + m*y = 1
看出什么来了吗?没错,当gcd(a , m) != 1 的时候是没有解的这也是 a*x + b*y = c 有解的充要条件: c % gcd(a , b) == 0
接着乘法逆元讲,一般,我们能够找到无数组解满足条件,但是一般是让你求解出最小的那组解,怎么做?我们求解出来了一个特殊的解 x0 那么,我们用 x0 % m其实就得到了最小的解了。为什么?
可以这样思考:
x 的通解不是 x0 + m*t 吗?
那么,也就是说, a 关于 m 的逆元是一个关于 m 同余的,那么根据最小整数原理,一定存在一个最小的正整数,它是 a 关于m 的逆元,而最小的肯定是在(0 , m)之间的,而且只有一个,这就好解释了。
可能有人注意到了,这里,我写通解的时候并不是 x0 + (m/gcd)*t ,但是想想一下就明白了,gcd = 1,所以写了跟没写是一样的,但是,由于问题的特殊性,有时候我们得到的特解 x0 是一个负数,还有的时候我们的 m 也是一个负数这怎么办?
当 m 是负数的时候,我们取 m 的绝对值就行了,当 x0 是负数的时候,他模上 m 的结果仍然是负数(在计算机计算的结果上是这样的,虽然定义的时候不是这样的),这时候,我们仍然让 x0 对abs(m) 取模,然后结果再加上abs(m) 就行了,于是,我们不难写出下面的代码求解一个数 a 对于另一个数 m 的乘法逆元:
文章图片
还有最小整数解之类的问题,但都是大同小异,只要细心的推一推就出来了,这里就不一一介绍了,下面给一些题目还有AC代码,仅供参考
ZOJ 3609 :http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=4712 求最小逆元
#include
#include
#include
#include
#include
#include
#include
#include
#include #define INF 0x7fffffff
#define EPS 1e-12
#define MOD 1000000007
#define PI 3.141592653579798
#define N 100000using namespace std;
typedef long long LL;
typedef double DB;
LL e_gcd(LL a,LL b,LL &x,LL &y)
{
if(b==0)
{
x=1;
y=0;
return a;
}
LL ans=e_gcd(b,a%b,x,y);
LL temp=x;
x=y;
y=temp-a/b*y;
return ans;
}LL cal(LL a,LL b,LL c)
{
LL x,y;
LL gcd=e_gcd(a,b,x,y);
if(c%gcd!=0) return -1;
x*=c/gcd;
b/=gcd;
if(b<0) b=-b;
LL ans=x%b;
if(ans<=0) ans+=b;
return ans;
}int main()
{
LL a,b,t;
scanf("%lld",&t);
while(t--)
{
scanf("%lld%lld",&a,&b);
LL ans=cal(a,b,1);
if(ans==-1) printf("Not Exist\n");
else printf("%lld\n",ans);
}
return 0;
}
ZOJ 3593 http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=3593 求最小的步数,处理特殊一点就过去了
#include
#include
#include
#include
#include
#include
#include
#include
#include #define INF 0x7fffffff
#define EPS 1e-12
#define MOD 100000007
#define PI 3.14159265357979823846
#define N 100005using namespace std;
typedef long long LL;
LL e_gcd(LL a,LL b,LL &x,LL &y)
{
if(b==0)
{
x=1;
y=0;
return a;
}
LL ans=e_gcd(b,a%b,x,y);
LL temp=x;
x=y;
y=temp-a/b*y;
return ans;
}
LL cal(LL a,LL b,LL L)
{
LL x,y;
LL gcd=e_gcd(a,b,x,y);
if(L%gcd!=0) return -1;
x*=L/gcd;
y*=L/gcd;
a/=gcd;
b/=gcd;
LL ans=((LL)INF)*((LL)INF), f;
LL mid=(y-x)/(a+b);
for(LL T=mid-1;
T<=mid+1;
T++)
{
if(abs(x+b*T)+abs(y-a*T)==abs(x+b*T+y-a*T))
f=max(abs(x+b*T),abs(y-a*T));
else
f=fabs(x-y+(a+b)*T);
ans=min(ans,f);
}
return ans;
}int main()
{
//freopen("in.in","r",stdin);
//freopen("out.out","w",stdout);
LL A,B,a,b,x,y;
int t;
scanf("%d",&t);
while(t--)
{
scanf("%lld%lld%lld%lld",&A,&B,&a,&b);
LL L=B-A;
LL ans=cal(a,b,L);
if(ans==-1) printf("-1\n");
else printf("%lld\n",ans);
}
return 0;
}
POJ 1061http://poj.org/problem?id=1061 青蛙的约会,裸的扩展欧几里得
#include
#include
#include
#include
#include
#include
#include
#include
#include #define INF 0x7fffffff
#define EPS 1e-12
#define MOD 1000000007
#define PI 3.141592653579798
#define N 100000using namespace std;
typedef long long LL;
typedef double DB;
LL e_gcd(LL a,LL b,LL &x,LL &y)
{
if(b==0)
{
x=1;
y=0;
return a;
}
LL ans=e_gcd(b,a%b,x,y);
LL temp=x;
x=y;
y=temp-a/b*y;
return ans;
}LL cal(LL a,LL b,LL c)
{
LL x,y;
LL gcd=e_gcd(a,b,x,y);
if(c%gcd!=0) return -1;
x*=c/gcd;
b/=gcd;
if(b<0) b=-b;
LL ans=x%b;
if(ans<=0) ans+=b;
return ans;
}int main()
{
LL x,y,m,n,L;
while(scanf("%lld%lld%lld%lld%lld",&x,&y,&m,&n,&L)!=EOF)
{
LL ans=cal(m-n,L,y-x);
if(ans==-1) printf("Impossible\n");
else printf("%lld\n",ans);
}
return 0;
}
HDU 1576http://acm.hdu.edu.cn/showproblem.php?pid=1576 做点处理即可
#include
#include
#include
#include
#include
#include
#include
#include
#include #define INF 0x7fffffff
#define EPS 1e-12
#define MOD 1000000007
#define PI 3.141592653579798
#define N 100000using namespace std;
typedef long long LL;
typedef double DB;
LL e_gcd(LL a,LL b,LL &x,LL &y)
{
if(b==0)
{
x=1;
y=0;
return a;
}
LL ans=e_gcd(b,a%b,x,y);
LL temp=x;
x=y;
y=temp-a/b*y;
return ans;
}LL cal(LL a,LL b,LL c)
{
LL x,y;
LL gcd=e_gcd(a,b,x,y);
if(c%gcd!=0) return -1;
x*=c/gcd;
b/=gcd;
if(b<0) b=-b;
LL ans=x%b;
if(ans<=0) ans+=b;
return ans;
}int main()
{
LL n,b,t;
scanf("%I64d",&t);
while(t--)
{
scanf("%I64d%I64d",&n,&b);
LL ans=cal(b,9973,n);
if(ans==-1) printf("Impossible\n");
else printf("%lld\n",ans);
}
return 0;
}
HDU 2669http://acm.hdu.edu.cn/showproblem.php?pid=2669 裸的扩展欧几里得
#include
#include
#include
#include
#include
#include
#include
#include
#include #define INF 0x7fffffff
#define EPS 1e-12
#define MOD 1000000007
#define PI 3.141592653579798
#define N 100000using namespace std;
typedef long long LL;
typedef double DB;
LL e_gcd(LL a,LL b,LL &x,LL &y)
{
if(b==0)
{
x=1;
y=0;
return a;
}
LL ans=e_gcd(b,a%b,x,y);
LL temp=x;
x=y;
y=temp-a/b*y;
return ans;
}LL cal(LL a,LL b,LL c)
{
LL x,y;
LL gcd=e_gcd(a,b,x,y);
if(c%gcd!=0) return -1;
x*=c/gcd;
b/=gcd;
if(b<0) b=-b;
LL ans=x%b;
if(ans<=0) ans+=b;
return ans;
}int main()
{
LL a,b;
while(scanf("%I64d%I64d",&a,&b)!=EOF)
{
LL ans=cal(a,b,1);
if(ans==-1) printf("sorry\n");
else printf("%I64d %I64d\n",ans,(1-ans*a)/b);
}
return 0;
}
推荐阅读
- HDU 5528【2015长春现场赛 B】 Count a * b
- 类欧几里得算法|[类欧几里得算法 数论] BZOJ 2987 Earthquake
- [数论] Codeforces 819D R #421 D.Mister B and Astronomers & 516E R #292 E. Drazil and His Happy Friends
- 模板 poj2947 Widget Factory 高斯消元
- 【扩展欧几里得】练习题
- 扩展欧几里得【数论
- 数论|hdu 5322 Hope(分治+NTT)
- HDU 5528 Count a × b
- 线性同余方程组
- 数论|AtCoder Beginner Contest 156 C.Rally