一個(gè)費(fèi)馬數(shù)分解算法的剖析與優(yōu)化

王玨

（廣東移遠(yuǎn)通訊技術(shù)有限公司，佛山528042）

分析一個(gè)費(fèi)馬數(shù)分解算法中的冗余步驟，給出相應(yīng)的優(yōu)化結(jié)果。針對(duì)相關(guān)文獻(xiàn)述及大費(fèi)馬數(shù)表示困難的問題，給出利用GMP大數(shù)運(yùn)算庫表示大費(fèi)馬數(shù)的一種方法。基于數(shù)學(xué)軟件Maple，進(jìn)行解析費(fèi)馬數(shù)小因數(shù)的試驗(yàn)。試驗(yàn)表明，優(yōu)化后的算法可提高計(jì)算效率。

整數(shù)分解；費(fèi)馬數(shù)；算法；密碼學(xué)

0 引言

算法，始終是程序設(shè)計(jì)人員關(guān)注的話題。作為一個(gè)程序員，工作之余了解一下國內(nèi)外與算法相關(guān)新聞、報(bào)道和事件，不失為另一類娛樂。最近看到Journal of Software（JSW）期刊上新出的一篇有關(guān)分解費(fèi)馬數(shù)的論文[1]，覺得很有意思，便仔細(xì)研讀了一下。該文介紹的是一個(gè)分解費(fèi)馬數(shù)的新方法，在理論上分解任何Fn（n>100001）的任何費(fèi)馬數(shù)。眾所周知，費(fèi)馬數(shù)是數(shù)論里面一類受人關(guān)注的數(shù)，國內(nèi)學(xué)術(shù)期刊也多有研究報(bào)道。例如，文獻(xiàn)[2-3]較為通俗地介紹了費(fèi)馬數(shù)的歷史與發(fā)展，文獻(xiàn)[4]較為系統(tǒng)地研究了費(fèi)馬數(shù)。也有學(xué)者對(duì)費(fèi)馬數(shù)的性質(zhì)進(jìn)行挖掘分析，如文獻(xiàn)[5-12]分別就費(fèi)馬數(shù)的因數(shù)性質(zhì)進(jìn)行了探討。在分解計(jì)算方面，文獻(xiàn)[13]介紹了利用橢圓曲線分解F7和F8的過程，文獻(xiàn)[14]則介紹了國外分解計(jì)算的情況。國外有專門的網(wǎng)站介紹分解費(fèi)馬數(shù)的進(jìn)展，如文獻(xiàn)[15]。從筆者不是十分專業(yè)的角度來看，除了目前看到的文獻(xiàn)[1]，國內(nèi)學(xué)者對(duì)于研究費(fèi)馬數(shù)分解算法的不多。

鑒于此，筆者對(duì)文獻(xiàn)[1]給出的算法進(jìn)行了認(rèn)真分析，發(fā)現(xiàn)這個(gè)算法雖然很巧妙，卻也存在可以改進(jìn)的地方。本文介紹筆者的分析和改進(jìn)結(jié)果。

1 問題的提出及解決方案

1.1 問題的提出

仔細(xì)研讀并按照文獻(xiàn)[1]給出的算法在Maple軟件編程時(shí)，會(huì)發(fā)現(xiàn)該文存在幾處可以優(yōu)化的地方。首先是該文在結(jié)論和展望部分提出的二個(gè)有待進(jìn)一步解決的問題：一是大數(shù)運(yùn)算工具的問題，二是確定算法中的參數(shù)c實(shí)現(xiàn)加速計(jì)算的問題。該文還列舉了m>23時(shí)，GMP大數(shù)庫和Maple都出現(xiàn)溢出的情況，可以推斷該文應(yīng)該是在利用目前常用的GMP大數(shù)運(yùn)算庫進(jìn)行了C/C++編程測(cè)試和利用Maple數(shù)學(xué)軟件測(cè)試時(shí)遇到了計(jì)算溢出。該文同時(shí)也指出NTL大數(shù)運(yùn)算庫和MIRACL大數(shù)運(yùn)算庫都不能完成運(yùn)算。于是該文作者顯得無可奈何地說，這些工具的開發(fā)已經(jīng)超出了他的能力范圍。誠然，開發(fā)大數(shù)運(yùn)算的工具，的確超越數(shù)學(xué)家們的工作范圍。

