扩展欧几里德和模线性性方程(总结)

扩展欧几里德和模线性性方程定理一：设 a, b, c 是三个不全为 0 的整数。如果 a= b* q+ c， 其中 q 是整数，则 (a,b)= (b,c).证明：设 d1= (a,b)， d2= (b,c). 有 d1| a+ (-q)* b，所以 d1| c，因而 d1 为 b, c 的公约数，所以 d1<= d2。同理，d2 是 a, b 的公因数，有 d2<= d1，故 d1== d2，定理成立。由定理一，可得到求最大公约数的欧几里德算法：int gcd( int a, int b ){ return !b? a: gcd( b, a% b ); } 有一结论: 对于正整数 a, b，(a,b)= d，则 d 一定可以表示成 d= ax+ by， x, y 为整数。故 gcd(a,b) 可以写成如下形式，gcd(a,b)= ax+ by, x, y 为整数，用扩展欧几里德算法可以求出 x, y。int extend_gcd( int a, int b, int &x, int &y ){ if( !b ){ x= 1, y= 0; return a; } int g= extend_gcd( b, a% b, x, y ); int t= x; x= y; y= t- a/ b* y; return g;}以上算法中，b== 0 时，返回 a, x0= 1, y0= 0, 满足 a= a*x0+ b* y0。b!= 0 时，算法先计算出 d= gcd( b, a% b ) 和 x, y，并且满足：d= b* x0+ (a% b )* y0所以， gcd(a,b)= gcd(b,a%b)= d= b* x0+ (a% b)* y0= b* x0+ ( a- [a/b]* b)* y0 = a* y0+ b* ( x0- [a/b]* y0 )= y0* a+ (x0- [a/b]* y0)* b因此更新 x= y0, y= (x0+ [a/b]* y0) 可以满足等式 d= ax+ by，这样就证明算法的正确性。定理二：同余方程 ax=b (mod n) 对于未知数 x 有解，仅当 gcd(a,n) | b。且方程有解时，方程有 gcd(a,n) 个解。证明比较复杂。不难得出，由这个定理可以推出以上的那个结论 gcd(a,b)= ax+ by 这样表示。特别的，如果 gcd(a,n)== 1，则方程只有唯一解。在这种情况下，如果 b== 1，同余方程就是 ax=1 (mod n ),gcd(a,n)= 1。这时称求出的 x 为 a 的对模 n 乘法的逆元。对于同余方程 ax= 1(mod n )， gcd(a,n)= 1 的求解就是求解方程ax+ ny= 1，x, y 为整数。这个可用扩展欧几里德算法求出，原同余方程的唯一解就是用扩展欧几里德算法得出的 x 。定理三：设 d= gcd(a,n)，假如整数 x 和 y，满足 d= ax+ ny(用扩展欧几里德得出)。如果 d| b，则方程ax= b( mod n ) 的一个解为 x0= x* (b/ d ) mod n，且方程的 d 个解分别为 xi= x0+ i* (n/ d ) {i= 0... d-1}。求解方程 ax= b (mod n ) 相当于求解方程 ax+ ny= b, (x, y为整数)a* x0+ n* y0= d， 方程两边乘以 b/ d，(因为 d|b，所以能够整除)，得到 a* x0* b/ d+ n* y0* b/ d= b。所以 x= x0* b/ d，y= y0* b/ d 为 ax+ ny= b 的一个解，所以 x= x0* b/ d 为 ax= b (mod n ) 的解。a* xi mod n= a( x0+ i* n/ d ) mod n = ( a* x0+ a* i* n/ d ) mod n = a* x0 mod n ( 因为 d| a ) = b。 int Modular_Linear(int a,int b,int n)//ax=b(mod n){int d,x,y,x0,i;d=extend_gcd(a,n,x,y);if(b%d==0){x0=(x*(b/d))%n;if(x0<n)x0+=n;for(i=0;i<d;i++)cout<<(x0+i*n/d)%n<<endl;}return 0;}求解模线性方程组(中国余数定理)语法：result=Modular_Expoent(int a,int b,int n); 参数: B[]、W[]： a=B[] (mod W[]) 的对应参数返回值：a 的值注意：其中W[],B[]已知，W[i]>0且W[i]与W[j]互质, 求a源程序：int ext_euclid(int a,int b,int &x,int &y) //求gcd(a,b)=ax+by{ int t,d; if (b==0) {x=1;y=0;return a;} d=ext_euclid(b,a %b,x,y); t=x; x=y; y=t-a/b*y; return d;}int China(int B[],int W[],int k){ int i; int d,x,y,a=0,m,n=1; for (i=0;i<k;i++) n*=W[i]; for (i=0;i<k;i++) { m=n/W[i]; d=ext_euclid(W[i],m,x,y); a=(a+y*m*B[i])%n; } if (a>0) return a; else return(a+n);}pku 1061 青蛙的约会由题意，我们设要跳 p 次才能相遇， 则有方程 (m- n)* p= y- x (mod L )，就是求解一个同余方程。代码：代码#include <iostream>#include <cstdio>using namespace std;typedef long long INT;INT x, y, m, n, l;INT gcd( INT a, INT b ){ return b== 0? a: gcd( b, a% b ); }INT extend_gcd( INT a, INT b, INT &x, INT &y ){ if( !b ){ x= 1, y= 0; return a; } INT g= extend_gcd( b, a% b, x, y ); INT t= x; x= y; y= t- a/ b* y; return g;}int main(){ while( cin>> x>> y>> m>> n>> l ){ INT g= gcd( m- n, l ); if( (y- x)% g!= 0 ){ cout<< "Impossible" << endl; continue; } INT p, q; extend_gcd( m- n, l, p, q ); INT ans= p* (y- x)/ g; cout << ( ans% l+ l )% l << endl; } return 0;}附加：TOJ2869 http://acm.tju.edu.cn/toj/showp2869.html大意：已知a,b,c,n,求出一组x,y,z（满足条件的任意一组）使ax+by+cz=n.若无解则输出-1.直接给代码：#include<iostream>using namespace std;int m1,n1;int ext_gcd(int a,int b){ if(b==0) { m1 = 1; n1 = 0; return a; } else { int dab = ext_gcd(b,a%b); int tm = m1; m1 = n1; n1 = tm - a/b*n1; return dab; }}int main(){ int a,b,c,n; while(cin>>a>>b>>c>>n) { int dab = ext_gcd(a,b); int ma,mb,mc; ma = m1; mb = n1; int dabc = ext_gcd(dab,c); if(n%dabc != 0){cout<<-1<<endl; continue;} mc = n/dabc*n1; ma *= (n/dabc*m1); mb *= (n/dabc*m1); int mabc = (a * b * c)/dabc; while(ma<0) { ma += mabc/a; mc -= mabc/c; } while(mb<0) { mb += mabc/b; mc -= mabc/c; } while(mc>0) { mc-=mabc/c; ma+=mabc/a; } mc = -mc; cout<<ma<<' '<<mb<<' '<<mc<<endl; } return 0;}hdu 1576 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)%9973。算法分析：扩展欧几里德算法，推导式为bx-9973y=n, 设二元一次不定方程ax+by=c有解，X,Y时它的一组解，则其所有解为 x=X+b/(a,b)*t y=Y+a/(a,b)*t t<-Z 那么现在的问题就是求出方程bx-9973y=n, 的任意一组解 x,y,G=9973/gcd(b,9973)则 (x%G+G)%G就是所求的答案 #include <iostream>#include <cstdio>#include <cstring>using namespace std;int mod = 9973;inline int fun(int a, int b, int & x, int & y){if(b == 0){x = 1;y = 0;return a;}int r = fun(b, a % b, x, y);int temp = x;x = y;y = temp - a / b * y;return r;}int main(){int t;scanf("%d", &t);while(t--){int n, b;int x, y;scanf("%d %d", &n, &b);int temp = fun(b, mod, x, y);x = x * (n / temp);//int ret = mod / temp; 这里不加也是AC的x = (x % mod + mod) % mod;printf("%d/n", x);}}TJU 3027 模线性方程组比较赤裸的模线性方程组模型了，给定的模都是两两互素的，可以用扩展欧几里德和孙子定理来解。第一次写这个，不知道健壮性如何，这个题目数据貌似很弱。。#include <cstdio>const int N = 8;void extended_gcd(int a, int b, int &x, int &y){ if (!b) { x = 1; y = 0; return; } extended_gcd(b, a % b, x, y); int tmp = x; x = y; y = tmp - a / b * y;}int modular_linear_equation_system(int b[N], int m[N], int n){ int M = 1, x, y, ret = 0; for (int i = 0; i < n; i++) M *= m[i]; for (int i = 0; i < n; i++) { extended_gcd(M / m[i], m[i], x, y); ret += M / m[i] * x * b[i]; } while (ret <= 0) ret += M; while (ret > M) ret -= M; return ret;}int main(){ int m[N], b[N], T, n, tmp, ans; const int num = 4, numc = 3; char str[N], hash[30] = " ABCDEFGHIJKLMNOPQRSTUVWXYZ "; scanf("%d", &T); while (T--) { scanf("%d", &n); for (int i = 0; i < num; i++) scanf("%d", &m[i]); for (int j = 0; j < n; j++) { scanf("%d", &tmp); for (int i = 1; i <= num; i++) { b[num-i] = tmp % 100; tmp /= 100; } ans = modular_linear_equation_system(b, m, num); for (int i = 1; i <= numc; i++) { str[numc-i] = hash[ans%100]; ans /= 100; } str[numc] = '/0'; if (j == n - 1) { int k = numc - 1; for (; k >= 0; k--) if (str[k] != ' ') break; str[k+1] = '/0'; } printf("%s", str); } putchar('/n'); } return 0;}感谢以下：http://hi.baidu.com/ve_bird/blog/item/cf585e3f04c82dcb9f3d6203.htmlhttp://www.cppblog.com/RyanWang/archive/2009/07/19/90525.htmlhttp://blog.163.com/kakarrot@yeah/blog/static/1201159252010413115321480/欧几里德算法(辗转相除法)是计算两个数最大公约数的传统算法，它无论从理论还是从效率上都是很好的。但是它有一个致命的缺陷，这个缺陷只有在大素数时才会显现出来。 Stein算法：1、如果A=0，B是最大公约数，算法结束2、如果B=0，A是最大公约数，算法结束3、设置A1=A、B1=B和C1=14、如果An和Bn都是偶数，则An+1=An/2，Bn+1=Bn/2，Cn+1=Cn*2(注意，乘2只要把整数左移一位即可，除2只要把整数右移一位即可)5、如果An是偶数，Bn不是偶数，则An+1=An/2，Bn+1=Bn，Cn+1=Cn (很显然啦，2不是奇数的约数) 6、如果Bn是偶数，An不是偶数，则Bn+1=Bn/2，An+1=An，Cn+1=Cn (很显然啦，2不是奇数的约数)7、如果An和Bn都不是偶数，则An+1=|An-Bn|，Bn+1=min(An,Bn)，Cn+1=Cn8、n加1，转1 考虑欧几里德算法，最恶劣的情况是，每次迭代a=2b-1,这样，迭代后，r=b-1。如果a小于2N，这样大约需要4N次迭代。而考虑Stein算法，每次迭代后，显然A(n+1)B(n+1)≤AnBn/2，最大迭代次数也不超过4N次。也就是说，迭代次数几乎是相等的。但是，需要注意的是，对于大素数，试商法将使每次迭代都更复杂，因此对于大素数Stein将更有优势。 //下面是根据该原理写的实现函数 unsigned MaxDivisor(unsigned a, unsigned b) { unsigned c = 0; while(1) { // 退出条件 if(a==0) return b << c; else if(b == 0) return a << c; // 为提高速度，避免用取模判断奇偶 if(!(a & 1) && !(b & 1)) //a,b 都是偶数 { a >>= 1; b >>= 1; ++c; } else if(!(a & 1) && (b & 1)) //a偶 b奇 { a >>= 1; } else if((a & 1) && !(b & 1)) //a奇 b偶 { b >>= 1; } else if((a & 1) && (b & 1)) //a,b都是奇数 { unsigned tmp = a>b?b:a; //取较小的一个 a = a>b?a-b:(b-a); //绝对差值 b = tmp; } }}int Stein(int a,int b) //最优的Stein算法{ if(a == 0)return b; if(b == 0)return a; if(a%2 == 0&&b == 0)return 2*gcd(a/2,b/2); else if(a%2 == 0)return gcd(a/2,b); else if(b%2 == 0)return gcd(a,b/2); else return gcd(abs(a - b),min(a,b));}