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

謝照林

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

0 引 言

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

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

求x=?

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

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

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

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

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

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

1) 經(jīng)典算法(中國余數(shù)定理)問題的提出:

給定同余式組x≡a1(modm1),x≡a2(modm2),…,x≡an(modmn),其中,m1,m2,m3,…,mn都是互質(zhì)的(互為質(zhì)數(shù),它們的最大公約數(shù)為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都是互質(zhì)的,求解x。以下稱n為該組的規(guī)模。杠桿算法的解題步驟是先解兩個同余式,然后把它的結(jié)果與第3個同余式求解,依此類推直至n個同余式都被用到。

2.1 兩個同余式的求解

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

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

q2m2+a2-a1=q1m1。

(1)

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

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

(2)

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

模數(shù)運算具有如下性質(zhì)[15-16]:

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

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

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

Δ= (q2m2+a2-a1) modm1=

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

(3)

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

Δ0=(a2-a1) modm1。

(4)

當(dāng)q2=1時,Δq2=1比Δq2=0多一項m2modm1,這是q2從0變?yōu)?時Δ的增量,由于q2的變化量是1,Δ的增量也就是函數(shù)Δ(q2)在起點的差分。由于Δ(q2) 是線性函數(shù),它也是在q2全域各點的差分。可以用符號δΔ把它定義為

δΔ=m2modm1,

(5)

由此,式(3)又可寫為

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

(6)

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

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

(7)

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

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

2) 式(7)等價于

q2*δΔ+Δ0=km1,

(8)

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

于是,可以得到

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試探,已經(jīng)得到x=q2m2+a2。按照中國余數(shù)定理,這個同余數(shù)是對于模數(shù)乘積m1m2而言的,即

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

如果使用新變量

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

(11)

把新加入的第3個同余式的下標(biāo)由3改為1,即

a1=a3,m1=m3,

(12)

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

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

(13)

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

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

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

δα=m1modδΔ。

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

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

(14)

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

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

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

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

(15)

即以p來求k。

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

同理,定義

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

(16)

可得

δβ=δΔmodδα,

定義

Δ000=Δ00modδα。

設(shè)定求解目標(biāo):β=0。由此得到第3級杠桿

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

(17)

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

同理,可得第4級杠桿

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

(18)

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

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

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

表1 各級杠桿的表達(dá)形式

由表1 可以看出,各級杠桿表達(dá)式的結(jié)構(gòu)是相同的,可以用一個通用的形式來概括,如表2 所示。

表2 杠桿迭代的通用形式

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

2) 迭代的實例:在計算機上進(jìn)行實例計算所得結(jié)果為

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)

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

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

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

下面討論δεi= 1 是否一定會出現(xiàn)。

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

2.5 杠桿算法的步驟

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

1) 對于規(guī)模在2以上的同余式組按照模數(shù)值從大到小進(jìn)行排序。這個排序?qū)τ诿恳坏李}目只做一次,目的是使各次外層迭代的運算時間盡可能均衡,從而減少總的運算時間。氣泡排序就是常見的有效算法。當(dāng)然,是否做模數(shù)排序?qū)τ嬎憬Y(jié)果沒有影響。

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

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

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

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

3 經(jīng)典算法和杠桿算法的比較

3.1 可處理的數(shù)值范圍的比較

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

表3 杠桿算法可處理的最大數(shù)值是經(jīng)典算法可處理最大數(shù)值的倍數(shù)

3.2 運算時間的比較

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

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

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

以上每一例杠桿算法運算時間都已經(jīng)包括了模數(shù)排序的運算時間。

4 結(jié)束語

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

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

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