快速素數生成方法綜述

李 峰,龔宗躍,2,雷翻翻,顧 申,高 鵬

1.大唐微電子技術有限公司,北京100094

2.北京航空航天大學,北京100191

3.北華大學 計算機科學技術學院,吉林132013

1 引言

素數在很多密碼算法和協議中都有使用,有些算法的安全性完全取決于素數的保密性[1–4],如RSA[1]中,取素數p和q,計算參數N=pq及私鑰d=e?1modφ(N)(e為選定的公鑰參數); 算法的理論安全性基于大數N分解的困難性,實際的安全性基于素數p和q的保密性及所使用p,q的安全性.

素數生成對很多密碼算法是至關重要的,而由于大素數分布的稀疏性,素數生成是一個很耗時的過程.因此,研究快速高效的素數生成方法是有必要的.

各種素數生成方法的算法流程雖然不完全一致,但是整體思路都是一樣的,流程如圖1所示.

圖1中可看到素數生成分為3 步,分別為選取初值、增量變換、素性檢測,對于快速素數生成算法的研究也從這三步著手進行.

初值p的選取有兩類主要算法,一類是直接選取隨機奇數,另一類是選取與小素數乘積互素的數.文獻[5]提到,確保一個隨機奇數不能被3、5和7 整除,即可在素數檢測前排除54% 的奇數; 確保隨機奇數與小于100 的素數均互素即可排除76% 的奇數.

進行增量變換時,增量函數f的選取是一個關鍵點,不同算法對f的選取不同,可以每次給初值加2、每次加一個素數或經過一個復雜的函數變換.

素性檢測是所有素數生成算法必不可少的一個環節,有兩大類檢測方法,一類是可證明素數檢測,另一類是概率素數檢測.實踐中概率素數檢測算法使用更為廣泛,經概率素數檢測算法通過的素數有時也被稱為“工業級素數”.

2 生成與小素數乘積互素數的算法

本節主要介紹幾種生成與小素數乘積互素數的算法及其相關理論基礎.素數生成算法中初值的生成常用這類算法.

2.1 查表法

定理1(中國剩余定理[5])若m1,m2,··· ,mn兩兩互素且都大于1,a1,a2,··· ,an是任意整數,則存在一個解x滿足同余方程組:

其解一定可以寫成以下形式:

且θi滿足同余方程組

則由中國剩余定理可得

故可得到以下等價條件:

對于一個固定的Π,可以計算得到θi,再選擇k個隨機數xi滿足這樣可以得到一個與Π的每個因子都互素的x.因此,可以得到算法1[6].

該算法需要存儲預計算得到的θi,因此需要較大的存儲空間.實際使用中通常對于選定Π的將θi進行預存,提高算法執行效率.

2.2 模搜索法

定義1(Carmichaelλ函數)對任意整數為素數.λ(N)的定義如下:

其中p為素數.

算法1 查表法Input:Π=∏ki=1pδii Output:與Π 互素的數c 1 預處理:2 計算θi:3 for i=1 to k do()[()?1]4 θi=ΠΠ pδi mod pδi i i pδi i 5 end 6 計算其他參數:t=C ·max(pδi)i ,pδi i 表示pδii 的比特數，C 通常取2.7 主運算:8 c=0,i=1 9 while i≤k do 10 取t 比特的隨機數x 11 if xδiθi mod Π=0 then 12 Goto Step 9 13 end 14 c=(c+xθi)mod Π 15 i=i+1 16 end 17 返回c

定理2(Carmichael 定理[7])設a和N都是正整數,且gcd(a,N)=1,則

其中λ(N)是Carmichael 函數.

由定理2易得推論1.

推論1設a和N都是正整數,若有aλ(N)≡1 modN成立,則有gcd(a,N)=1,其中λ(N)是Carmichael 函數.

由推論1,可得算法2.該算法不需要任何預存值,但要多次計算指數為λ(Π)的模冪運算,算法執行效率較低.

算法2 模搜索法Input:Π=∏k i=1pδii Output:與Π 互素的數c 1 取n 比特隨機數c,使其滿足c < Π 2 while cλ(Π) mod Π=1 do 3 c=c+1 4 end 5 返回c