一個(gè)值得注意的事情是，文獻(xiàn)[1]在JSW上正式發(fā)表以前有個(gè)Preprint版本，公開在文獻(xiàn)[16]。在這個(gè)Preprint版本中，作者哀嘆說所遇到的困難是，在m>2w-1時(shí)無法用GMP大數(shù)庫來表達(dá)Fn導(dǎo)致不能計(jì)算大費(fèi)馬數(shù)的分解，這里w是計(jì)算機(jī)的字長。而在JSW正式發(fā)表的版本里，作者已經(jīng)沒有這個(gè)哀嘆了。說明該作者已經(jīng)利用GMP大數(shù)庫解決了大費(fèi)馬數(shù)的表示問題了，反而是在mpz_gcd這個(gè)庫函數(shù)上遇到了新的問題。筆者經(jīng)過分析，也找到了一個(gè)表示費(fèi)馬數(shù)的方法。鑒于文獻(xiàn)[1]沒有列出它采用的方法，本文在后面小節(jié)給出筆者的方法。

對(duì)于第二個(gè)問題，即確定算法中的參數(shù)c實(shí)現(xiàn)加速計(jì)算的問題，該文作者沒有進(jìn)行多的說明，僅僅呼吁更多年輕人來參與解決所述問題。分析文[1]對(duì)c的幾處說明，發(fā)現(xiàn)這是一個(gè)減少搜索次數(shù)的參數(shù)，這是因?yàn)樗鼭M足b=■ ■log2N-χ，且在幾個(gè)搜索過程中都出現(xiàn)了for ifromb downto3 do的語句。

鑒于文獻(xiàn)[1]沒有詳細(xì)介紹這個(gè)參數(shù)的選擇而是把它當(dāng)成一個(gè)問題提出，筆者通過研究并經(jīng)Maple編程測(cè)試，發(fā)現(xiàn)了一些規(guī)律，后文將給予介紹。

除了前述兩個(gè)問題之外，文獻(xiàn)[1]還存在可以挖掘的地方，那就是具體分解費(fèi)馬數(shù)的Procedure 4。筆者經(jīng)研究發(fā)現(xiàn)，這個(gè)過程存在兩個(gè)可以優(yōu)化的地方。一個(gè)是該過程中的參數(shù)a的上限，二是該過程中的基本計(jì)算。詳見后文介紹。

1.2 解決方案

本小節(jié)給出大費(fèi)馬數(shù)的表示方案、文獻(xiàn)[1]中Procedure4的優(yōu)化方法以及參數(shù)c的選擇方案。

（1）大費(fèi)馬數(shù)的表示

針對(duì)前小節(jié)的第一個(gè)問題，筆者經(jīng)仔細(xì)分析NTL大數(shù)運(yùn)算庫和MIRACL大數(shù)運(yùn)算庫的文檔發(fā)現(xiàn)，這兩個(gè)庫的確無法表示大費(fèi)馬數(shù)。經(jīng)分析GMP大數(shù)庫的文檔發(fā)現(xiàn)，利用mpz_ui_pow_ui和mpz_pow_ui兩個(gè)庫函數(shù)，可以解決大費(fèi)馬數(shù)的表示問題。mpz_ui_pow_ui和mpz_pow_ui兩個(gè)庫函數(shù)的原型與功能分別如表1所示。

表1函數(shù)原型與功能

利用這兩個(gè)函數(shù)，可以實(shí)現(xiàn)大費(fèi)馬數(shù)的表示。具體思路如下：

對(duì)于任意無符號(hào)長整型的數(shù)m，利用Euclid除法將m按照計(jì)算機(jī)的字長w表示成

式中s≥0，0≤r＜w均為整數(shù)。

