謝照林

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

0 引 言

一元線性同余式組是數(shù)論的一個(gè)組成部分。早在公元5世紀(jì),中國(guó)的數(shù)學(xué)著作《孫子算經(jīng)》就提出了一個(gè)“物不知數(shù)”問(wèn)題:“有物不知其數(shù),三三數(shù)之剩二,五五數(shù)之剩三,七七數(shù)之剩二。問(wèn)物幾何?”即,一個(gè)整數(shù)除以三余二,除以五余三,除以七余二,求這個(gè)整數(shù)。這個(gè)問(wèn)題用現(xiàn)代數(shù)學(xué)語(yǔ)言可以表達(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ù)書九章》的《大衍類》中對(duì)“物不知數(shù)” 問(wèn)題作出了完整系統(tǒng)的解答[1-8]。在國(guó)外,從6世紀(jì)到19世紀(jì)諸多印度和歐洲的科學(xué)家都提出了對(duì)同余式組的算法[9-10]。1801年,德國(guó)科學(xué)家高斯(Carl Friedrich Gauss) 提出了一個(gè)和秦九韶一致的算法,這是歐洲最早的完整系統(tǒng)算法[11]。1852年,英國(guó)傳教士偉烈亞力(Alexander Wylie)把《數(shù)書九章》翻譯成英文[12],它很快又被翻譯成德文和法文,并在西方廣泛傳播。19 世紀(jì)末20世紀(jì)初,西方國(guó)家把秦九韶和高斯對(duì)一元線性同余式組的算法命名為中國(guó)余數(shù)定理。中國(guó)余數(shù)定理具有極高的理論和實(shí)用價(jià)值,已成為解一元線性同余式組算法的經(jīng)典,本文稱它為經(jīng)典算法,并將在第1節(jié)對(duì)其進(jìn)行簡(jiǎn)要介紹。

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

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

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

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

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

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

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

① 計(jì)算乘積M=m1m2m3…mn;

② 對(duì)于i=1,2,…,n,依次計(jì)算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ī)模。杠桿算法的解題步驟是先解兩個(gè)同余式,然后把它的結(jié)果與第3個(gè)同余式求解,依此類推直至n個(gè)同余式都被用到。

2.1 兩個(gè)同余式的求解

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ù)。

約定兩個(gè)同余式中模數(shù)大的同余式稱為x≡a2(modm2),另一個(gè)稱為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ù)運(yùn)算具有如下性質(zhì)[15-16]:

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

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

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

Δ= (q2m2+a2-a1) modm1=

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

(3)

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

Δ0=(a2-a1) modm1。

(4)

當(dāng)q2=1時(shí),Δq2=1比Δq2=0多一項(xiàng)m2modm1,這是q2從0變?yōu)?時(shí)Δ的增量,由于q2的變化量是1,Δ的增量也就是函數(shù)Δ(q2)在起點(diǎn)的差分。由于Δ(q2) 是線性函數(shù),它也是在q2全域各點(diǎn)的差分?？梢杂梅?hào)δΔ把它定義為

δΔ=m2modm1,

(5)

由此,式(3)又可寫為

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

(6)

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

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

(7)

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

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

2) 式(7)等價(jià)于

q2*δΔ+Δ0=km1,

(8)

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

于是,可以得到

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

(9)

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

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

(10)

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

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

2.2 三個(gè)以上同余式的求解

無(wú)論用q2試探還是k試探,已經(jīng)得到x=q2m2+a2。按照中國(guó)余數(shù)定理,這個(gè)同余數(shù)是對(duì)于模數(shù)乘積m1m2而言的,即

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

如果使用新變量

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

(11)

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

a1=a3,m1=m3,

(12)

那么,可以得到以下兩個(gè)新的同余式

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

(13)

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

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

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

δα=m1modδΔ。

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

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

(14)

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

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

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

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

(15)

即以p來(lái)求k。

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

同理,定義

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

(16)

可得

δβ=δΔmodδα,

定義

Δ000=Δ00modδα。

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

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

(17)

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

同理,可得第4級(jí)杠桿

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

(18)

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

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

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

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

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

表2 杠桿迭代的通用形式

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

2) 迭代的實(shí)例:在計(jì)算機(jī)上進(jìn)行實(shí)例計(jì)算所得結(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)

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

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

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

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

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

2.5 杠桿算法的步驟

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

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

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

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

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

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

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

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

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

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

3.2 運(yùn)算時(shí)間的比較

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

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

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

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

4 結(jié)束語(yǔ)

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

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

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