2.3 改進的模搜索法

定理3令整數k,r∈ZΠ,若gcd(k,r,Π)=1,則有

可以利用定理3對算法2 進行優化,得算法3[8].

算法3 改進的模搜索法Input:Π=∏k i=1pδii ,λ(Π)Output:與Π 互素的數c 1 取隨機數c ∈ZΠ 2 U=(1?cλ(Π))mod Π 3 while U=0 do 4 取隨機數r ∈ZΠ 5 c=c+rU 6 U=(1?cλ(Π))mod Π 7 end 8 返回c

該算法運算量和算法2基本相當,但循環次數有所減少,從而縮短運算時間.

3 素數生成主算法

本節介紹素數生成中最重要的一個步驟,增量函數f的設計.下文描述中使用T(q)表示對q進行素性檢測,當q是素數時返回True,否則返回False.

3.1 原生算法

最原始的生成素數方法就是取一個隨機奇數,進行素性檢測,如果不是素數則重新取一個隨機數,直到找到一個素數.

算法4 原生素數生成算法Input:要生成素數的比特位數n Output:素數q 1 取n 比特隨機奇數q 2 while !T(q)do 3 取n 比特隨機奇數q 4 end 5 返回q

接下來討論理論上該算法所需素數檢測次數.

定義2(素數分布函數[9])素數分布函數π(x)表示小于或等于x的素數個數,表達式可以寫為且p是素數1,其中x是大于1 的正實數.

定理4(素數定理)π(x)是的漸進,即

定義3(Li(x)函數)Li(x)[10]是一個比更好的π(x)的逼近,令x是大于1 的正實數,則

由于公式(2)需計算極限,計算不是很方便,在有限誤差范圍內,Li(x)可以近似寫成

該式與式(1)僅相差一個常數Li(2)≈1.05.

定理5π(x)是Li(x)的漸進,即

根據定理4和定理5,得到算法4生成滿n比特隨機數時平均素數檢測次數的兩個近似值,分別為

和

表1 算法4平均素數檢測次數Table 1 Average prime test times of Algorithm 4

由表1可知算法4所需素數檢測的次數較多,生成素數平均時間較長.

3.2 增量素數生成算法

算法4效率較低,算法5對其進行改進.該算法由Jorgen Brandt 等[11]提出并在文獻[12]中對其性能進行了全面分析.

算法5 增量素數生成算法Input:要生成素數的比特位數n Output:素數q 1 取n 比特隨機奇數q 2 while !T(q)do 3 q=q+2 4 end 5 返回q

由于素數分布的不均勻性,因此從理論上無法計算出該算法平均需要進行的素數檢測次數,除非已知了指定范圍內所有素數的分布.經過實驗分析,該算法平均素數檢測次數比算法4略有減少,算法5的另一優勢在于能提高素數檢測效率,將在4.2節進行討論.

3.3 改進的增量素數生成算法

對算法5進行兩方面改進,一方面取初值q滿足gcd(q,Π)=1,其中Π是小素數乘積,通常Π=2×3×5×···×29; 另一方面,迭代過程由q=q+2 改進為q=q+Π,這樣可以保證每個q均沒有這些小素因子.

詳細算法描述見算法6.

算法6 改進的增量素數生成算法Input:要生成素數的比特位數n,前k 個小素數乘積Π=∏ki=1pi Output:素數q 1 取n 比特隨機奇數q,使之滿足gcd(q,Π)=1 2 while !T(q)do 3 q=q+Π 4 end 5 返回q

與算法5類似,由于素數分布的非均勻性,該算法需要進行的素數檢測次數也無法得到理論上的確定值,但可以確定該算法的素數檢測次數與參數k直接相關,在一定范圍內,k取值越大所需素數檢測次數越少.步驟1可使用算法1–3實現.

3.4 M-J 素數生成算法

M-J 素數生成算法由Marc Joye 等在文獻[6]中提出.該算法對算法6進行了改進,減少了所需素數檢測的次數.

素數生成一般指定一個生成范圍[wmin,wmax],在密碼算法中一般都是生成固定n比特的素數,通常取wmin=2n?1+1,wmax=2n?1.

