gcd及扩展gcd

其实这题用欧几里德扩展原理可以很快的解决，先来看下什么是欧几里德扩展原理：

　　欧几里德算法又称辗转相除法，用于计算两个整数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)的公约数是一样的，其最大公约数也必然相等，得证

　　欧几里德算法就是根据这个原理来做的，其算法用C++语言描述为：

　　int Gcd(int a, int b)
　　{
　　if(b == 0)
　　return a;
return Gcd(b, a % b);
　　}

　　当然你也可以写成迭代形式：

　　int Gcd(int a, int b)
　　{
　　while(b != 0)
　　{
　　int r = b;
　　 b = a % b;
　　 a = r;
　　}
　　return a;
　　}

　　本质上都是用的上面那个原理。

　　补充: 扩展欧几里德算法是用来在已知a, b求解一组x，y使得a*x+b*y=Gcd(a,b)(解一定存在，根据数论中的相关定理)。扩展欧几里德常用在求解模线性方程及方程组中。下面是一个使

用C++的实现：

　　int exGcd(int a, int b, int &x, int &y)
　　{
　　if(b == 0)
　　{
　　x = 1;
　　y = 0;
　　 return a;
　　}
　　int r = exGcd(b, a % b, x, y);
　　int t = x;
　　x = y;
　　y = t - a / b * y;
　　 return r;
　　}

把这个实现和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).

求解 x，y的方法的理解

　　设 a>b。

　　1，显然当 b=0，gcd（a，b）=a。此时 x=1，y=0；

　　2，ab<>0 时

　　设 ax1+by1=gcd(a,b);

　　bx2+(a mod b)y2=gcd(b,a mod b);

　　根据朴素的欧几里德原理有 gcd(a,b)=gcd(b,a mod b);

　　则:ax1+by1=bx2+(a mod b)y2;

　　即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+bx2-(a/b)*by2;

　　根据恒等定理得：x1=y2; y1=x2-(a/b)*y2;

　　这样我们就得到了求解 x1,y1 的方法：x1，y1 的值基于 x2，y2.

　　上面的思想是以递归定义的，因为 gcd 不断的递归求解一定会有个时候 b=0，所以递归可以

　结束。

　　在网上看了很多关于不定方程方程求解的问题，可都没有说全，都只说了一部分，看了好多之后才真正弄清楚不定方程的求解全过程，步骤如下：

　　求a * x + b * y = c的整数解。

　　1、先计算Gcd(a,b)，若n不能被Gcd(a,b)整除，则方程无整数解；否则，在方程两边同时除以Gcd(a,b)，得到新的不定方程a' * x + b' * y = c'，此时Gcd(a',b')=1;

2、利用上面所说的欧几里德算法求出方程a' * x + b' * y = 1的一组整数解x0,y0，则c' * x0,c' * y0是方程a' * x + b' * y = c'的一组整数解；

　　3、根据数论中的相关定理，可得方程a' * x + b' * y = c'的所有整数解为：

　　　　其实我们求得的解只是一组，

　　　　a*x0+lcm(a,b)+b*y0-lcm(a,b)=1;

　　　　a*x +b*y =1;

　　　　x=x0+b/gcd(a,b);y=y0-a/gcd(a,b);

　　　　a/gcd(a,b)*x'+b/gcd(a,b)*y'=c/gcd(a,b);

　　　　x'=c/gcd(a,b)*x0+b/gcd(a,b);y'=c/gcd(a,b)*y0-a/gcd(a,b);

x = c' * x0 + b' * t
y = c' * y0 - a' * t
(t为整数)

　　　 上面的解也就是a * x + b * y = n 的全部整数解。

欧几里德算法：

1 long long  gcd(long long  a,long long  b)
2 {
3     if(b==0)
4     {
5         return a;
6 
7     }
8     else return gcd(b,a%b);
9 }

