擴充歐幾裡德和模線性性方程(總結)

擴充歐幾裡德和模線性性方程<br />定理一：設 a, b, c 是三個不全為 0 的整數。如果 a= b* q+ c， 其中 q 是整數，則 (a,b)= (b,c).<br />證明：設 d1= (a,b)， d2= (b,c). 有 d1| a+ (-q)* b，所以 d1| c，因而 d1 為 b, c 的公約數，所以 d1<= d2。<br />同理，d2 是 a, b 的公因數，有 d2<= d1，故 d1== d2，定理成立。<br />由定理一，可得到求最大公約數的歐幾裡德演算法：<br />int gcd( int a, int b ){<br /> return !b? a: gcd( b, a% b );<br />}<br />有一結論: 對於正整數 a, b，(a,b)= d，則 d 一定可以表示成 d= ax+ by， x, y 為整數。<br />故 gcd(a,b) 可以寫成如下形式，gcd(a,b)= ax+ by, x, y 為整數，用擴充歐幾裡德演算法可以求出 x, y。<br />int extend_gcd( int a, int b, int &x, int &y ){<br /> if( !b ){<br /> x= 1, y= 0; return a;<br /> }<br /> int g= extend_gcd( b, a% b, x, y );<br /> int t= x; x= y; y= t- a/ b* y;<br /> return g;<br />}<br />以上演算法中，b== 0 時，返回 a, x0= 1, y0= 0, 滿足 a= a*x0+ b* y0。<br />b!= 0 時，演算法先計算出 d= gcd( b, a% b ) 和 x, y，並且滿足：<br />d= b* x0+ (a% b )* y0<br />所以， gcd(a,b)= gcd(b,a%b)= d= b* x0+ (a% b)* y0= b* x0+ ( a- [a/b]* b)* y0<br /> = a* y0+ b* ( x0- [a/b]* y0 )= y0* a+ (x0- [a/b]* y0)* b<br />因此更新 x= y0, y= (x0+ [a/b]* y0) 可以滿足等式 d= ax+ by，這樣就證明演算法的正確性。<br />定理二：同餘方程 ax=b (mod n) 對於未知數 x 有解，僅當 gcd(a,n) | b。且方程有解時，方程有 gcd(a,n) 個解。<br />證明比較複雜。不難得出，由這個定理可以推出以上的那個結論 gcd(a,b)= ax+ by 這樣表示。<br />特別的，如果 gcd(a,n)== 1，則方程只有唯一解。在這種情況下，如果 b== 1，同餘方程就是 ax=1 (mod n ),<br />gcd(a,n)= 1。這時稱求出的 x 為 a 的對模 n 乘法的逆元。<br />對於同餘方程 ax= 1(mod n )， gcd(a,n)= 1 的求解就是求解方程<br />ax+ ny= 1，x, y 為整數。這個可用擴充歐幾裡德演算法求出，原同餘方程的唯一解就是用擴充歐幾裡德演算法得出的 x 。<br />定理三：設 d= gcd(a,n)，假如整數 x 和 y，滿足 d= ax+ ny(用擴充歐幾裡德得出)。如果 d| b，則方程<br />ax= b( mod n ) 的一個解為 x0= x* (b/ d ) mod n，且方程的 d 個解分別為 xi= x0+ i* (n/ d ) {i= 0... d-1}。<br />求解方程 ax= b (mod n ) 相當於求解方程 ax+ ny= b, (x, y為整數)<br />a* x0+ n* y0= d， 方程兩邊乘以 b/ d，(因為 d|b，所以能夠整除)，得到 a* x0* b/ d+ n* y0* b/ d= b。<br />所以 x= x0* b/ d，y= y0* b/ d 為 ax+ ny= b 的一個解，所以 x= x0* b/ d 為 ax= b (mod n ) 的解。<br />a* xi mod n= a( x0+ i* n/ d ) mod n<br /> = ( a* x0+ a* i* n/ d ) mod n<br /> = a* x0 mod n ( 因為 d| a )<br /> = b。</p><p>int Modular_Linear(int a,int b,int n)//ax=b(mod n)<br />{<br />int d,x,y,x0,i;<br />d=extend_gcd(a,n,x,y);<br />if(b%d==0)<br />{<br />x0=(x*(b/d))%n;<br />if(x0<n)<br />x0+=n;<br />for(i=0;i<d;i++)<br />cout<<(x0+i*n/d)%n<<endl;<br />}<br />return 0;<br />}</p><p>求解模線性方程組(中國餘數定理)<br />文法：result=Modular_Expoent(int a,int b,int n);<br />參數: B[]、W[]： a=B[] (mod W[]) 的對應參數<br />傳回值：a 的值<br />注意：其中W[],B[]已知，W[i]>0且W[i]與W[j]互質, 求a<br />來源程式：<br />int ext_euclid(int a,int b,int &x,int &y) //求gcd(a,b)=ax+by<br />{<br /> int t,d;<br /> if (b==0) {x=1;y=0;return a;}<br /> d=ext_euclid(b,a %b,x,y);<br /> t=x;<br /> x=y;<br /> y=t-a/b*y;<br /> return d;<br />}<br />int China(int B[],int W[],int k)<br />{<br /> int i;<br /> int d,x,y,a=0,m,n=1;<br /> for (i=0;i<k;i++)<br /> n*=W[i];<br /> for (i=0;i<k;i++)<br /> {<br /> m=n/W[i];<br /> d=ext_euclid(W[i],m,x,y);<br /> a=(a+y*m*B[i])%n;<br /> }<br /> if (a>0) return a;<br /> else return(a+n);<br />}<br />pku 1061 青蛙的約會<br />由題意，我們設要跳 p 次才能相遇， 則有方程 (m- n)* p= y- x (mod L )，就是求解一個同餘方程。<br />代碼：<br />代碼<br />#include <iostream><br />#include <cstdio><br />using namespace std;<br />typedef long long INT;<br />INT x, y, m, n, l;<br />INT gcd( INT a, INT b ){<br /> return b== 0? a: gcd( b, a% b ); }<br />INT extend_gcd( INT a, INT b, INT &x, INT &y ){<br /> if( !b ){<br /> x= 1, y= 0; return a;<br /> }<br /> INT g= extend_gcd( b, a% b, x, y );<br /> INT t= x; x= y; y= t- a/ b* y;<br /> return g;<br />}<br />int main(){<br /> while( cin>> x>> y>> m>> n>> l ){<br /> INT g= gcd( m- n, l );<br /> if( (y- x)% g!= 0 ){<br /> cout<< "Impossible" << endl;<br /> continue;<br /> }</p><p> INT p, q;<br /> extend_gcd( m- n, l, p, q );<br /> INT ans= p* (y- x)/ g;</p><p> cout << ( ans% l+ l )% l << endl;<br /> }</p><p> return 0;<br />}<br />附加：TOJ2869 http://acm.tju.edu.cn/toj/showp2869.html<br />大意：已知a,b,c,n,求出一組x,y,z（滿足條件的任意一組）使ax+by+cz=n.<br />若無解則輸出-1.<br />直接給代碼：<br />#include<iostream><br />using namespace std;<br />int m1,n1;<br />int ext_gcd(int a,int b)<br />{<br /> if(b==0)<br /> {<br /> m1 = 1;<br /> n1 = 0;<br /> return a;<br /> }<br /> else<br /> {<br /> int dab = ext_gcd(b,a%b);<br /> int tm = m1;<br /> m1 = n1;<br /> n1 = tm - a/b*n1;<br /> return dab;<br /> }<br />}<br />int main()<br />{<br /> int a,b,c,n;<br /> while(cin>>a>>b>>c>>n)<br /> {<br /> int dab = ext_gcd(a,b);<br /> int ma,mb,mc;<br /> ma = m1;<br /> mb = n1;<br /> int dabc = ext_gcd(dab,c);<br /> if(n%dabc != 0){cout<<-1<<endl; continue;}<br /> mc = n/dabc*n1;<br /> ma *= (n/dabc*m1);<br /> mb *= (n/dabc*m1);<br /> int mabc = (a * b * c)/dabc;<br /> while(ma<0)<br /> {<br /> ma += mabc/a;<br /> mc -= mabc/c;<br /> }<br /> while(mb<0)<br /> {<br /> mb += mabc/b;<br /> mc -= mabc/c;<br /> }<br /> while(mc>0)<br /> {<br /> mc-=mabc/c;<br /> ma+=mabc/a;<br /> }<br /> mc = -mc;<br /> cout<<ma<<' '<<mb<<' '<<mc<<endl;<br /> }<br /> return 0;<br />}<br />hdu 1576 http://acm.hdu.edu.cn/showproblem.php?pid=1576<br />題目大意：要求(A/B)%9973，但由於A很大，我們只給出n(n=A%9973)(我們給定的A必能被B整除，且 gcd(B,9973) = 1)。 對應每組資料輸出(A/B)%9973。<br />演算法分析：擴充歐幾裡德演算法，推導式為bx-9973y=n,<br /> 設二元一次不定方程ax+by=c有解，X,Y時它的一組解，則其所有解為<br /> x=X+b/(a,b)*t<br /> y=Y+a/(a,b)*t t<-Z 那麼現在的問題就是求出方程bx-9973y=n, 的任意一組解 x,y,<br />G=9973/gcd(b,9973)<br />則 (x%G+G)%G就是所求的答案<br />#include <iostream><br />#include <cstdio><br />#include <cstring><br />using namespace std;<br />int mod = 9973;<br />inline int fun(int a, int b, int & x, int & y){<br />if(b == 0){<br />x = 1;<br />y = 0;<br />return a;<br />}<br />int r = fun(b, a % b, x, y);<br />int temp = x;<br />x = y;<br />y = temp - a / b * y;<br />return r;<br />}<br />int main(){<br />int t;<br />scanf("%d", &t);<br />while(t--){<br />int n, b;<br />int x, y;<br />scanf("%d %d", &n, &b);<br />int temp = fun(b, mod, x, y);<br />x = x * (n / temp);<br />//int ret = mod / temp; 這裡不加也是AC的<br />x = (x % mod + mod) % mod;<br />printf("%d/n", x);<br />}<br />}<br />TJU 3027 模線性方程組<br />比較赤裸的模線性方程組模型了，給定的模都是兩兩互素的，可以用擴充歐幾裡德和孫子定理來解。第一次寫這個，不知道健壯性如何，這個題目資料貌似很弱。。<br />#include <cstdio><br />const int N = 8;<br />void extended_gcd(int a, int b, int &x, int &y)<br />{<br /> if (!b)<br /> {<br /> x = 1;<br /> y = 0;<br /> return;<br /> }<br /> extended_gcd(b, a % b, x, y);<br /> int tmp = x;<br /> x = y;<br /> y = tmp - a / b * y;<br />}<br />int modular_linear_equation_system(int b[N], int m[N], int n)<br />{<br /> int M = 1, x, y, ret = 0;<br /> for (int i = 0; i < n; i++)<br /> M *= m[i];<br /> for (int i = 0; i < n; i++)<br /> {<br /> extended_gcd(M / m[i], m[i], x, y);<br /> ret += M / m[i] * x * b[i];<br /> }<br /> while (ret <= 0) ret += M;<br /> while (ret > M) ret -= M;<br /> return ret;<br />}<br />int main()<br />{<br /> int m[N], b[N], T, n, tmp, ans;<br /> const int num = 4, numc = 3;<br /> char str[N], hash[30] = " ABCDEFGHIJKLMNOPQRSTUVWXYZ ";<br /> scanf("%d", &T);<br /> while (T--)<br /> {<br /> scanf("%d", &n);<br /> for (int i = 0; i < num; i++)<br /> scanf("%d", &m[i]);<br /> for (int j = 0; j < n; j++)<br /> {<br /> scanf("%d", &tmp);<br /> for (int i = 1; i <= num; i++)<br /> {<br /> b[num-i] = tmp % 100;<br /> tmp /= 100;<br /> }<br /> ans = modular_linear_equation_system(b, m, num);<br /> for (int i = 1; i <= numc; i++)<br /> {<br /> str[numc-i] = hash[ans%100];<br /> ans /= 100;<br /> }<br /> str[numc] = '/0';<br /> if (j == n - 1)<br /> {<br /> int k = numc - 1;<br /> for (; k >= 0; k--)<br /> if (str[k] != ' ')<br /> break;<br /> str[k+1] = '/0';<br /> }<br /> printf("%s", str);<br /> }<br /> putchar('/n');<br /> }<br /> return 0;<br />}</p><p>感謝以下：</p><p>http://hi.baidu.com/ve_bird/blog/item/cf585e3f04c82dcb9f3d6203.html</p><p>http://www.cppblog.com/RyanWang/archive/2009/07/19/90525.html</p><p>http://blog.163.com/kakarrot@yeah/blog/static/1201159252010413115321480/</p><p>歐幾裡德演算法(輾轉相除法)是計算兩個數最大公約數的傳統演算法，它無論從理論還是從效率上都是很好的。但是它有一個致命的缺陷，這個缺陷只有在大素數時才會顯現出來。<br />Stein演算法：<br />1、如果A=0，B是最大公約數，演算法結束<br />2、如果B=0，A是最大公約數，演算法結束<br />3、設定A1=A、B1=B和C1=1<br />4、如果An和Bn都是偶數，則An+1=An/2，Bn+1=Bn/2，Cn+1=Cn*2(注意，乘2隻要把整數左移一位即可，除2隻要把整數右移一位即可)<br />5、如果An是偶數，Bn不是偶數，則An+1=An/2，Bn+1=Bn，Cn+1=Cn (很顯然啦，2不是奇數的約數)<br />6、如果Bn是偶數，An不是偶數，則Bn+1=Bn/2，An+1=An，Cn+1=Cn (很顯然啦，2不是奇數的約數)<br />7、如果An和Bn都不是偶數，則An+1=|An-Bn|，Bn+1=min(An,Bn)，Cn+1=Cn<br />8、n加1，轉1 考慮歐幾裡德演算法，最惡劣的情況是，每次迭代a=2b-1,這樣，迭代後，r=b-1。如果a小於2N，這樣大約需要4N次迭代。而考慮Stein演算法，每次迭代後，顯然A(n+1)B(n+1)≤AnBn/2，最大迭代次數也不超過4N次。也就是說，迭代次數幾乎是相等的。但是，需要注意的是，對於大素數，試商法將使每次迭代都更複雜，因此對於大素數Stein將更有優勢。<br />//下面是根據該原理寫的實現函數<br />unsigned MaxDivisor(unsigned a, unsigned b)<br />{<br /> unsigned c = 0;<br /> while(1)<br /> {<br /> // 允出準則<br /> if(a==0)<br /> return b << c;<br /> else if(b == 0)<br /> return a << c;<br /> // 為提高速度，避免用模數判斷奇偶<br /> if(!(a & 1) && !(b & 1)) //a,b 都是偶數<br /> {<br /> a >>= 1; b >>= 1; ++c;<br /> }<br /> else if(!(a & 1) && (b & 1)) //a偶 b奇<br /> {<br /> a >>= 1;<br /> }<br /> else if((a & 1) && !(b & 1)) //a奇 b偶<br /> {<br /> b >>= 1;<br /> }<br /> else if((a & 1) && (b & 1)) //a,b都是奇數<br /> {<br /> unsigned tmp = a>b?b:a; //取較小的一個<br /> a = a>b?a-b:(b-a); //絕對差值<br /> b = tmp;<br /> }<br /> }<br />}</p><p>int Stein(int a,int b) //最優的Stein演算法<br />{<br /> if(a == 0)return b;<br /> if(b == 0)return a;<br /> if(a%2 == 0&&b == 0)return 2*gcd(a/2,b/2);<br /> else if(a%2 == 0)return gcd(a/2,b);<br /> else if(b%2 == 0)return gcd(a,b/2);<br /> else<br /> return gcd(abs(a - b),min(a,b));<br />}<br />