取參數Π滿足Π

且使不等式兩端的值盡量接近,即

為下文描述方便,記ψ=εmaxΠ,τ=εminΠ.M-J 算法流程見算法7.

算法7 M-J 素數生成算法Input:要生成素數的比特位數n,ψ和τ Output:素數q 1 取隨機數c,使c ∈Z?ψ 2 計算q=c+τ 3 while !T(q)do 4 c=fa(c)5 q=c+τ 6 end 7 返回q

Marc Joye 在文獻[6]中給出生成n比特素數時算法7的時間復雜度為O(n4/lnn),但未給出完整的證明.

考慮算法實現的效率,可以令步驟4中a=2,即ci=2ci?1modψ,乘2 的運算可以通過左移位快速實現,同時乘2 后的取模運算可以轉換為大數減運算,即ci=2ci?1或ci=2ci?1?ψ.由于有限制條件故取a=2時需要限制gcd(ψ,2)=1,則Π中不能包括因子2,εmax為奇數.因為ψ是奇數,所以步驟5得到的候選值qi可能是偶數.當qi為偶數時,給其加Π保證進行素數檢測的qi全部為奇數.由此得到算法7的一個特例,流程見算法7A[6].

算法7A M-J 素數生成算法的特例Input:要生成素數的比特位數n,ψ和τ和Π(不含因子2)Output:素數q 1 取隨機奇數c,使c ∈Z?ψ 2 計算q=c+τ 3 if q mod 2=0 then 4 q=q+Π 5 if !T(q)then 6 c=2c mod ψ 7 Goto Step 2 8 返回q

實際中常使用算法7A,可有效縮短生成素數的時間,提高效率.

3.5 改進的M-J 素數生成算法

給定目標素數范圍[wmin,wmax],首先生成算法參數.選一個小數0 <ε≤1(一般取10?2或者10?3),選擇小素數的乘積則存在整數t,u,v滿足以下條件:

(1)1?ε<(uΠ?1)/(wmax?wmin)≤1;

(2)vΠ+t≥wmin;

(3)(u+v)Π+t?1≤wmax;

(4)Π包括盡量多的不同的素因子且Π

根據參數(wmin,wmax,ε)得到參數(Π,t,u,v),即可進行素數生成,見算法8.

令ψ=uΠ,τ=vΠ,取隨機數且a=1.

算法8 改進的M-J 素數生成算法Input:t,ψ,τ和a Output:素數q 1 取隨機數k ∈Z?ψ 2 計算q=[(k?t)mod ψ]+t+τ 3 if !T(q)then 4 k=a·k mod ψ 5 Goto Step 2 6 返回q

算法8生成素數實際的分布范圍為[vΠ+t,(u+v)Π+t?1]?[wmin,wmax],是目標素數范圍的一個子集.該子集與目標范圍的接近程度由參數ε決定.

該算法可以保證進行素數檢測的每一個q均滿足gcd(q,Π)=1.該算法通用性較好,輸出素數的分布也比算法7更好.文獻[8]給出了算法8生成n比特素數時所需素數檢測次數為n·ln 2·?(Π)/Π=O(n/lnn).

令u=1,這樣一次取值的區間長度就變成了Π?1,為了使取隨機數的范圍盡可能接近目標范圍,將t的取值由固定值改為Π的隨機倍數,即t=bΠ,每次取數時b在變化.此外,令a=2,步驟4只需要進行移位和減運算,速度有明顯提升;a=2時必須使才能保證成立此時,上述的4 個條件將變成如下形式:

(1)1?ε<((bmax?bmin+1)Π?1)/(wmax?wmin)≤1;

(2)(v+bmin)Π≥wmin;

(3)(v+1+bmax)Π?1≤wmax;

(4)Π包括盡量多的不同的素因子且Π

此時就得到算法8的一個特例,算法8A.

