擴充套件歐幾里德演算法求解線性同餘方程
阿新 • • 發佈：2019-01-03
歐幾里德演算法
歐幾里德演算法又稱輾轉相除法，用於計算兩個整數a,b的最大公約數。其計算原理依賴於下面的定理：
定理：gcd(a,b) = gcd(b,a mod b)
證明：a可以表示成a = kb + r，則r = a mod b
假設d是a,b的一個公約數，則有
d|a, d|b，而r = a - kb，因此d|r
因此d是(b,a mod b)的公約數
假設d是(b,a mod b)的公約數，則
d | b , d |r，但是a = kb +r
因此d也是(a,b)的公約數
因此(a,b)和(b,a mod b)的公約數是一樣的，其最大公約數也必然相等，得證
歐幾里德演算法就是根據這個原理來做的，其演算法用C++語言描述為：
- int GCD(int a,int b)
- {
- while (b!=0) { int k=b; b=a%b; a=k; }
- return a;
- }
（我做一下補充：C++裡面有一個庫函式可以直接求最大公約數，即__gcd(int, int)，前面是兩個連續的下劃線，這樣就不必要每次都自己寫了。。。）
擴充套件歐幾里德演算法
擴充套件歐幾里德演算法是用來在已知a, b求解一組p，q使得p * a+q * b = Gcd(a, b) (解一定存在，根據數論中的相關定理)。
演算法描述為：
- intextended_gcd(int a,int b,int &x,int &y)
- {
- int ans,t;
- if (b==0) { x=1; y=0; return a; }
- else { ans=extended_gcd(b,a%b,x,y); t=x; x=y; y=t-(a/b)*y;}
- return ans;
- }
把這個實現和Gcd的遞迴實現相比，發現多了下面的x,y賦值過程，這就是擴充套件歐幾里德演算法的精髓。
可以這樣思考:
對於a' = b, b' = a % b而言，我們求得x, y使得a'x + b'y = Gcd(a', b')
由於b' = a % b = a - a / b * b (注：這裡的/是程式設計語言中的除法)
那麼可以得到:
a'x + b'y = Gcd(a', b') ===>
bx + (a - a / b * b)y = Gcd(a', b') = Gcd(a, b) ===>
ay +b(x - a / b*y) = Gcd(a, b)
因此對於a和b而言，他們的相對應的p，q分別是y和(x-a/b*y)。
線性同餘方程
對於方程 a*x+b*y=n；有整數解得充分必要條件是（n %（a，b）==0），這個定理這裡就不證明了，數論書上都有。
所以方程 a*x+b*y=n；我們可以先用擴充套件歐幾里德演算法求出一組x0，y0。也就是a*x0+b*y0=（a，b）；然後兩邊同時除以（a，b），再乘以n。這樣就得到了方程a*x0*n/（a，b）+b*y0*n/（a，b）=n；我們也就找到了方程的一個解。
還有一個定理：若（a，b）=1，且x0，y0為a*x+b*y=n的一組解，則該方程的任一解可表示為：x=x0+b*t，y=y0-a*t；且對任一整數t，皆成立。（這個證明比較簡單，就不寫了）
這樣我們就可以求出方程的所有解了，但實際問題中，我們往往被要求去求最小整數解，所以我們就可以將一個特解x，t=b/（a，b），x=(x%t+t）%t；就可以了。
方程組的情形（中國剩餘定理）
對於同餘方程組：
x=a1 （mod m1）； 1
x=a2 （mod m2）； 2
方程組有一個小於m（m1，m2的最小公倍數）的非負整數解的充分必要條件是（a1-a2）%（m1，m2）==0 ，同樣利用擴充套件歐幾里德演算法。
兩式聯立：a1+m1*y=a2+m2*z。
則：a1-a2=m2*z-m1*y; 這樣就可以瞭解出z和y，則：x=a2+m2*z；
現在我們將其推廣到一般情形：（設m1,m2,···,mk兩兩互素）
x=a1（mod m1）;
x=a2（mod m2）;
···
x=ak（mod mk）;其在M=m1*m2*···*mk;中有唯一整數解。
記Mi=M/mi;因為（Mi,mi）=1,故有兩整數pi,qi滿足Mi*pi+mi*qi=1，如果記ei=Mi*pi;那麼：ei=0 （mod mj）,j!=i; ei=1（mod mj）,j=i;
很明顯，e1*a1+e2*a2+···+ek*ak就是方程的一個解，加減M倍後就可以得到最小非負整數解了。
如果m1,m2,···,mk不互素，那隻能兩個兩個求了。
x=a1 （mod m1）；
x=a2 （mod m2）；
解完後，a=x； m=m1和m2的最小公倍數。即可。
暫時就寫到這裡，多元的高次的我還沒看，這裡推薦幾道題目：pku2115；pku2891；pku1061；pku1006；pku2142；強烈推薦sgu106；題目的結題報告我也會貼在部落格裡。