解一元線性同余式組的一個新算法

謝照林

(美國布魯克斯自動化公司,加州 佛利蒙市 94538)

0 引 言

一元線性同余式組是數論的一個組成部分。早在公元5世紀,中國的數學著作《孫子算經》就提出了一個“物不知數”問題:“有物不知其數,三三數之剩二,五五數之剩三,七七數之剩二。問物幾何?”即,一個整數除以三余二,除以五余三,除以七余二,求這個整數。這個問題用現代數學語言可以表達為

x=2(mod 3),x=3(mod 5),x=2(mod 7),

求x=?

孫子的答案為x=2*70+3*21+2*15-2*105=23,用現代的同余式表達為x≡23(mod 105)。

1247年,宋朝數學家秦九韶在《數書九章》的《大衍類》中對“物不知數” 問題作出了完整系統的解答[1-8]。在國外,從6世紀到19世紀諸多印度和歐洲的科學家都提出了對同余式組的算法[9-10]。1801年,德國科學家高斯(Carl Friedrich Gauss) 提出了一個和秦九韶一致的算法,這是歐洲最早的完整系統算法[11]。1852年,英國傳教士偉烈亞力(Alexander Wylie)把《數書九章》翻譯成英文[12],它很快又被翻譯成德文和法文,并在西方廣泛傳播。19 世紀末20世紀初,西方國家把秦九韶和高斯對一元線性同余式組的算法命名為中國余數定理。中國余數定理具有極高的理論和實用價值,已成為解一元線性同余式組算法的經典,本文稱它為經典算法,并將在第1節對其進行簡要介紹。

近年來,國內外學者對一元線性同余式組問題的研究還在持續不斷[13-19],研究的重點是它在各個領域的應用。除此之外,在互聯網上對于擴展的同余式或同余式組的各種不同解法的討論也很活躍,但是有關中國余數定理提出的基本的一元線性同余式組問題則無論在互聯網上或是出版文獻上都極少見到有人提出新的解法。偶爾見到的有窮舉法、解不定方程法、化為相同除數法,等等。這些解法擴大了人們的思路,也各有其適用的場合,不足的是求解過程往往因題而異,依賴人工判斷,也未見到這些解法各自的可用于計算機計算的通用算法。

經典算法的計算步驟清晰且通用,在它形成的年代計算機尚未出現,但是它卻成為一個現成的可用于計算機計算的通用算法。本文提出解一元線性同余式組新算法的初衷,是想探討能否在經典算法之外找出另一種更高效的可用于計算機計算的通用算法。在現在這樣一個需要快速處理大量大數值數據的時代,多一種算法供使用者選擇總是有益無害的。

本文提出的新算法的思路和經典算法完全不同,在文中被稱為杠桿算法。計算機對比計算的結果表明,與經典算法相比,杠桿算法可以處理數值更大的數據以及具有更快的運算速度。盡管如此,杠桿算法也僅僅是一個算法,而不是一個理論,完全無法與經典算法相比擬,可以把它看作一個不同的解題計算工具。

1 一元線性同余式組的經典算法(中國余數定理)

1) 經典算法(中國余數定理)問題的提出:

給定同余式組x≡a1(modm1),x≡a2(modm2),…,x≡an(modmn),其中,m1,m2,m3,…,mn都是互質的(互為質數,它們的最大公約數為1),求解x。

2) 解題步驟可歸納如下:

① 計算乘積M=m1m2m3…mn;

② 對于i=1,2,…,n,依次計算ui=M/mi=m1m2…mi-1mi+1…mn;

2 一元線性同余式組的新算法(杠桿算法)

給定同余式組x≡a1(modm1),x≡a2(modm2),…,x≡an(modmn),其中,m1,m2,m3,…,mn都是互質的,求解x。以下稱n為該組的規模。杠桿算法的解題步驟是先解兩個同余式,然后把它的結果與第3個同余式求解,依此類推直至n個同余式都被用到。

2.1 兩個同余式的求解

x≡a1(modm1) 可以寫為x=q1m1+a1。其中,q1是x除以m1的商,a1為余數,m1為模數。x≡a2(modm2) 可以寫為x=q2m2+a2。其中,q2是x除以m2的商,a2為余數,m2為模數。

約定兩個同余式中模數大的同余式稱為x≡a2(modm2),另一個稱為x≡a1(modm1),兩式相減可得

q2m2+a2-a1=q1m1。

(1)

式(1)表明,(q2m2+a2-a1) 必須是m1的整數倍,可表達為

Δ=(q2m2+a2-a1) modm1=0。

(2)

求解目標1是找出滿足Δ=0 的q2值,則x=q2m2+a2就是求得的解。

模數運算具有如下性質[15-16]:

(A+B) modC=(AmodC+BmodC) modC,

(A*B) modC=(AmodC*BmodC) modC,

可以看出,只要括號外有modC,那么對于括號內的各項,modC可以任意增減。因此,有

Δ= (q2m2+a2-a1) modm1=

(q2(m2modm1)+(a2-a1) modm1)modm1。

(3)