扩展欧几里德：

 1 long long  exgcd(long long  a,long long  b,long long  &x,long long  &y)//返回值为最大公约数
 2 {
 3     if(b==0)
 4     {
 5         x=1;
 6         y=0;
 7          return  a;//最大公约数
 8     }
 9     ll d=exgcd(b,a%b,x,y);
10     long long  t=x;
11     x=y;
12     y=t-(a/b)*y;
13     return d;
14 }

扩展欧几里德算法，主要用于求解模线性方程 ax%b=c%b

则可化为 ax - by = c

这就可以用扩展欧几里德求

对于任意的方程 ax+by=c 求其整数解，若 c % gcd(a,b) ==0 则存在整数解，否则不存在整数解

首先我们先利用扩展欧几里德求 ax+by = gcd（a,b）；

然后得到的解 x *（c/gcd(a,b)）,则为 ax+by=c 的解

注意的地方：

1：扩展欧几里德求的是 ax+by = gcd（a,b）的解

2：若 x0 y0 ，是方程的一个解则，x0 +b ，y0 -a 也是方程的解

3：求解的到的 x 值只是方程的一个解，可能为正可能为负，，对于一般我们要将其转化为整数（x%mod+mod）%mod （x 可能为负的很大，所以先%mod ）

4：扩展欧几里德函数的参数，不用考虑其正负如求解 ax - by = c 直接 exgcd（a,b,x,y）;

hdu 1576 A/B

http://acm.hdu.edu.cn/showproblem.php?pid=1576

要求(A/B)%9973，但由于A很大，我们只给出n(n=A%9973)(我们给定的A必能被B整除，且gcd(B,9973) = 1)。

设 a/b = x; a=bx; bx%9973=n

bx-9973*y=n

扩展欧几里德做就行了

View Code

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<iostream>
 4 #include<algorithm>
 5 #include<set>
 6 #include<map>
 7 #include<queue>
 8 #include<vector>
 9 #include<string>
10 #define Min(a,b) a<b?a:b
11 #define Max(a,b) a>b?a:b
12 #define CL(a,num) memset(a,num,sizeof(num));
13 #define maxn  2100
14 #define inf 9999999
15 #define mx 1<<60
16 #define  mod   9973
17 using namespace std;
18 void exgcd(int a,int b,int &x,int &y)
19 {
20     if(b==0)
21     {
22         x=1;
23         y=0;
24         return ;
25     }
26     else
27     {
28         exgcd(b,a%b,x,y);
29         int t = x;
30         x=y;
31         y=t-(a/b)*y;
32     }
33 }
34 int main()
35 {
36     int t,x,y,n,b;
37     scanf("%d",&t);
38     while(t--)
39     {
40         scanf("%d%d",&n,&b);
41         exgcd( b,mod,x,y);
42         x*=n;
43         printf("%d\n",(x%mod+mod)%mod);
44     }
45 }

hdu 2669 Romantic

http://acm.hdu.edu.cn/showproblem.php?pid=2669

题意：输入 a b 求解 ax + by =1 的x 的最小正数解

View Code

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<iostream>
 4 #include<algorithm>
 5 #include<set>
 6 #include<map>
 7 #include<queue>
 8 #include<vector>
 9 #include<string>
10 #define Min(a,b) a<b?a:b
11 #define Max(a,b) a>b?a:b
12 #define CL(a,num) memset(a,num,sizeof(num));
13 #define maxn  2100
14 #define inf 9999999
15 #define mx 1<<60
16 #define  mod   9973
17 using namespace std;
18 int gcd(int a,int b)
19 {
20     if(b==0)
21     {
22         return a;
23 
24     }
25     else return gcd(b,a%b);
26 }
27 void  exgcd(int a,int b,int &x,int &y)
28 {
29     if(b==0)
30     {
31         x=1;
32         y=0;
33          return ;
34     }
35     exgcd(b,a%b,x,y);
36     int t=x;
37     x=y;
38     y=t-(a/b)*y;
39 }
40 int main()
41 {
42     int a,b,x,y;
43     while(scanf("%d%d",&a,&b)!=EOF)
44     {
45          int k=gcd(a,b);
46          if(1%k)printf("sorry\n");
47          else
48          {
49              exgcd(a,b,x,y);
50              x=x*(1/k);
51              y=y*(1/k);
52              while(x<0) x+=b,y-=a;
53 
54              printf("%d %d\n",x,y);
55          }
56     }
57 }

