【数论】[总结]扩展欧几里德算法
2012-12-08 21:47
232 查看
NOIP2012 Day2 T1 的 mod 一题会用到扩展欧几里德,本人理解能力有限,花了两三天终于理解了,这里总结一下
(以下的所有符号以C/C++中的符号为准,特别是 / 在下面的情况下相当于Pascal的div)
欧几里德的辗转相除法
现在再看这样一个等式ax+by=c,如果 c % gcd(a,b) =0 的话,那么方程就必定有整数解
证明很容易,ax % gcd(a,b) =0,而且 bx % gcd(a,b) =0,那么 (ax+by) % gcd(a,b) =0,又 c % gcd(a,b) =0 ,所以ax+by=c必定有整数解
那么 ax+by=1 也必定有整数解(因为 c%1=0 是显然成立的 ),所以如果我们能求出 ax+by=1 的一组整数解x0,y0,那么 ax+by=c 的整数解就可以求得 x=c*x0,y=c*y0
现在我们求得了一组整数解x,y,但是这个x和y是其中的一组,比如我们要找最小的x该怎么办呢?我们可以由原方程ax+by得到 ax+by=a(x+b)+b(y-a) ,那么我们求的x最小整数解就可以通过+b来求得了,但是时间一下子又慢了(如果求得的x是一个很小的负数),所以我们可以根据负数取余的特性(绝对值取余后再加上符号,例如 -13%3 = -(13%3) )快速得到最小正整数解:(x%b+b)%b <由于只求x,所以就没有处理y>
现在的问题就转化成了如何求 ax+by=1 的一组整数解了
欧几里得的算法是 gcd(a,b)=gcd(b,a%b) ,即 ax+by=gcd(a,b)成立,gcd(a,b)=gcd(b,a%b) ,所以bx+(a%b)y=gcd(b,a%b)=gcd(a,b),那么同理 bx+(a%b)y=gcd(b,a%b) 成
立。为了表示方便,我们直接把gcd(a,b)写成1来表示(这里gcd(a,b)的值就是1)
由数学上的知识很容易得到 a%b = a-a/b*b
上面方程bx+(a%b)y=1就可以写成 bx+(a-a/b*b)y=1 化简得到 ==> bx+ay-(a/b*b)y=1
好了关键就在这里,后面的a和b和前面ax+by=1的a和b是相同的,但是x和y不同,那么x和y之间有什么关系呢?我们把它们整理一下
ax+by=1
ay+b(x-a/b*y)=1
也就是说在进行了一次gcd后,以前的x变成了y,而y变成了x-a/b*y 。。。。。。。。。这个结论先放着,等会再用
gcd进行到边界,也就是b=0的时候返回a,这个时候我们看看方程ax+by=1,就成了ax=1,我们再想想这个时候的a是多少?对了,a肯定等于1!因为这时候的a=gcd(a,b)=1 ,所以很自然的我们解得 x=1,y=0
现在就把刚才的结论用到这里开始进行倒推,最后求得一组x,y的整数解
这个时候就可以把扩展欧几里德的代码放上来了
当然,如果你的目的仅仅是为了求得一组x,y的话,那么就可以定义为void,而不用返回值了,当然定义成int也没错。。。
(以下的所有符号以C/C++中的符号为准,特别是 / 在下面的情况下相当于Pascal的div)
欧几里德的辗转相除法
int gcd(int a,int b) { return (b)?gcd(b,a%b):a; }
现在再看这样一个等式ax+by=c,如果 c % gcd(a,b) =0 的话,那么方程就必定有整数解
证明很容易,ax % gcd(a,b) =0,而且 bx % gcd(a,b) =0,那么 (ax+by) % gcd(a,b) =0,又 c % gcd(a,b) =0 ,所以ax+by=c必定有整数解
那么 ax+by=1 也必定有整数解(因为 c%1=0 是显然成立的 ),所以如果我们能求出 ax+by=1 的一组整数解x0,y0,那么 ax+by=c 的整数解就可以求得 x=c*x0,y=c*y0
现在我们求得了一组整数解x,y,但是这个x和y是其中的一组,比如我们要找最小的x该怎么办呢?我们可以由原方程ax+by得到 ax+by=a(x+b)+b(y-a) ,那么我们求的x最小整数解就可以通过+b来求得了,但是时间一下子又慢了(如果求得的x是一个很小的负数),所以我们可以根据负数取余的特性(绝对值取余后再加上符号,例如 -13%3 = -(13%3) )快速得到最小正整数解:(x%b+b)%b <由于只求x,所以就没有处理y>
现在的问题就转化成了如何求 ax+by=1 的一组整数解了
欧几里得的算法是 gcd(a,b)=gcd(b,a%b) ,即 ax+by=gcd(a,b)成立,gcd(a,b)=gcd(b,a%b) ,所以bx+(a%b)y=gcd(b,a%b)=gcd(a,b),那么同理 bx+(a%b)y=gcd(b,a%b) 成
立。为了表示方便,我们直接把gcd(a,b)写成1来表示(这里gcd(a,b)的值就是1)
由数学上的知识很容易得到 a%b = a-a/b*b
上面方程bx+(a%b)y=1就可以写成 bx+(a-a/b*b)y=1 化简得到 ==> bx+ay-(a/b*b)y=1
好了关键就在这里,后面的a和b和前面ax+by=1的a和b是相同的,但是x和y不同,那么x和y之间有什么关系呢?我们把它们整理一下
ax+by=1
ay+b(x-a/b*y)=1
也就是说在进行了一次gcd后,以前的x变成了y,而y变成了x-a/b*y 。。。。。。。。。这个结论先放着,等会再用
gcd进行到边界,也就是b=0的时候返回a,这个时候我们看看方程ax+by=1,就成了ax=1,我们再想想这个时候的a是多少?对了,a肯定等于1!因为这时候的a=gcd(a,b)=1 ,所以很自然的我们解得 x=1,y=0
现在就把刚才的结论用到这里开始进行倒推,最后求得一组x,y的整数解
这个时候就可以把扩展欧几里德的代码放上来了
int gcd(int a,int b,int &x,int &y)//一个函数,返回gcd(a,b) { if(b==0) { x=1;y=0;//边界时候的赋值 return a; } else{ int ans=gcd(b,a%b); int t=x;x=y;//这一句和后一句就是刚才推到的求法 y=t-a/b*y; return ans; } }
当然,如果你的目的仅仅是为了求得一组x,y的话,那么就可以定义为void,而不用返回值了,当然定义成int也没错。。。