欧几里德算法及扩展欧几里德 中科大 ACM 1209

欧几里德算法
欧几里德算法又称辗转相除法，用于计算两个整数a,b的最大公约数。其计算原理依赖于下面的定理：

定理：gcd(a,b) = gcd(b,a mod b)

证明：a可以表示成a = kb + r，则r = a mod b
假设d是a,b的一个公约数，则有
d|a, d|b，而r = a - kb，因此d|r
因此d是(b,a mod b)的公约数

假设d 是(b,a mod b)的公约数，则
d | b , d |r ，但是a = kb +r
因此d也是(a,b)的公约数

因此(a,b)和(b,a mod b)的公约数是一样的，其最大公约数也必然相等，得证

 

可以递归或者迭代写出程序：

 

int gcd(int n,int m)
{
 int max,min;
 int r;
 n>m ? (max=n,min=m):(max=m,min=n);
 while(min)
 {
  r=max%min;
  max=min;
  min=r;
 }
 return max;
}

 

int gcd_(int n, int m)
{
 int max,min;
 n>m ? (max=n,min=m):(max=m,min=n);
 if(min==0)
  return max;
 return gcd_(min,max%min);
}

 

 

扩展欧几里德算法不但能计算(a，b）的最大公约数，而且能计算a模b及b模a的乘法逆元，扩展欧几里德算法是用来在已知a, b求解一组p，q使得p * a+q  * b = Gcd(a, b)  (解一定存在，根据数论中的相关定理)。扩展欧几里德常用在求解模线性方程及方程组中。

 

 

把这个实现和Gcd的递归实现相比，发现多了下面的x,y赋值过程，这就是扩展欧几里德算法的精髓。
可以这样思考:
对于a' = b, b' = a % b 而言，我们求得 x, y使得 a'x + b'y = Gcd(a', b')
由于b' = a % b = a - a / b * b (注：这里的/是程序设计语言中的除法)
那么可以得到:
a'x + b'y = Gcd(a', b')  ===>
bx + (a - a / b * b)y = Gcd(a', b') = Gcd(a, b)  ===>
ay +b(x - a / b*y) = Gcd(a, b)
因此对于a和b而言，他们的相对应的p，q分别是 y和(x-a/b*y)

 

 

 

扩展欧几里德解一般二元一次方程a * x + b * y = c

从上面的过程可以看到，ax + by = gcd(a ,b)一定有解，再看一般形式ax + by=c，只有当且仅当c是gcd(a, b)的整数倍时才有解，否则无解。因为如果我们两边同除gcd(a, b)，会出现这样一种情况a'x + b'y=c / gcd(a, b)，如果c不是gcd(a,b)的整数倍，那么左边是整数，右边是分数，显然无解。

      至此，到网上搜很多博客都会发现这句话“ax+by=c的求解可以先求出ax  + by = gcd(a, b)，然后将x y扩大c / gcd(a,b)倍就可以了“，而我自己无法证明这句话的正确性，反而觉得这句话是错的，在前面已经得出ax + by = gcd(a, b)的解是：

      x = x0 + b / gcd(a, b) * k
      y = y0 - a / gcd(a, b) * k  （k为任意整数）

如果将x，y扩大c / gcd(a, b)倍后就成了

      x = x0 * c / gcd(a, b) + b * c * k / (gcd(a, b) ^ 2)
      y = y0 * c / gcd(a, b) - a * c * k / (gcd(a, b) ^ 2)（k为任意整数）

而我解出的通解是

      x = x0 * c / gcd(a, b) + b / gcd(a, b) * k
      y = y0 * c / gcd(a, b) - a / gcd(a, b) * k  （k为任意整数）

上面两组通解明显不等价，下面的一组通解包含了上面的一组通解，通过代入原方程ax + by = c验证，下面的一组通解才是正确的。所以，在求通解的过程中，只需要在特解x0 y0上扩大c / gcd(a, b)倍就可以了。

/*      Author: bcegkmqsw*/#include<stdio.h>#include<string.h>#include<stdlib.h>typedef __int64 LL;LL exgcd(LL a, LL b, LL &x, LL &y){    LL d, t;    if(b == 0)    {        x = 1, y = 0;        return a;    }    d = exgcd(b, a % b, x, y);    t = x, x = y, y = t - (a / b) * x;    return d;}int main(){    LL x, y, m, n, L, d, t, a, b, c;    while(scanf("%I64d%I64d%I64d%I64d%I64d", &x, &y, &m, &n, &L) != EOF)    {        a = m - n, b = L, c = y - x;        if(a < 0) a = -a, c = -c; // a > 0 => d > 0 => b > 0        d = exgcd(a, b, x, y);        if(m == n || c % d != 0)            printf("Impossible\n");        else        {            c /= d, t = c * x;            printf("%I64d\n", (t % b + b) % b);        }    }    return 0;}

不定方程ax+by=c 现在终于到了本文重点解二元一次不定方程。看起来扩展Euclid算法是不定方程的一种特殊情况实际上呢不定方程却是用Euclid算法解的。 对于不定方程ax+by=c设gcd(a,b)=d如果ax+by=c有解则d|c这也是许多奥数题的切入点。所以一旦d不是c的约数那么 ax+by=c一定无解。当d|c时先求出ax’+by’=d=gcd(a,b)的x'和y'由于已经有ax’+by’=d要求ax’+by’=c将整个式子同乘c/d倍即可。则x=x'*c/dy=y'*c/d。由上一段可知 只要ax+by=c有一个解它就有无数个解。 Euclid算法还可以求解同余方程ax≡b(mod m)及其最小x。这其实和不定方程ax+my=b没有区别。不定方程和同余方程一般都有范围限制这其实也很容易解决就不说了扩展欧几里得解一般二元一次方程 最小的正整数解不定方程ax+by=c

现在终于到了本文重点解二元一次不定方程。看起来扩展Euclid算法是不定方程的一种特殊情况实际上呢不定方程却是用Euclid算法解的。

 对于不定方程ax+by=c设gcd(a,b)=d如果ax+by=c有解则d|c这也是许多奥数题的切入点。所以一旦d不是c的约数那么 ax+by=c一定无解。

当d|c时先求出ax’+by’=d=gcd(a,b)的x'和y'由于已经有ax’+by’=d要求ax’+by’=c将整个式子同乘c/d倍即可。则x=x'*c/dy=y'*c/d。由上一段可知

只要ax+by=c有一个解它就有无数个解。 Euclid算法还可以求解同余方程ax≡b(mod m)及其最小x。这其实和不定方程ax+my=b没有区别。

不定方程和同余方程一般都有范围限制这其实也很容易解决就不说了

 http://blog.csdn.net/zjsxzjb/article/details/6262667  这个blog写的不错

 

 

Problem 1209

科大的浪漫

Time Limit: 1000ms
Memory Limit: 65536kb

Description

科大的浪漫，是春天微茫的细雨，
男孩和女孩同撑伞，
——三教上自习；
科大的浪漫，是夏季傍晚的凉意，
男孩和女孩肩并肩，
——讨论物化题；
科大的浪漫，是秋日晴朗的天际，
男孩和女孩面对面，
——做数学分析；
科大的浪漫，是冬天绵软的雪地，
男孩和女孩手牵手，
——一起去答疑。
（转自瀚海星云）

男孩和女孩正在三教上自习。他们有着共同的爱好，却有着不同的习惯。
男孩每自习a分钟会起身休息b分钟，而女孩每自习c分钟会起身休息d分钟。
如果他们同时起身休息，那么他们就会情不自禁地结束自习，收拾书包一同去吃饭。
男孩和女孩同时开始自习，他们想知道，多久以后他们才能结束自习？

Input

输入包含多组数据。
每组数据占一行，包括四个整数a, b, c, d (1≤a, b, c, d≤10⁹)
输入数据以a=b=c=d=0结束。

Output

对每组数据，输出多少分钟后男孩和女孩才能结束自习。如果他们永远不会结束自习，那么输出”Never”。
请注意结果可能会超过32位整型范围。虽然这很荒谬，但是他们的确如此强悍。这就是爱的力量。

Sample Input

2 1 3 11 1 2 21 1 1 10 0 0 0

Sample Output

11Never1

 

 

 

 

高效方法 扩展欧几里德算法
#include<stdio.h>
typedef struct
{
 long long d;
 long long x;
 long long y;
}Solve;
Solve exgcd(long long m,long long n)
{
 Solve s1={m,1,0},s2;
 if(n==0) return s1;
 s1=exgcd(n,m%n);
 s2.d=s1.d;
 s2.x=s1.y;
 s2.y=s1.x-m/n*s1.y;
 return s2;
}
int main(void)
{
 long long a,b,c,d;
 Solve solve;
 while(scanf("%lld%lld%lld%lld",&a,&b,&c,&d) && a)
 {
  solve=exgcd(a+b,c+d);
  if((c-a)%solve.d!=0)
  {
   printf("Never\n");
   continue;
  }
  solve.y*=(a-c)/solve.d;
  solve.y=solve.y%((a+b)/solve.d);
  if(solve.y<0)
   solve.y+=(a+b)/solve.d;
  printf("%lld\n",solve.y*(c+d)+c);
 }
 return 0;
}

 

 

#include <iostream>
#include <cstdio>
using namespace std;

long long gcd(long long a, long long b)
{
 if(b==0)
  return a;
 return gcd(b, a%b);
}

void egcd(long long a, long long b, long long c, long long & x, long long & y)
{
 if(b==0)
 {
  x=c/a;
  y=0;
 }
 else
 {
  long long xx, yy;
  egcd(b, a%b, c, xx, yy);
  x=-yy;
  y=-xx-a/b*yy;
 }
}

int main()
{
 long long a,b,c,d;
 while(cin>>a>>b>>c>>d, a)
 {
  long long c1=a+b, c2=c+d, c3=b-d, g;
  g=gcd(c1, c2);
  if(c3%g!=0)
  {
   cout<<"Never"<<endl;
   continue;
  }
  long long x, y;
  egcd(c1, c2, c3, x, y);
  if(x<=0)
  {
   x+=(-x)/(c2/g)*(c2/g)+c2/g;
  }
  else
  {
   x=x%(c2/g);
  }
  cout<<x*(a+b)-b<<endl;
 }
 return 0;
}