poj 2115 C Looooops

http://poj.org/problem?id=2115

求解 a + cy= b % 2 的k次幂

View Code

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<iostream>
 4 #include<algorithm>
 5 #include<set>
 6 #include<map>
 7 #include<queue>
 8 #include<vector>
 9 #include<string>
10 #define Min(a,b) a<b?a:b
11 #define Max(a,b) a>b?a:b
12 #define CL(a,num) memset(a,num,sizeof(num));
13 #define maxn  2100
14 #define inf 9999999
15 #define mx 1<<60
16 #define  mod   9973
17 
18 using namespace std;
19 long long  gcd(long long  a,long long  b)
20 {
21     if(b==0)
22     {
23         return a;
24 
25     }
26     else return gcd(b,a%b);
27 }
28 void  exgcd(long long  a,long long  b,long long  &x,long long  &y)
29 {
30     if(b==0)
31     {
32         x=1;
33         y=0;
34          return ;
35     }
36     exgcd(b,a%b,x,y);
37     long long  t=x;
38     x=y;
39     y=t-(a/b)*y;
40 }
41 int main()
42 {
43     long long  a,b,c,y,x;
44     long long  k;
45     while(scanf("%I64d %I64d %I64d %I64d",&a,&b,&c,&k),a+b+c+k)
46     {
47          k = (long long )1<<k;
48          long long  m=gcd(c,k);
49 
50          long long  s=b-a;
51          if(s%m)printf("FOREVER\n");
52          else
53          {
54               exgcd(c,k,x,y);
55               x*=s/m;
56               k=k/m;
57               if(x<0)
58                 x=x%k+k;
59 
60               printf("%lld\n",x%k);
61 
62          }
63 
64 
65     }
66 }

poj 1061 青蛙的约会

http://poj.org/problem?id=1061

View Code

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<iostream>
 4 #include<algorithm>
 5 #include<set>
 6 #include<map>
 7 #include<queue>
 8 #include<vector>
 9 #include<string>
10 #define Min(a,b) a<b?a:b
11 #define Max(a,b) a>b?a:b
12 #define CL(a,num) memset(a,num,sizeof(num));
13 #define maxn  2100
14 #define inf 9999999
15 #define mx 1<<60
16 #define  mod   9973
17 #define ll  long long
18 using namespace std;
19 long long  gcd(long long  a,long long  b)
20 {
21     if(b==0)
22     {
23         return a;
24 
25     }
26     else return gcd(b,a%b);
27 }
28 long long  exgcd(long long  a,long long  b,long long  &x,long long  &y)
29 {
30     if(b==0)
31     {
32         x=1;
33         y=0;
34          return  a;
35     }
36     ll d=exgcd(b,a%b,x,y);
37     long long  t=x;
38     x=y;
39     y=t-(a/b)*y;
40     return d;
41 }
42 int main()
43 {
44 
45     ll x ,y ,n ,m ,l,k;
46     while(scanf("%lld %lld %lld %lld %lld",&x,&y,&m,&n,&l)!=EOF)
47     {
48 
49           ll a = (m-n);
50           ll b = l;
51           ll c = ( y- x);
52           ll d = exgcd(a,b,x,y);
53 
54           if( c % d|| n==m)printf("Impossible\n");
55           else
56           {
57                c=c/d;
58 
59 
60 
61               x=((x*c)%b+b)%b;// 变为最小正数，x 有可能 < 0
62 
63               
64 
65                printf("%lld\n",x);
66           }
67     }
68 }