算法8A 改進的M-J 素數生成算法的特例Input:Π=∏i pi(pi=2),τ=vΠ,bmin,bmaxOutput:素數q 1 取隨機數k ∈Z?Π 2 取隨機數b ∈[bmin,bmax],計算t=bΠ 3 計算q=k+t+τ 4 if q%2=0 then 5 q=Π?k+t+τ 6 if !T(q)then 7 k=2·k mod Π 8 Goto Step 2 9 返回q

算法8A中步驟1可以使用第2節中的3 個算法進行運算,步驟7可以通過移位和減實現,避免使用模乘,保證算法的高效性.應用中經常使用算法8A.

算法8中一些參數取特定值時還可以進一步簡化,令參數u=1,t=0,a=65537,步驟1中k=65537,則算法可以簡化為

其中α,v是滿足要求的隨機數,計算得到q后進行素數檢測,如果不通過則重新選擇α和v.根據文獻[13]該方法是某主流密碼設備廠商使用的RSALib 中的方法,文獻[16]提出一種針對該算法的攻擊——RoCA 攻擊,使RSA 的密鑰在較短時間內被恢復.

4 素數檢測算法

素數檢測算法主要分為可證明素數檢測和概率素數檢測.

可證明素數檢測是基于完整的數學推導,理論上可保證通過檢測的肯定是素數,常用的可證明素數檢測方法有n?1 分解法、Jacobi和測試法、基于橢圓曲線的素數檢測法等,這些方法運算量都較大,耗時較長,在一般生成用于密碼算法密鑰參數時使用較少,本文不詳述這類算法,可參考文獻[10]和[13].

概率素數檢測算法在實際中使用更為廣泛,這類算法相比于可證明素數檢測算法計算量大大降低,但是其給出的結果有一定的可能性是錯誤的,即輸入一個素數,這類算法給出的結果肯定是素數; 輸入一個合數,算法檢測后可能認為該數是素數.

4.1 概率素數檢測算法

常見概率素數檢測類算法有Fermat 檢測、Solovay-Strassen 檢測和Miller-Rabin 檢測三種算法.

Fermat 檢測是基于費馬小定理,但是該類檢測算法會將費馬偽素數誤認為是素數,而且需要進行模冪運算,運算量較大,實際中基本已經不再使用.

Solovay-Strassen 檢測的理論基礎是歐拉準則,但是該檢測算法會將歐拉偽素數誤認為是素數,同時計算中需要計算“雅可比記號”(Jacobi Symbol),計算較為復雜.由于歐拉偽素數比費馬偽素數少很多,因此該檢測算法的正確率比Fermat 高.但是在Miller-Rabin 算法提出后該算法的使用少了很多.

算法9 Miller-Rabin(n,t)Input:待檢測奇數n ≥3,檢測迭代輪數t Output:n 是否是概率素數1 將n?1 寫成n?1=2sr,其中r 是奇數2 for i=1 to t do 3 取隨機整數a ∈[2,n?2]4 計算y=ar mod n 5 if y=1 且y=n?1 then 6 j=1 7 while j≤s?1 且y=n?1 do 8 y=y2 mod n 9 if y=1 then 10 返回“n 是合數”11 end 12 j ++13 end 14 if y=n?1 then 15 返回“n 是合數”16 end 17 end 18 end 19 返回“n 可能是素數”

Miller-Rabin 的理論依據是定理6.

定理6取素數n,令n?1=2sr其中r是奇數.對任意整數a滿足gcd(a,n)=1,均有ar≡1 modn或者a2/r≡?1 modn對某些j,0≤j≤s?1.

Miller-Rabin 檢測算法對于一些強偽素數[11]可能給出錯誤的判斷.對于任一合數,經過Miller-Rabin(n,t)檢測,認為是素數的概率不超過(1/4)t.

Miller-Rabin 檢測相對于Fermat 檢測、Solovay-Strassen 檢測效率要高很多,但是仍然需要進行模冪運算,因此實際使用中還有一些技巧提高素數檢測的效率.

4.2 素數檢測效率的提高

(1)進行小素數試除

生成的候選值n可能含有小的素因子,因此在進行Miller-Rabin 檢測前可以先進行小素數試除,分別計算Ri=nmodpi其中pi是選定的小素數,若存在Ri=0 則表明n是合數,通過小素數試除的候選值再進行Miller-Rabin 檢測.

