1. 程式人生 > >【數學專題】(一) 中國剩餘定理

【數學專題】(一) 中國剩餘定理

中國剩餘定理 ,求的是模線性方程組的通解,如圖1所示,其中(mi, ai)都是已知量,x是未知量。需要求的就是x,讓它滿足以下n個等式:
圖1 舉個最簡單的例子,如圖2所示,x滿足三個同餘方程,我們可以將三個方程進行一個轉化。 圖2 很明顯,圖2和圖3的方程組等價,圖3表示的是4個未知數,3個方程。由於未知數個數大於方程數,所以這個方程組要麼無解,要麼有無窮多解。
圖3 樸素的做法就是列舉x,然後對每個方程進行試除,令x=(0,1,2,3,...),那麼我們發現當x=26時可以滿足所有三個方程。然而事情並沒有那麼簡單,除數的不同,或者方程數的增多,都可能導致x的值變得很大。 如圖4所示的這種情況,僅僅三個方程,x就已經是10^8的數量級了,所以列舉這種做法顯然是行不通的。 圖4 現在介紹一種O(N)的演算法,其中N為方程組的個數。首先,我們需要定義一些向量,向量a[]代表除數集合,向量m[]代表餘數集合,向量x[]則表示未知數集合。如圖5所示,其中K為需要求的數,我們將n個方程從0開始按順序編號。
圖5 然後,我們合併(0)和(1)兩個等式,得到: 圖6 代換後A,B,C均為常數,X,Y為未知數,然後求A和B的最大公約數G,如果G不能整除C,則方程組無解;否則方程兩邊A,B,C同時除以G,等式不變。接著,就可以通過擴充套件歐幾里得演算法求解Y了,需要注意的是,擴充套件歐幾里得求解時需要保證A和B都大於等於零,所以我們可以做如下處理:
        if(A < 0) {
                // 將A轉化成正數,但還需保證等式性質,所以A,B,C同時變號
        	A = -A, B = -B, C = -C;
                // 將B轉化成正數,C++的取模%可能出現負數,所以需要模完後加上模數再取模
        	B = (B % A + A) % A;
	}
擴充套件歐幾里得演算法求解二元一次方程的解法可以參考這篇文章: 初等數論基礎演算法           Y的通解表示如下(其中Y0已經求出為常數):
圖7 Y就是(1)方程中的x[1],我們將Y代入方程(1),得到等式如下: 圖8 這時候,我們發現兩個等式合併以後的等式,仍然可以表示成 同餘方程 的形式。這說明可以將這個等式和(2)繼續合併,合併完的等式,繼續和(3)合併,直到最後只剩一個等式,令最後這個等式為: K=ax+b,由於x可以取任意整數,所以我們可以保證a和b都大於等於0,並且a >= b,那麼這時候K=b就是方程組的最小非負整數解,K=ax+b就是方程組的通解。 這就是中國剩餘定理的O(N)解法。 最後來看實際演算法流程: 0)列舉 i = 1 to n-1 1) 合併等式(0)和(i), a[0]*x[0] + m[0]= a[i]*x[i] + m[i]; a[0]*x[0] - a[i]*x[i] = m[i] - m[0]; A * X + B * Y = C; 2) 根據擴充套件歐幾里得,求得x[i] = Y = Y0 + A*t; (t為任意整數) 3) 將x[i]代入等式(i),得到: K = a[0]*a[i]*t + (m[i]+a[i]*Y0); 4) 更新等式(0), a[0] = a[0]*a[i] m[0] = (m[i]+a[i]*Y0) % a[0]; 5) n-1次迭代完畢,m[0]就是最小非負整數解。 C++程式碼實現地址:中國剩餘定理O(N)演算法

中國剩餘定理題集解題報告