當q2=0時,式(3)成為Δq2=0=(a2-a1)modm1,這是Δ的初值,可以用符號Δ0把它定義為

Δ0=(a2-a1) modm1。

(4)

當q2=1時,Δq2=1比Δq2=0多一項m2modm1,這是q2從0變為1時Δ的增量,由于q2的變化量是1,Δ的增量也就是函數Δ(q2)在起點的差分。由于Δ(q2) 是線性函數,它也是在q2全域各點的差分?？梢杂梅枽摩ぐ阉x為

δΔ=m2modm1,

(5)

由此,式(3)又可寫為

Δ=(q2*δΔ+Δ0) modm1。

(6)

參數a1,m1,a2,m2是一組同余式的固有參數。Δ0和δΔ包含了這4個參數,所以Δ0和δΔ就成為這一組同余式的特征參數,是杠桿算法的出發點。求解目標2是找出q2值使式(7)成立,

Δ=(q2*δΔ+Δ0) modm1= 0。

(7)

得到滿足條件的q2值的方法有兩種:

1) 逐一試探q2=1,2,3,…,直到式(7)成立。

2) 式(7)等價于

q2*δΔ+Δ0=km1,

(8)

式中:k是一個最小可能的正整數。

于是,可以得到

q2= (km1-Δ0)/δΔ。

(9)

逐一試探k=1,2,3,…,直到

α=(km1-Δ0) modδΔ=0,

(10)

把找到的k值代入式(9),就可以求出q2。

由于δΔ=m2modm1,所以δΔ1。所以,由任何k值所計算出的q2都比k大。以k求q2,就是以小值求大值。這就是下面要討論的杠桿。

2.2 三個以上同余式的求解

無論用q2試探還是k試探,已經得到x=q2m2+a2。按照中國余數定理,這個同余數是對于模數乘積m1m2而言的,即

x=(q2m2+a2) (mod (m1m2))。

如果使用新變量

a2=(q2m2+a2),m2=(m1m2),

(11)

把新加入的第3個同余式的下標由3改為1,即

a1=a3,m1=m3,

(12)

那么,可以得到以下兩個新的同余式

x≡a1(modm1),x≡a2(modm2)。

(13)

在解出3個同余式(13)以后,重復上述式(11)～式(13)的方法,總共進行n-1次就可以完成n個同余式的求解。

2.3 杠桿算法(減少試探次數)

k試探比q2試探高效,但隨著模數數值的增大,k試探也越來越耗時,必須找到加速k試探的途徑?；仡櫴?10),如同處理由Δ表達式得出δΔ一樣,可以得到

δα=m1modδΔ。

定義Δ00=Δ0modδΔ,則式(10) 可以寫為

α=(k*δα-Δ00) modδΔ=0,

(14)

也就是(k*δα-Δ00) 必須是δΔ的整數倍,這里需要的是最小的整數倍。為此,

k*δα-Δ00=p*δΔ,

式中:p是一個最小可能的正整數。故

k=(p*δΔ+Δ00)/δα,

(15)

即以p來求k。

δΔ/δα>1,所以任何p值都能產生一個比它大的k值。p和k的關系又是一個杠桿,這是第2級杠桿。

同理,定義

β=(p*δΔ+Δ00) modδα,

(16)

可得

δβ=δΔmodδα,

定義

Δ000=Δ00modδα。

設定求解目標:β=0。由此得到第3級杠桿

p=(s*δα-Δ000)/δβ,

(17)

式中:s是一個最小可能的正整數。這是以s求p。

同理,可得第4級杠桿

s=(t*δβ+Δ0000)/δγ,

(18)

式中:t是一個最小可能的正整數。這是以t求s。

2.4 杠桿迭代(徹底免除試探)

1) 迭代的形式:匯總以上用過的各個表達式,如表1 所示。

表1 各級杠桿的表達形式

由表1 可以看出,各級杠桿表達式的結構是相同的,可以用一個通用的形式來概括,如表2 所示。

表2 杠桿迭代的通用形式

表2 中Δi的下標i代表迭代次序[i]。在迭代[1]中Δi是Δ0,在迭代[2]中Δi是Δ00,依此類推; 符號Δi等效于Δ的下標為i個0,每次迭代中它都是該次迭代運算的初值。表2 中各表達式的形式更便于進行迭代運算,稱這樣的迭代為杠桿迭代。其中,第1步是正向迭代:按照[0],[1],[2],[3],…的順序自上而下計算第1列和第2列; 第2步是反向迭代:按照相反次序自下而上計算第3列,直到算出k1。對照表1 和表2 可知,k1就是q2。

2) 迭代的實例:在計算機上進行實例計算所得結果為

Solve(4 801,9 739),(5 107,10 009):

[1]:D0=306,dEps0=10 009,dEps1=9 739,dEps2=270,bSignOfD=1

[2]:D0=36,dEps0=9 739,dEps1=270,dEps2=19,bSignOfD=0

[3]:D0=17,dEps0=270,dEps1=19,dEps2=4,bSignOfD=1

[4]:D0=1,dEps0=19,dEps1=4,dEps2=3,bSignOfD=0