(2)針對3.2節方法的進一步優化

若素數生成算法使用算法5,候選值依次加2,則在小素數試除過程中,計算Ri=nmodpi可以進一步提速.對于初始候選值n計算一次Ri=nmodpi將結果進行保存,則n=n+ 2時只需計算Ri=(Ri+2)modpi,可以將大數取模運算變成一次模加運算,有效提高效率.

(3)Miller-Rabin 算法步驟4的優化

算法9中步驟4需要進行底數為a的模冪運算,模冪耗時較長; 但取a=2時,模冪運算可以簡化,2r可以通過移位快速實現.使用Miller-Rabin 檢測時需要取t個不同的a,因此在第一輪迭代中固定a=2可以有效提高整體檢測效率.

5 實驗結果

本節對上述素數生成算法進行實驗測試,對比分析了各算法的效率.測試使用Sage Math 編程,運行在Intel Xeon E5-1620@ 3.60 GHz 處理器的工作站.為了比較不同平臺對算法性能的影響,將部分算法移植到智能卡進行了測試.

5.1 生成與小素數乘積互素數算法性能對比

該小節對算法1–3的性能進行對比,算法生成一個輸出的耗時與小素數個數k強相關,實驗結果如圖2所示.

圖2 生成與小素數乘積互素數算法性能對比Figure 2 Performance comparison of generating a number co-prime with a product of small primes

圖2可看出當k< 60時算法3有著明顯性能優勢,當k> 60時算法1的性能更好一些.但考慮到算法1要存儲預計算值且一般取k<60,后續實驗需要用到生成小素數乘積互素數時都選用算法3.

選擇k=60 的情況在智能卡上進行實現,所需預計算在卡外進行.為了準確分析算法執行時長,采集功耗曲線進行分析(采集功耗曲線可以去除數據傳輸的時間,7816 接口數據傳輸較慢),算法1–3的功耗曲線如圖3所示.

圖3 智能卡上實現算法1–3的功耗曲線對比Figure 3 Comparison of power traces of Algorithm 1–3 on smart card

從圖3中可以明顯看出算法1調用模乘的次數,算法2–3調用模冪的次數.大量采集各算法執行時的功耗曲線,量取各次執行的時長,得到平均執行時長.算法1的平均時長7.5 ms,算法2的平均時長10.5 ms,算法3的平均執行時長6.7 ms,3 個算法相對執行時長與工作站的實驗結果一致.

5.2 素數生成主算法性能對比

本小節對比第3節中介紹的素數生成主算法的性能,分別比較平均生成一個素數所需素數檢測次數和平均時長.

由于算法6的性能與輸入參數中選擇素數的個數強相關,因此先分析素數個數對算法6性能的影響,如圖4所示.

圖4 算法6輸入參數k 的不同取值算法性能對比(生成512 位)Figure 4 Performance comparison of parameter k’s value of Algorithm 6(generate 512 bits prime)

圖4是生成512 bit 素數時,素數個數k取不同值對算法6性能(包括平均素數檢測次數和平均耗時兩方面)的影響,圖中可以看出輸入的素數個數越多,算法性能越好,因此在條件允許時使用盡量多的小素數.但是選擇越多的素數,意味著需要越多的預存空間,因此實際使用時需要對時間和空間進行平衡和取舍.本文選擇性能優先,在算法6性能盡量好的情況下,對比算法4–6、7A和8A五個算法的性能,結果如表2所示,表中數據是10 萬次運算得到的平均值.

從表2可知算法4和算法5性能較差,平均時長和平均素數檢測次數都沒有優勢; 算法7A和算法8A性能基本相當,文獻[8]中提到原理上算法8A生成素數的分布均勻性比算法7A好,因此二者中優先考慮算法8A.表中算法6的素數檢測次數比算法7A、8A多一些,但是平均執行時長更短,分析原因是算法6 的gcd()方法使用Sage Math 的系統函數,效率較高,而算法7A、8A中需要調用算法3,其中模冪使用二進制展開法實現,比Sage Math 的系統函數pow()1分析Python 中函數pow()的底層C 源碼,指數長度不同時算法不同.性能低,此外算法8A中步驟7直接使用的系統模乘函數,未使用移位提速.因此在專用密碼設備上使用相同的底層算子實現算法時,算法8A和算法6時間基本相當,甚至會更快一些,對于生成較短的素數時算法8A的平均素數檢測次數更少.