1 Strange Way to Express Integers 題意:給定N組數(ai, ri)。求一個最小的非負整數K,滿足所有的K mod ai = ri。 題解:中國剩餘定理的模板題。
2 Biorhythms 題意:慢慢讀題吧,讀到最後就是N=3的中國剩餘定理。 題解:直接上模板。
3 X問題 題意:給定N和n組數(ai, ri),求滿足所有的X mod ai = ri,且X<=N的正整數X的個數。 題解:首先套模板求出最小解x,並且求LCM = ai的最小公倍數;然後分情況: 1) x>N 或 無解,輸出0;(我的模板無解返回的是-1,結果因為這個WA了很多次,無語) 2)x=0,輸出N/LCM 3)x>0, 輸出(N-x)/LCM+1 4 Circle 題意:給定長度為n(n <= 20)的約瑟夫環,按順序報數,每次報到k的那個人出列;然後繼續從下一個人報數,每次報到k的那個人出列。現在給定出隊序列,求最小的滿足條件的k。 題解:hash+中國剩餘定理。出隊序列已經給出,那麼可以按順序模擬,每次往環上走一步,如果遇到本次出隊的人,將出隊的人hash掉,記錄上一次到本次的步數。以剩餘的人個數為模數,步數為餘數建立一個方程。N次建立N個方程然後求一次中國剩餘即可。 5 Chinese Remainder Theorem Again 題意:給定陣列Mi和數字a,求一個最小的數K,滿足K mod Mi = (Mi - a)。 題解:看起來赤裸裸的中國剩餘定理,如果直接套模板會WA,因為資料範圍的問題,可能會超int64,所以由於這題的特殊性,有更簡單的解法。還是按照中國剩餘定理的思路,合併等式發現,最後的解其實是所有Mi的最小公倍數再減去a。 6 AndNow,a Remainder from Our Sponsor 題意:按照同餘來加解密的題,有點意思。 題解:模擬 + 中國剩餘定理。
7 HelloKiki 題意:給定N組數(ai, ri)。求一個最小的正整數K,滿足所有的K mod ai = ri。 題解:中國剩餘定理的模板題。注意最小正整數,如果結果為0,需要輸出ai的最小公倍數。
8 A Poor Officer 題意:N(N <= 10000)個傘兵編號遞增排列(即12345...),給定一種重排規則A[N],A[i]表示將第A[i]個傘兵放到第i個位置(1<=i<=N)。然後給定一個目標狀態,問從起始狀態到達目標狀態,最少需要經過多少步。 圖9 題解:置換群 + 中國剩餘定理。首先,對於任意的重排序列,經過有限次重排後必然回到初始狀態,如圖10所示,重排陣列A=[2,1,4,5,7,6,3];目標陣列T=[2,1,7,3,4,6,5]。 圖10 然後對起始狀態和目標狀態分別進行一次重排,就能算出各自的置換群了,如果兩者的置換群不一致,顯然狀態不可達,直接輸出-1即可;如圖11所示,兩者的置換群是一致的,因為(3457)可以通過迴圈左移經過(4573)、(5734)變成(7345),同理,(12)可以經過一步變成(21),所以兩者置換群一致。 圖11 然後,計算起始狀態和目標狀態的置換群中每個迴圈可達需要的最少步數,例如(12)到(21)需要1步,(3457)到(7345)需要3步,(6)到(6)需要0步;然後就可以列出三個方程,即總步數K滿足: 圖12 成功轉換成中國剩餘定理求解。 9 Shuffling 解法同 A Poor Officer 10 Mysterious For 題意:給定n(n <= 10^6),表示迴圈最多執行n次,再給定兩種型別的迴圈,定義如下: 現在總共有m(m <= 10^5)個巢狀,【型別1】的迴圈不大於16個。問所有巢狀的迴圈一共執行多少次。由於次數很大, 需要將答案模上 364875103
題解:Lucas定理 + 中國剩餘定理。 首先,以【型別1】為斷點,如果a[i]是【型別1】的迴圈,那麼它和a[i-1]必然沒關係。圖中XX()的執行次數代表一個【型別1】迴圈嵌套了兩個【型別2】的迴圈。 我們用S(n,t)表示1個【型別1】迴圈和t個【型別2】迴圈的總執行次數。那麼S(n,0)=n, S(n,1)=1+2+...+n。利用數學歸納法,我們總結出S(n,t)的遞推公式:S(n,t)=S(n-1,t)+S(n,t-1)。 然後將n和t記錄到表格中,每個格子的值 = 它上面格子的值 + 它左邊格子的值,如圖所示:

觀察可得,它是一個斜著的楊輝三角,並且可以轉化成組合數求解:S(n,t)=C(n+t,t+1)。於是轉變成了求組合數C(a,b) mod 364875103。 C(a,b) mod 364875103的求法採用Lucas定理 + 中國剩餘定理。 由於364875103 = 97 * 3761599,為兩個素數的乘積,我們可以先求出x = C(a,b) mod 97,和y = C(a, b) mod 3761599,那麼令K = C(a,b) mod 364875103,K必然滿足: 如果x和y都是已知的,那麼就可以通過中國剩餘定理來求K了。 合數C(a,b)%p(其中p為素數)可以採用Lucas定理求解。Lucas定理的遞迴形式如下:
Lucas(n, m, p) {
	if(m == 0) return 1;
	return Lucas(n/p, m/p, p) * Comb(n%p, m%p, p) % p;
}
其中Comb(n,m,p)表示求C(n,m)%p (n<p && m<p && p為素數) 那麼現在就是要求n! mod p的值,其中n的範圍<=1100000,直接預處理存到陣列中,O(1)獲取即可。n!^(p-2)mod p可以採用二分快速冪取模。
11 The Lastest Math Theory Problem 最後推薦一道我出的題^_^ 題意:給定n(n <= 6)和陣列a[n](a[i] < 10000),要求求一個最小的N,滿足n個等式(其中x[i]為未知量): 由於N可能很大,所以只需要輸出N mod 19880502即可;如果不存在這樣的N,則輸出-1。 題解:中國剩餘定理 + 二分快速冪。任何數都可以表示成素數的乘積,又由於N滿足以上六個等式,所以N的素因子一定是a[i]的素因子,假設N的表示如下:

任何一個p[i]必定小於10000,那麼我們只需要列舉所有小於10000的素因子,然後想辦法求出指數e[i]即可。對於每個小於10000的素數p[i],列舉所有a[j], 檢查a[j]這個數中有多少p的素因子,記為cnt[j],則e[i] = k*a[j] + cnt[j](k>=0);那麼對於每個e[i]利用中國剩餘定理求n個同餘方程,得到的e[i]利用二分快速冪求p[i]^e[i] mod 19880502即可。最後所有素數分量連乘再取模就是最後的答案了。 這題需要注意一點,求中國剩餘定理的時候,由於值的範圍在int64,所以乘法可能溢位。所以如果兩個數相乘需要採用二進位制法。這裡提供一種非遞迴版本的a*b mod c,全程只有加法和取模。

LL Product_Mod(LL a, LL b, LL c) {
	LL sum = 0;
	while(b) {
		if(b & 1) sum = (sum + a) % c;
		a = (a + a) % c;
		b >>= 1;
	}
	return (sum + c) % c;
}