模线性方程(递归版+迭代版)& 扩展欧几里德

  线性方程:设ab是两个整数,g = gcd(a,b)ab的最大公约数。求满足方程 a*x + b*y = g 的整数解xy

  递归版:扩张欧几里德

  在用欧几里德算法算ab的最大公约数时,我们依次得到:

  a = q(1) * b + r(1)

  b = q(2) * r(1) + r(2)

  r1 = q(3) * r(2) + r(3)

    ......

  r(n-2) = q(n) * r(n-1) + r(n)

  当r(n)0r(n-1)便是ab的最大公约数g。我们在Euclid递归返回的过程中,构造解。

  设第k-1层有:

    r(k-3) = q(k-1) * r(k-2) + r(k-1)

  第k层有:

    r(k-2) = q(k) * r(k-1) + r(k)

  递归时,设第k层的解已经构造好了,即我们已经知道了:x(k)y(k),满足:x(k) * r(k-2) + y(k) * r(k-1) = g.

  则第k-1层为:x(k-1) * r(k-3) + y(k-1) * r(k-2) = g,即x(k-1)*( q(k-1)*r(k-2) + r(k-1) ) + y(k-1)*r(k-2) = g。整理得:x(k-1)*r(k-1) + (x(k-1)*q(k-1) + y(k-1))*r(k-2) = g

  得到系数的转移方程:

    x(k-1) = y(k)

    y(k-1) = x(k) - y(k)*q(k-1)

  扩展欧几里德代码:

 1 //zzy2012.7.17
 2 #include<cstdio>
 3 #include<iostream>
 4 
 5 using namespace std;
 6 
 7 int Extended_Euclid(int a,int b,int &x,int &y){
 8     int q,temp,g;
 9     if(b==0){
10         x = 1;
11         y = 0;
12         return a;
13     }
14     g = Extended_Euclid(b,a%b,x,y);
15     q = a/b;
16     temp = x;
17     x = y;
18     y = temp - q*y;
19     return g;
20 
21 }
22 
23 int main()
24 {
25     int a,b,x,y,g;
26     cin>>a>>b;
27     g = Extended_Euclid(a,b,x,y);
28     cout<<a<<'*'<<x<<'+'<<b<<'*'<<y<<'='<<g<<endl;
29     return 0;
30 }

 

  迭代版

  考察第k-1层:

    r(k-3) = q(k-1) * r(k-2) + r(k-1)

  得r(k-1) = r(k-3) - q(k-1) * r(k-2),得系数递推关系:

    x(k) = x(k-2) - q(k) * x(k-1)

    y(k) = y(k-2) - q(k) * y(k-1)

  迭代版代码:

 1 //zzy2012.2.21
 2 //求解形如a*x+b*y=gcd(a,b)的线性方程
 3 #include<cstdio>
 4 #include<iostream>
 5 
 6 using namespace std;
 7 
 8 typedef struct
 9 {
10     int x,y;
11 }node;
12 
13 node r[3];
14 
15 int solve(int a,int b)
16 {
17     int f,q,t;//f指向r元素的下标,q每次迭代算出的商,t保存余数,r[]数组保存每次计算出的系数
18     f=2;
19     r[0].x=1;
20     r[0].y=0;
21     r[1].x=0;
22     r[1].y=1;
23     while(b!=0)
24     {
25         q=a/b;
26         t=a%b;
27         a=b;
28         b=t;
29         r[f%3].x=r[(f-2)%3].x-q*r[(f-1)%3].x;
30         r[f%3].y=r[(f-2)%3].y-q*r[(f-1)%3].y;
31         f++;
32     }
33     return (f-2)%3;
34 }
35 
36 int main()
37 {
38     int a,b,f;
39     cin>>a>>b;
40     if(a!=0 || b!=0)
41     {
42         f=solve(a,b);
43         printf("%d %d\n",r[f].x,r[f].y);
44     }
45     return 0;
46 }

 

   同余式a*x = c (mod m)可以转换为求线性方程。

  代码:

 1 /*-----ЭЌгрЪН-----2011.7.10-----Library-----zzy-----*/
 2 #include<cstdio>
 3 #include<cmath>
 4 
 5 int arg[2];
 6 
 7 void ini(int a,int b)
 8 {
 9     arg[0]=0;
10     arg[1]=1;
11     arg[2]=1;
12     arg[3]=-a/b;
13 }
14 
15 int lineeq(int a,int b)
16 {
17     int x,q;
18     if(b==0)
19         return a;
20     q=-a/b;
21     x=arg[0]-q*arg[1];
22     arg[0]=arg[1];
23     arg[1]=x;
24     return lineeq(b,a%b);
25 }
26 
27 int main()
28 {
29     int a,c,m,x,g,i;
30     scanf("%d %d %d",&a,&c,&m);
31     if(m==0)
32     {
33         printf("illegal m.\n");
34         return 0;
35     }
36     if(a==0)
37     {
38         if(c%m==0)
39             printf("x can be any integer.\n");
40         else
41             printf("no solution.\n");
42         return 0;
43     }
44     ini(a,-m);
45     g=lineeq(-m,a%(-m));
46     g=abs(g);
47     if(c%g!=0)
48     {
49         printf("no solution.\n");
50         return 0;
51     }
52     x=arg[0]*c/g;
53     for(i=0; i<g; i++)
54     {
55         printf("%d + k * %d\n",x+i*m/g,m);
56     }
57     return 0;
58 }

 

posted on 2012-07-17 20:09  Lattexiaoyu  阅读(315)  评论(11编辑  收藏  举报

导航