扩展欧几里得就是已知a, b求解一组x,y,使它们满足等式: ax+by = gcd(a, b) =c,它的解是一定存在的,只要c是gcd的倍数。
在说扩展欧几里得算法之前需要知道朴素欧几里得算法,即gcd(a,b) = gcd(b,a%b),(最大公约数Greatest Common Divisor),证明也很简单,我采用的是反证法,这里不再赘述。
突然发现一篇博客写的非常非常好,这里直接转载了。(侵删)
当我真正理解了扩展欧几里得定理
首先、扩展欧几里得定理:对于两个不全为0的整数a、b,必存在一组解x,y,使得ax+by==gcd(a,b);
实现如下:
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))。很好理解吧~
再深入一点,就解出这么一组解其实一般来说是解决不了什么问题的。没有哪个ACM的题这么简单吧。。。比如我们现在要得到所有的解,那么这所有的解究竟是什么呢?
直接说吧,假设d=gcd(a,b). 那么x=x0+b/d*t; y=y0-a/d*t;其中t为任意常整数。
这个是怎么推导出来的,说实话我也不知道,就先这么记着吧!
好了,说了这么多,光说不练也无济于事
先做一个最简单的题目,pku 1061青蛙的约会
先说一下大概题意:有两只青蛙,一只在坐标x,另一直在坐标y,青蛙x一次跳跃可以前进m单位距离,青蛙y一次跳跃可以前进n单位的距离,两青蛙都在同一纬度,该纬度长度为L。两只青蛙同方向同时跳啊跳,问你最少跳多少次,它们才可以相遇,如果不能相遇,输出impossble
分析:假设跳了T次以后,青蛙1的坐标便是x+m*T,青蛙2的坐标为y+n*T。它们能够相遇的情况为(x+m*T)-(y+n*T)==P*L,其中P为某一个整数,变形一下
得到(n-m)*T+P*L==x-y 我们设a=(n-m),b=L,c=x-y,T=x,P=y.于是便得到ax+by==c。激动啊,这不就是上面一样的式子吗~
直接套用扩展欧几里得函数,得到一组解x,y。由于问题是问最少跳多少次,于是只有x是我们需要的信息。那么再想,x是最小的吗?
答案是可能不是!那么如何得到最小解呢? 我们考虑x的所有解的式子: x=x0+b/d*t。x0是我们刚刚求到的,很显然右边是有个单调函数,当t为某一个与x正负性质相反的数时,可以得到最小的x。 令x的正负性质为正,那么x=x0-b/d*t1 (t1==-t)。令x==0,那么t=x0*d/b,最小的x等于x0减去t*b/d。这里得到的x可能是负数,如果是负数,我们再为它加上一个b/d即是所求答案了!
#include<iostream>
#include<string>
#include<cmath>
#include<algorithm>
usingnamespace std;
__int64 x,y,a,b,c,d;
__int64 n,m,X,Y,L;
__int64 gcd(__int64 a,__int64 b)
{
__int64 t,d;
if(b==0)
{
x=1;
y=0;
return a;
}
d=gcd(b,a%b);
t=x;
x=y;
y=t-(a/b)*y;
return d;
}
int main()
{
while(scanf("%I64d%I64d%I64d%I64d%I64d",&X,&Y,&m,&n,&L)==5)
{
a=n-m;
b=L;
c=X-Y;
d=gcd(a,b);
if(c%d!=0)
{
printf("Impossible\n");
continue;
}
x=x*(c/d);
y=y*(c/d);
/*通解:
x1=x+b/d*t;
y1=y-a/d*t;
t为任意整数
*/
//找最小的x1,即求x+b/d*t最小,那么只有t为某一个数时才最小
//显然t必须与x正负相反才有最小,那么就看做x-b/d*t,这个式子的最小值便是t=x/(b/d)时,注意这是整型除法
__int64 k=x*d/b;
k=x-k*b/d;
if(k<0)
k+=b/d;
printf("%I64d\n",k);
}
return0;
}
再来看一题难一点的 pku 2142 the balance http://poj.org/problem?id=2142
还是先说一下题意吧:一个家伙有一种天平,这种天平只有两种重量的砝码a和b,现在要称出重量为c的物品,问你至少需要多少a和b,答案需要满足a的数量加上b的数量和最小,并且他们的重量和也要最小。(两个盘都可以放砝码)
分析:
这题我刚刚开始还以为是动规或者搜索,也算是碰了一鼻子的灰吧。
假设a砝码我们用了x个,b砝码我们用了y个。那么天平平衡时,就应该满足ax+by==c。x,y为正时表示放在和c物品的另一边,为负时表示放在c物品的同一边。
于是题意就变成了求|x|+|y|的最小值了。x和y是不定式ax+by==c的解。
刚刚上面已经提到了关于x,y的所以解的同式,即
x=x0+b/d*t
y=y0-a/d*t
你是不是下意识的想要穷举所有解,取|x|+|y|最小的?显然是行不通的,仔细分析:|x|+|y|==|x0+b/d*t|+|y0-a/d*t|,我们规定a>b(如果a<b,我们便交换a b),从这个式子中,我们可以轻易的发现:|x0+b/d*t|是单调递增的,|y0-a/d*t|是单调递减的,而由于我们规定了a>b,那么减的斜率边要大于增的斜率,于是整个函数减少的要比增加的快,但是由于绝对值的符号的作用,最终函数还是递增的。也就是说,函数是凹的,先减小,再增大。那么什么时候最小呢?很显然是y0-a/d*t==0的时候,于是我们的最小值|x|+|y|也一定是在t=y0*d/a附近了,我是在t点左右5个点的范围内取最小的(据说左右一个点都可以,不过我试了一下wa了)。
#include<iostream>
#include<string>
#include<queue>
#include<algorithm>
usingnamespace std;
#define min(a,b) a<b?a:b
#define max(a,b) a>b?a:b
int x,y;
int a,b,c;
int gcd(int a,int b)
{
int d;
if(b==0)
{
x=1;y=0;
return a;
}
else
{
d=gcd(b,a%b);
int t=x;
x=y;
y=t-(a/b)*y;
}
return d;
}
int main()
{
int d;
//freopen("e:\\in.txt","r",stdin);
while(scanf("%d%d%d",&a,&b,&c)!=EOF)
{
if(!a &&!b &&!c)
break;
int flag=0;
if(a<b)
{
flag=1;
swap(a,b);
}
d=gcd(a,b);
x=x*c/d; y=y*c/d;
int t=y*d/a,ans=0xfffffff;
int x1,y1,x2,y2;
for(int i=t-5;i<=t+5;i++)
{
x2=x+(b/d)*i;
y2=y-(a/d)*i;
if(abs(x2)+abs(y2)<ans)
{
ans=abs(x2)+abs(y2);
x1=x2;
y1=y2;
}
}
if(!flag)
printf("%d %d\n",abs(x1),abs(y1));
else
printf("%d %d\n",abs(y1),abs(x1));
}
return0;
}
嘿嘿,每次搞懂新问题都很开心~
再来个更难一点的 PKU 2891 Strange Way to Express Integers http://poj.org/problem?id=2891
大概题意是,给你很多同余关系,要你求一个数,使得这个数c mod a1==b1 ,c mod a2==b2,,,,,, c mod ak==bk(a1,a2,a3....ak之间可能不互质)
这个其实就是中国剩余定理的欧几里得解法,如果(a1,a2,,,,,ak)两两互质的话,那么就可以直接套用中国剩余定理的公式了,可是这题可能不互质。那么就只能用扩展欧几里得来解了,对于头两个关系式,c % a1==b1, c % a2==b2,可以联立方程组 a1*x+a2*y==b2-b1 ,很熟悉,直接套用扩展欧几里得算法得到一组解x,y,再取x的最小值,带入式子c%a1==b1中,便得到了答案 c 。
这是只有两个式子的,很简单,那么对于剩下的第三个第四个怎么办呢?刚刚得到的x已经满足了头两个同余式,对于新来的第三个式子,他们得到的数必须要前面两个式子的要求,同时还要满足第三个新式子的要求,那么我们便将前两个式子合并成一个式子,用这个新式子与第三个式子联立求解,然后依次类推,知道解出了所有式子的解。
那么如何合并呢?
新式子公式是这样的:c1 mod lcm(a1,a2) ==c 这里c1是要求的新答案,c是前两个式子求的的一个解,lcm(a1,a2)是a1,a2的最小公倍数
你可以理解为新的答案c1即要mod a1,又要mod a2,所以就mod lcm(a1,ba)。而余数c 是刚前面求得的答案,也就是说c1不管怎么折腾,都还会保留一个c,即不会破坏前面的成果!
#include<iostream>
#include<string>
usingnamespace std;
__int64 b,b1,c,c1,d,x,y,A,B,C;
__int64 gcd(__int64 a,__int64 b)
{
__int64 d,t;
if(!b)
{
x=1;
y=0;
return a;
}
d=gcd(b,a%b);
t=x;
x=y;
y=t-(a/b)*y;
return d;
}
int main()
{
int k,i,flag;
freopen("in.txt","r",stdin);
while(scanf("%d",&k)!=EOF)
{
flag=0;
scanf("%I64d%I64d",&b,&c);
if(k==1)
{
printf("%I64d\n",b+c);
continue;
}
for(i=1;i<k;i++)
{
scanf("%I64d%I64d",&b1,&c1);
if(flag)
continue;
A=b; B=b1; C=c1-c;
d=gcd(A,B);
if(C%d!=0)
{
flag=1;
continue;
}
x=x*C/d; //取得x的最小值
x=x-(x*d/B)*(B/d);
if(x<0)
x+=B/d;
c=x*b+c; b=(b*b1)/d; //按前面的分析头两合并式子
}
if(!flag)
{
printf("%I64d\n",c);
}
else
{
printf("-1\n");
}
}
return0;
}