扩展欧几里德和模线性性方程 定理一:设 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.html http://www.cppblog.com/RyanWang/archive/2009/07/19/90525.html http://blog.163.com/kakarrot@yeah/blog/static/1201159252010413115321480/ 欧几里德算法(辗转相除法)是计算两个数最大公约数的传统算法,它无论从理论还是从效率上都是很好的。但是它有一个致命的缺陷,这个缺陷只有在大素数时才会显现出来。 Stein算法: 1、如果A=0,B是最大公约数,算法结束 2、如果B=0,A是最大公约数,算法结束 3、设置A1=A、B1=B和C1=1 4、如果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=Cn 8、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)); }