令P=2r，則P是可以表示的。再利用mpz_ui_pow_ui（M，2，P）計(jì)算一個(gè)大整數(shù)M=2P。

顯然，若s=0則m=r且M=2P=22r即為所計(jì)算的費(fèi)馬數(shù)。

若s>0，采用先記錄M=22r，比如記錄到T，即T=22r。下面置P=2w，則P也是可以表示的。此時(shí)可以 利用mpz_pow_ui（U，T，P）計(jì)算出 一 個(gè)U=TP=(22r)2w=22r+w。可見，在每次將U置于T循環(huán)s次就能得到U=22r+sw=22m。

到此為止，利用GMP大數(shù)庫在理論上（GMP庫的有效范圍內(nèi)）可以表示任意費(fèi)馬數(shù)了。

（2）長字長系統(tǒng)方案

文獻(xiàn)[1]表明它的試驗(yàn)在Windows XP平臺(tái)用Maple和C/C++測(cè)試止于F23。根據(jù)前小節(jié)的推算，這是極有可能的。因?yàn)閃indows XP是32位操作系統(tǒng)。筆者在64位Win7系統(tǒng)用Maple軟件測(cè)試發(fā)現(xiàn)，計(jì)算F29時(shí)是可行的，但是F30就出現(xiàn)溢出。由此可以推斷，利用長字長OS系統(tǒng)可以計(jì)算更多的費(fèi)馬數(shù)。當(dāng)然這種方法只能取決于計(jì)算機(jī)系統(tǒng)的發(fā)展，與目前分解費(fèi)馬數(shù)采用的網(wǎng)格計(jì)算[17]是不可比的。

（3）分解算法的優(yōu)化方案

文獻(xiàn)[1]給出的分解費(fèi)馬數(shù)過程Procedure4是建立在該文求對(duì)一般形如2au+1因數(shù)的算法上的。分析Procedure 4，可以發(fā)現(xiàn)有些冗余。其實(shí)對(duì)費(fèi)馬數(shù)而言，輸入?yún)?shù)就是n和c，根本不用計(jì)算n=■l og2log2(N-1)■和b=■l og2(N-1)■-χ。這樣一來c就是一個(gè)關(guān)鍵。這也是文獻(xiàn)[1]所指出的問題之一。另外，文獻(xiàn)[1]發(fā)現(xiàn)了局部變量a是有界，因此它給出a=m+20的估算。筆者認(rèn)為該作者沒有注意到文獻(xiàn)[15]里面列出的表2。

表2 參數(shù)a取值與頻率

可見，文獻(xiàn)[1]建議a=m+20可以改為a=m+13。（4）c的選擇方案

根據(jù)文獻(xiàn)[1]的陳述，c是算法里面的一個(gè)輸入?yún)?shù)。前面小節(jié)已經(jīng)述及，這是減少搜索次數(shù)的參數(shù)。文獻(xiàn)[1]有幾個(gè)地方介紹了這個(gè)參數(shù)的來由：

并且指出，若事先能估值u的范圍，那么c原則上應(yīng)滿足u>2χ。

根據(jù)這個(gè)原則，結(jié)合文獻(xiàn)[15]所列舉的各費(fèi)馬數(shù)數(shù)的u值，筆者建議對(duì)于小費(fèi)馬數(shù)，c的上限取值如表3所示（表中除去了F14、F20、F24三個(gè)小于F29缺失u之范圍的數(shù)）。

限于篇幅這里不一一列舉其他。按照這個(gè)規(guī)律以及文獻(xiàn)[15]所介紹u針對(duì)m的分布范圍，選擇合適的c是不成問題的。筆者在Maple里面的實(shí)驗(yàn)結(jié)果也證實(shí)了這個(gè)方案的有效性，詳見后文。

表3 各費(fèi)馬數(shù)c的上限建議值

2 算法與實(shí)現(xiàn)

本小節(jié)根據(jù)前述優(yōu)化方案，給出分解費(fèi)馬數(shù)的優(yōu)化算法，并給出在64位Maple下的實(shí)驗(yàn)結(jié)果。

2.1 算法

