poj2115 扩展欧几里德算法小结

原创网址：http://blog.csdn.net/non_cease/article/details/7364092

以下是学习扩展欧几里德算法的资料的整合，有的转自别处，如百度百科。

扩展欧几里德算法源于欧几里德算法。

欧几里德算法：gcd（a，b）= gcd（b，a%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)的公约数是一样的，其最大公约数也必然相等，得证

拓展欧几里德是用来求二元一次不定方程a*x+b*y=gcd(a,b),---我们这样想 对于a'=b, b'=a%b=a-(a/b*b),

 a*x+b*y=gcd(a,b)

a'*x0+b'*y0=gcd(a',b')=gcd(a,b)=a*x+b*y

b*x0+(a-a/b*b)y0=a*x+b*y

a*y0+b*(x0-a/b*y0)=a*x+b*y

所以：  x=y0   ; y=x0-a/b*y0

扩展欧几里德算法如下：

[cpp]  view plain  copy

1. __int64 x, y; //x,y是a*x+b*y=gcd(a,b)的解  
2.   
3. __int64 Extended_Euclid(__int64 a, __int64 b) {//返回值为a，b最大公约数  
4.       if (b == 0) {  
5.          x = 1;  
6.          y = 0;  
7.          return a;  
8.       }  
9.       __int64 m = Extended_Euclid(b, a % b);  
10.       __int64 tmp = x;  
11.       x = y;  
12.       y = tmp - a / b * y;  
13.       return m;  
14. }  

还有就是纠结过一点为什么递归到最后x=1，y=0，因为最后是求方程a*x+0*y=a（gcd（a，0）==a），此时只要满足a=1，y为任意值就行，对于每一个y值，最后求出的都是方程的一组解。

例如：求20*x+12*y=4       注意下面是根据上面求解得到的递推公式x=y0   y=x0-a/b*y0由下往上递推得到的

a   b      ：      20                      12

x   y      ：      3*y0-x0              2*x0-5*y0

a   b      ：      12                      8

x   y      ：      x0-2*y0              3*y0-x0

a   b      ：      8                        4

x   y        ：     y0                      x0-2*y0

a   b      ：       4                        0

x   y       ：      x0                       y0

所以对于x0 = 1，y0没取一个值，就将得到原方程20*x+12*y=4的一组解，简单起见，y0取0即可。

其实算法就是用递归的方式完成递推，用4*x+0*y=4的解通过公式x=y0 ; y=x0-a/b*y0递推求出8*x+4*y=4,

一直递推直到求出20*x+12*y=4的解。

但是往往我们需要求解的方程并不是a*x+b*y=gcd(a,b)而是更常规的a*x+b*y=n

然后可以通过对a*x+b*y=gcd(a,b)方程的求解，转化成为求二元一次不定方程a*x+b*y=n

若gcd（a，b）不能整除n，这个方程无整数解，反之，若解得a*x+b*y=gcd(a,b)的解为x0，y0，

则方程两边同乘n再除以gcd（a，b）得a*（n/gcd(a,b)*x0)+b*(n/gcd(a,b)*y0)=n

所以方程解为x=n/gcd(a,b)*x0，y=n/gcd(a,b)*y0。

更通常的是：我们需要求解方程的最小整数解

若我们已经求得x0,y0为方程中x的一组特解，那么

x=x0+b/gcd(a,b)*t，y=y0-a/gcd(a,b)*t(t为任意整数）也为方程的解

且b/gcd(a,b)，a/gcd(a,b)分别为x，y的解的最小间距，所以x在0~b/gcd(a,b)区间有且仅有一个解，

同理y在0~a/gcd(a,b)同样有且仅有一个解，这个解即为我们所需求的最小正整数解。

为什么b/gcd(a,b)，a/gcd(a,b)分别为x，y的解的最小间距？

解：假设c为x的解的最小间距，此时d为y的解的间距，所以x=x0+c*t，y=y0-d*t（x0，y0为一组特解，t为任意整数）

      带入方程得：a*x0+a*c*t+b*y0-b*d*t=n,因为a*x0+b*y0=n，所以a*c*t-b*d*t=0，t不等于0时，a*c=b*d

      因为a,b,c,d都为正整数，所以用最小的c，d，使得等式成立，ac，bd就应该等于a，b的最小公倍数a*b/gcd(a,b),

      所以c=b/gcd(a,b)，d就等于a/gcd(a,b)。

所以，若最后所求解要求x为最小整数，那么x=(x0%（b/gcd(a,b))+b/gcd(a,b))%(b/gcd(a,b))即为x的最小整数解。

x0%（b/gcd(a,b))使解落到区间-b/gcd(a,b)~b/gcd(a,b)，再加上b/gcd(a,b)使解在区间0~2*b/gcd(a,b)，

再模上b/gcd(a,b)，则得到最小整数解（注意b/gcd(a,b)为解的最小距离，重要）

poj2115

[cpp]  view plain  copy

1. //poj2115  
2.   
3. #include <iostream>  
4. #include <cstdio>  
5. using namespace std;  
6.   
7. __int64 x, y;  
8.   
9. __int64 Extended_Euclid(__int64 a, __int64 b) {  
10.       if (b == 0) {  
11.          x = 1;  
12.          y = 0;  
13.          return a;  
14.       }  
15.       __int64 m = Extended_Euclid(b, a % b);  
16.       __int64 tmp = x;  
17.       x = y;  
18.       y = tmp - a / b * y;  
19.       return m;  
20. }  
21.   
22. int main()  
23. {  
24.     __int64 a, b, c, k;  
25.   
26.     while (scanf ("%I64d%I64d%I64d%I64d", &a, &b, &c, &k)) {  
27.         if (!a && !b && !c && !k) break;  
28.         __int64 d = (__int64)1<<k; //此处先讲1转化为64位，否则系统默认1为int型，右移32位会溢出  
29.         __int64 e = Extended_Euclid(c, d);  
30.         if ((b-a) % e) { printf ("FOREVER\n"); continue; }  
31.         __int64 t = d / e;  
32.         x = (x  * (b - a) / e % t + t) % t;  
33.         printf ("%I64d\n", x);  
34.     }  
35.     return 0;  
36. }