表2 素數生成主要算法性能對比Table 2 Performance comparison of main step of prime generation

此外需要注意,表中的平均執行時長與處理器結構有關,甚至算法間的相對時長也有差異.使用Intel Core i5 8250 4 核8 線程處理器執行同樣的程序,生成1024 比特的素數時,算法8A的執行時長為5.61 ms而算法7A的平均執行時長為5.80 ms.因此,平均素數檢測次數更值得關注.

6 應用

素數在密碼學領域有著大量的應用,常見的非對稱密碼算法RSA 中使用了大素數.RSA 算法基本原理如下.

選擇兩個大素數p和q,計算n=pq,隨機選取加密密鑰e,使p與(p?1)(q?1)互素,然后計算解密密鑰d=e?1mod((p?1)(q?1)).e和n作為公鑰,d是私鑰,p,q不再需要,通過安全方式將其丟棄.加密消息m時計算c=memodn,解密時計算m=cdmodn.可以通過CRT 方式提高運算速度.為了獲得最大程度的安全性,一般取p和q長度相等,同時還要保證p和q不是很接近.

應用中,通常使用相同的公鑰e(常取3 或65 537),因此生成算法參數p和q時就要保證(p?1)(q?1)與e互素.本節選擇對算法3和算法8進行改進,使之生成的素數p能夠滿足gcd(e,p?1)=1.

(1)e是一個素數.只需在p的素性檢測后判斷pmode是否等于1,若等于1 則重新生成素數.

(2)e的所有素因子ei均滿足ei|Π.

取整數α使之滿足gcd(α,Π)=1 且αei?1≡1 modei對所有ei均成立.記e+=gcd(e,Π).令算法3步驟1的初值c≡αmode+,可以取c=α+re,r為隨機數.再令算法8步驟4中a=α2,因為e+|Π,所以可保證算法8生成p的序列始終滿足p≡α2j+1mode+,其中j為整數.因此對所有的ei均有p≡1 modei成立,故可保證gcd(e,p?1)=1.

(3)e不是素數且存在素因子ei?Π.由(2)可知當ei|Π時,通過對算法3和算法8的部分參數進行特定取值可以保證p≡1 modei,故只需對ei?Π驗證p≡1 modei是否成立.可以使用定理2進行判斷,記只需判斷(p?1)λ(e?)≡1 mode?是否成立,若不成立則重新生成素數.

綜上所述,只要算法3步驟1取c=α+re,算法8步驟4中a=α2,再加一步額外的判斷即可使用算法3和算法8直接生成滿足要求的RSA 的參數p和q.

不同應用中對素數有不同的要求,可以根據相應的要求對算法進行適當的改造,使算法能夠直接生成滿足要求的素數.文獻[6]和[8]給出了生成安全素數、應用于ANSI X9.31、生成強素數等多種不同場景下算法的改造.

7 結論

本文對近年來提出的一些快速素數生成算法進行歸納、總結和對比.將素數生成過程分為初值生成、增量變換和素性檢測3 個過程,分別介紹了每部分的算法,給出了各算法必要的理論基礎,分析了算法的演進過程,并給出一些算法實現中的技巧,提到了部分算法的安全性.

本文通過實驗驗證了各算法的性能,給出了各算法性能的對比數據,由實驗結果可知改進的增量素數生成算法(算法6),改進的M-J 素數生成算法的特例(算法8A)分別結合改進的模搜索法(算法3)這兩種方案整體性能較好.本文還給出了素數生成算法在RSA 算法中的應用.

本文僅給出了各算法的基本形式和性能分析,實際用于密碼設備時還需要考慮算法實現[17]的安全性,考慮是否有側信道信息的泄漏[18],考慮結果分布是否均勻等,用于RSA 算法時還需要保證使用的素數是強素數[19].