在文獻(xiàn)[1]給出的Procedure 4過程中，輸入的參數(shù)是N=22n+1與c，然后通過公式n=l og2log2(N-1)、b=l og2(N-1)-χ和a=n+20來計(jì)算出n及至b和a。事實(shí)上，對(duì)于費(fèi)馬數(shù)而言，這些計(jì)算顯得冗余。另外，原過程中的計(jì)算也可以直接優(yōu)化為l=22n-i-1,r=22n-i+1+1。如此優(yōu)化，還避免了浮點(diǎn)計(jì)算。這里給出優(yōu)化后的算法如下。

算法：解析費(fèi)馬數(shù)小因數(shù)的算法

2.2 Maple實(shí)現(xiàn)的效果

基于前述改進(jìn)的算法，在一臺(tái)配置Intel E5645@2.4GHz CPU、8G內(nèi)存、Win7 64位臺(tái)式機(jī)上進(jìn)行實(shí)驗(yàn)，選擇Maple-64位軟件編程，參照文獻(xiàn)[15]的結(jié)果，一些費(fèi)馬數(shù)尋找其小因數(shù)得到表4的運(yùn)行結(jié)果。

表4 若干費(fèi)馬數(shù)及其分解情況

3 結(jié)語

對(duì)比文獻(xiàn)[1]的表1、文獻(xiàn)[15]的清單以及本文的表4可得到以下結(jié)論：

（1）文獻(xiàn)[1]的算法的確能夠找到費(fèi)馬數(shù)的小因子，對(duì)于較小的u也是非常快的算法，但是也存在優(yōu)化的空間。本文進(jìn)行的優(yōu)化也證明了這一點(diǎn)。例如，文獻(xiàn)[1]花了50496步搜索找到了F6的因子274177而經(jīng)本文優(yōu)化后僅化了28步，文獻(xiàn)[1]花了4034686步搜索找到了F10的因子45592577而經(jīng)本文優(yōu)化后僅化了1474步，文獻(xiàn)[1]花了1606849步搜索找到了F19的因子70525124609而經(jīng)本文優(yōu)化后僅化了435步。

（2）大數(shù)表示和計(jì)算的工具實(shí)為非常緊迫的需要。文獻(xiàn)[1]表示利用GMP大數(shù)庫運(yùn)算在32位系統(tǒng)只能計(jì)算到F22。筆者在64位系統(tǒng)利用Maple也只能計(jì)算到F29。須知，Maple也是采用GMP大數(shù)庫作為其內(nèi)核運(yùn)算工具的。因此，完善大數(shù)表示和計(jì)算的工具實(shí)為非常緊迫的需要。

（3）更多更高效的算法仍需開發(fā)。文獻(xiàn)[1]的算法和本文優(yōu)化的算法，總體上還是適合較小的u。當(dāng)u較大時(shí)，如F7的u為116503103764643，所需O(0.25×116503103764643×214)=O(238598356509988864)的計(jì)算量。采用2.4GHz主頻的計(jì)算機(jī)需要運(yùn)行約99415982秒（27615小時(shí)?3年）才能算出結(jié)果。其次，文獻(xiàn)[1]及本文算法的時(shí)間復(fù)雜度為O(0.125u(log2N)2)。注意到對(duì)于費(fèi)馬數(shù)Fm而言，其計(jì)算量為O(22m-3u)因?yàn)閘og2Fm=2m，依然是巨大的計(jì)算量。筆者沒有研究國際上目前采用的網(wǎng)格計(jì)算方法，但是從其幾個(gè)月甚至1年多才產(chǎn)生一個(gè)新結(jié)果的進(jìn)度來看，優(yōu)秀算法依然可期。

綜上所述，盡管從理論上，文獻(xiàn)[1]和本文算法可以有效析出費(fèi)馬數(shù)的小因數(shù)，更快的算法還是值得研究的。筆者考慮是否在并行和超級(jí)計(jì)算方面能否開發(fā)出合適的算法。限于知識(shí)的制約，希望未來能看到相應(yīng)的成果。