[5]:D0=1,dEps0=4,dEps1=3,dEps2=1,bSignOfD=1

dEps2=1,sosetk[6]=1

[5R]:k[6]=1,k=2

[4R]:k[5]=2,k=3

[3R]:k[4]=3,k=10

[2R]:k[3]=10,k=144

[1R]:k[2]=144,k=5 193

Finally,k=144,q2=5 193,x=51 981 844

Thecongruenceisx=51 981 844 (mod97 477 651)

其中,方括號里的數字代表迭代的次序; 字母R代表反向;D0在各次迭代中分別代表Δ0,Δ00,Δ000,Δ0000,…;dEps0,dEp1s,dEps2 在迭代[2]中代表δε0,δε1和 δε2,在迭代[3]中代表δε1,δε2和 δε3,依此類推; 布爾變量bSignOfD代表計算ki值時Δi應有的正負號; 1代表應該用負號,0代表應該用正號。從該實例可以看出,5次迭代計算就得到 k=144,q2=5 193,比k試探,尤其是比q2試探,效率提高了很多倍。

3) 正向迭代結束的條件:在正向的每一次迭代中都有對δεi=δεi-2modδεi-1的計算,所以,當δεi= 1時,正向迭代就完成了。理由如下:

表2 的每一行都有ki=(ki+1*δεi-1±Δi)/δεi,當δεi= 1時,取任何整數ki+1都可以得到ki為整數。這里要取最小的正整數ki+1=1。所以,正向迭代結束的條件就是δεi=1。此時,取ki+1=1。

下面討論δεi= 1 是否一定會出現。

表2 中正向迭代[2]是從δε2= δε0modδε1開始的,也就是從δε2= m1mod(m2modm1) 開始的。其后δε3= δε1modδε2就是δε3= (m2modm1)mod(m1mod(m2modm1)),依此類推。這就是m2和m1的輾轉相除,即求m2和 m1的最大公約數(GCD)的過程。而m2和 m1是互質的,即它們的GCD是1。因此,δεi= 1一定會出現,至于δεi= 1在第幾次迭代中出現,取決于題目數據m2和 m1的值。

2.5 杠桿算法的步驟

杠桿算法用到兩層迭代。外層迭代是先解兩個同余式,然后其結果又和第3個同余式一起按照解兩個同余式的過程求解,其結果又和第4個同余式一起求解,依此類推,總共需要n-1次,n是同余式組的規模。內層迭代就是杠桿迭代,它是在每一次外層迭代中解兩個同余式的過程中都要使用的。設每一次外層迭代的已知數據是 (a1,m1),(a2,m2),杠桿算法的解題步驟如下:

1) 對于規模在2以上的同余式組按照模數值從大到小進行排序。這個排序對于每一道題目只做一次,目的是使各次外層迭代的運算時間盡可能均衡,從而減少總的運算時間。氣泡排序就是常見的有效算法。當然,是否做模數排序對計算結果沒有影響。

2) 計算常數δε0=m1,Δ0=(a2-a1)modm1,δε1= m2modm1。

3) 按照表2逐次正向迭代計算Δi和δεi。

4) 當δεi= 1,令ki+1= 1,開始反向迭代,直到算出k1。這個k1是q2。

5) 計算a2=q2m2+ a2,m2=m1m2。把下一個同余式命名為(a1,m1),進入下一次外層迭代的步驟2)。

3 經典算法和杠桿算法的比較

3.1 可處理的數值范圍的比較

一個64位(二進制位數)計算機能接受的最大正整數是Z=9 223 372 036 854 775 807,即十進制19位數。

表3 杠桿算法可處理的最大數值是經典算法可處理最大數值的倍數

3.2 運算時間的比較

解一道一元線性同余式組問題的運算時間是微秒級的,在通常條件下無法直接測得。在以下的實例中每道題都重復計算了4 000萬次,然后算出運算一次的平均時間,計算結果如圖1 所示。

(a) 運算時間隨最大模數值的變化

圖1(a) 中每一個算例都是一個規模為3的同余式組。最大模數值是指3個同余式中最大的模數值。由圖1(a) 可知,運算時間隨最大模數值的增大而增大。但是,對于最大模數值的跨度101～1 000 003 而言,增長率可視為0。杠桿算法運算時間的平均值是0.87 μs,比經典算法的1.33 μs減少35%。由圖1(b)可知,杠桿算法運算時間的平均值是1.53 μs,比經典算法的2.08 μs減少27%。

以上每一例杠桿算法運算時間都已經包括了模數排序的運算時間。

4 結束語

1) 單純從求一元線性同余式組解題答案的角度來看,杠桿算法和經典算法是兩個并行的算法,二者殊途同歸。二者計算任何題目所得答案完全相同。但是,從數學學科的角度來看,經典算法具有嚴謹的數學思維,是數學殿堂中的一件瑰寶,而杠桿算法只作為計算工具。

2) 與經典算法相比,在相同的計算機條件下,杠桿算法可以處理數值更大的數據以及具有更短的運算時間。

4) 今后可以嘗試在解擴展的同余式和同余式組問題方面進行高效的計算機算法的探討。