线性方程:设a和b是两个整数,g = gcd(a,b)是a和b的最大公约数。求满足方程 a*x + b*y = g 的整数解x和y。
递归版:扩张欧几里德
在用欧几里德算法算a和b的最大公约数时,我们依次得到:
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)为0时r(n-1)便是a和b的最大公约数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 }