關(guān)沛冬, 羅玉琴, 張方國(guó), 田海博
1. 中山大學(xué)計(jì)算機(jī)學(xué)院, 廣州 510006
2. 廣東省信息安全技術(shù)重點(diǎn)實(shí)驗(yàn)室, 廣州 510006
有限域上基于橢圓曲線的公鑰密碼學(xué)是1985 年由Koblitz[1]和Miller[2]提出的. 從那時(shí)起, 有限域上的橢圓曲線被應(yīng)用在許多密碼系統(tǒng)和密碼協(xié)議中, 如Diffie-Hellman 密鑰協(xié)議方案[3]、橢圓曲線數(shù)字簽名算法[4]等. 橢圓曲線密碼體制的安全性是基于橢圓曲線離散對(duì)數(shù)問(wèn)題(ECDLP) 這一困難問(wèn)題. 解決ECDLP 意味著破解了ECC, 因此這一困難問(wèn)題受到了各界的廣泛關(guān)注.
1997 年, Certicom[5]提出了ECC 挑戰(zhàn), 其本質(zhì)是求解ECDLP. 有限域大小小于100 位的挑戰(zhàn)被當(dāng)作練習(xí)被快速解決, 有限域大小大于100 位的挑戰(zhàn)則分為兩種級(jí)別. 第一級(jí)是109 位和131 位的挑戰(zhàn); 第二級(jí)是163 位、191 位、239 位和359 位的挑戰(zhàn). 迄今為止, ECC2K-108、ECC2-109、ECCp-109 已經(jīng)被解決. Bailey 等人[6]曾通過(guò)網(wǎng)上大量征集計(jì)算資源嘗試攻擊ECC2K-130, 但是最后沒(méi)有給出成功攻擊的結(jié)果. ECC2-131 等更困難的挑戰(zhàn)則迄今沒(méi)有被解決.
Pollard[7]在1978 年基于生日攻擊提出了一種求解一般循環(huán)群離散對(duì)數(shù)的算法, 這種算法被稱為rho 算法. 設(shè)G是一個(gè)循環(huán)群, 則Pollard rho 算法的迭代函數(shù)是一個(gè)函數(shù)F:G →G, 而且此函數(shù)從一個(gè)初始點(diǎn)P0出發(fā), 依次產(chǎn)生Pi+1=F(Pi), 其中i= 0,1,···. 其中通過(guò)映射函數(shù)產(chǎn)生的元素序列為P0,P1,···. 由于Pi ∈G, 而循環(huán)群的階|G| 是有限的, 因此該序列一定會(huì)產(chǎn)生一個(gè)元素與前面的某個(gè)元素相同. 這種現(xiàn)象稱為碰撞或者匹配. 當(dāng)碰撞產(chǎn)生后, 將有很大概率求解出離散對(duì)數(shù)問(wèn)題的解.
對(duì)于求解ECDLP, 由Gallant 等人[8]和Wiener 等人[9]所修改的兩種Pollard rho 算法是迄今最有效的通用求解算法. Van Oorschot 等人[10]提出了一種并行化的Pollard rho 算法, 這種算法通過(guò)并行計(jì)算的方法能線性提高求解離散對(duì)數(shù)問(wèn)題的效率, 而且適用于現(xiàn)代計(jì)算機(jī)的計(jì)算結(jié)構(gòu), 能高效地用程序?qū)崿F(xiàn). 分布式Pollard rho 算法是目前大家公認(rèn)的最有效的解決一般循環(huán)群離散對(duì)數(shù)的算法.
自從分布式Pollard rho 算法提出以來(lái), 多位密碼學(xué)家在該基礎(chǔ)上提出了各種優(yōu)化算法和改進(jìn)方案.而這些工作都主要都是對(duì)Pollard rho 算法的迭代函數(shù)進(jìn)行優(yōu)化. Teske[11]通過(guò)在迭代函數(shù)中使用更多的點(diǎn)加運(yùn)算, 提出了高效的r-加游走迭代函數(shù). Teske 還提出并分析了r+q-混合游走迭代函數(shù), 他在r-加游走迭代函數(shù)中引入了倍點(diǎn)運(yùn)算. 然而,r+q-混合游走會(huì)降低一定的隨機(jī)性, 并且在橢圓曲線的仿射坐標(biāo)下, 倍點(diǎn)運(yùn)算的效率不比點(diǎn)加運(yùn)算的效率高. Bessalov[12]首先在ECDLP 中使用了半分的思想, 但是他沒(méi)有對(duì)該算法進(jìn)行詳細(xì)分析. Zhang 和Wang[13]使用橢圓曲線點(diǎn)半分運(yùn)算為Pollard rho 算法提出了一個(gè)新的迭代函數(shù)r+h-混合游走, 但是在分析的過(guò)程中沒(méi)有考慮橢圓曲線點(diǎn)加運(yùn)算的高效實(shí)現(xiàn)對(duì)他們所提算法的影響.
本文針對(duì)ECC2-131, 對(duì)基于r-加游走,r+q-混合游走和r+h-混合游走這三種不同迭代函數(shù)的分布式Pollard rho 算法進(jìn)行快速實(shí)現(xiàn), 并對(duì)其算法性能進(jìn)行比較, 發(fā)現(xiàn)迭代函數(shù)為r-加游走的分布式Pollard rho 算法在理論與程序?qū)崿F(xiàn)上具有優(yōu)勢(shì). 本文還給出了基于ECC2-131 對(duì)三種分布式Pollard rho 算法的軟件實(shí)現(xiàn), 以及給出了針對(duì)有限域F2131的最優(yōu)不可約多項(xiàng)式. 使用該不可約多項(xiàng)式進(jìn)行編程實(shí)現(xiàn)能對(duì)有限域F2131模運(yùn)算的運(yùn)算效率提高11.29%, 對(duì)乘法運(yùn)算的運(yùn)算效率提高11.23%. 最后本文給出了基于r-加游走的分布式Pollard rho 算法在工作站和在天河二號(hào)超級(jí)計(jì)算機(jī)中運(yùn)行的效率, 發(fā)現(xiàn)在超級(jí)計(jì)算機(jī)中需要消耗約1296.41×106核時(shí), 花費(fèi)約1.296 億元才能攻破ECC2-131. 由此可見(jiàn)當(dāng)前ECC2-131 仍然是個(gè)無(wú)法有效解決的問(wèn)題.
令Fq是有q個(gè)元素的有限域, 定義在有限域Fq上的橢圓曲線E是由一個(gè)非奇異的Weierstrass 方程所定義的. 該非奇異Weierstrass 方程為:
設(shè)E是一條定義在有限域上的橢圓曲線. 設(shè)點(diǎn)P是一個(gè)階為n的點(diǎn). 定義群G是由點(diǎn)P生成的n階循環(huán)群. 對(duì)于任意的Q ∈G, 存在整數(shù)s,0≤s <n, 使得Q=sP成立, 其中sP表示s個(gè)點(diǎn)P相加.因此, 給定橢圓曲線E以及其中滿足Q=sP的點(diǎn)P,Q, 求出整數(shù)s就是橢圓曲線離散對(duì)數(shù)問(wèn)題, 簡(jiǎn)稱ECDLP.
Pollard rho 算法的工作原理是,首先定義一個(gè)周期性循環(huán)的元素序列,然后在序列中尋找兩個(gè)“碰撞”的元素. 而分布式Pollard rho 算法則是在Pollard rho 算法的基礎(chǔ)上篩選一些具有某些性質(zhì)的可區(qū)分元素, 然后只在這些可區(qū)分元素上面去尋“碰撞”. 分布式Pollard rho 算法的兩個(gè)關(guān)鍵思想是生成序列化的迭代函數(shù)和“碰撞” 檢測(cè)算法. 本節(jié)將介紹分布式Pollard rho 算法以及三個(gè)比較具有代表性的迭代函數(shù).
若gcd(bi-bj,n)=1, 則s=(aj-ai)(bi-bj)-1modn, 從而求解出ECDLP.
Pollard rho 算法在迭代的過(guò)程中迭代函數(shù)會(huì)產(chǎn)生大量的點(diǎn). 為了減少迭代點(diǎn)的存儲(chǔ), 只存儲(chǔ)可區(qū)分點(diǎn)是一種常見(jiàn)的手法, 存儲(chǔ)過(guò)程中通過(guò)篩選滿足一定條件的點(diǎn), 能減少大量迭代點(diǎn)的存儲(chǔ). 定義群G中的一個(gè)子集D, 其中D包含了所有的滿足特殊條件的點(diǎn)的集合, 例如具有固定數(shù)量的前導(dǎo)零的元素. 集合D中的點(diǎn)被稱為可區(qū)分點(diǎn). 分布式Pollard rho 算法的核心思想在于多個(gè)程序同時(shí)執(zhí)行迭代函數(shù), 實(shí)現(xiàn)快速大量的隨機(jī)游走以加速碰撞的發(fā)生.
接下來(lái)將會(huì)介紹在Pollard rho 算法的三種改進(jìn)的迭代函數(shù), 分別是Teske[14]提出的r-加游走和r+q-混合游走, 還有Zhang[13]等人提出的r+h-混合游走.
2.2.2r-加游走迭代函數(shù)
根據(jù)Teske[14]提出的r-加游走(r-adding walk) 迭代函數(shù), 隨機(jī)選取2 個(gè)整數(shù)序列, 每個(gè)序列有r個(gè)整數(shù):
2.2.4r+h-混合游走迭代函數(shù)
根據(jù)Zhang[13]等人提出的r+h-混合游走(r+h-mixed walk) 迭代函數(shù), 隨機(jī)選取2 個(gè)整數(shù)序列,每個(gè)序列有r個(gè)整數(shù):
ECC2-131 是Certicom[5]公司設(shè)立的橢圓曲線離散對(duì)數(shù)的挑戰(zhàn), 直至目前該挑戰(zhàn)還沒(méi)有被攻破. 本文將根據(jù)此挑戰(zhàn)設(shè)計(jì)出最優(yōu)的程序?qū)崿F(xiàn)方案.
ECC2-131 的詳細(xì)參數(shù)見(jiàn)文獻(xiàn)[16], 它所在的有限域F2131的不可約多項(xiàng)式為f(z) =z131+z13+z2+z+1.為了令分布式Pollard rho 算法在求解ECC2-131 的過(guò)程中達(dá)到最高效率,實(shí)現(xiàn)分布式Pollard rho 算法的程序使用SIMD 指令對(duì)有限域F2131的基本運(yùn)算進(jìn)行優(yōu)化. 在有限域元素?cái)?shù)據(jù)結(jié)構(gòu)的設(shè)計(jì)上,每個(gè)元素由2 個(gè)128 位的整數(shù)組成. 若e ∈F2131, 把e寫(xiě)成多項(xiàng)式的形式
記為(e130,e129,··· ,e0). 因此使用一個(gè)128 位整數(shù)存儲(chǔ)(e127,e126,··· ,e0), 使用另一個(gè)128 位整數(shù)存儲(chǔ)(0,··· ,0,e130,e129,e128), 其中e130前有125 個(gè)0.
有限域F2131的快速模運(yùn)算參考文獻(xiàn)[17] 中所使用的算法, 實(shí)現(xiàn)模運(yùn)算的時(shí)候把有限域的元素使用多個(gè)64 位的整數(shù)表示. 與此同時(shí), 模運(yùn)算僅是乘法運(yùn)算的一個(gè)部分, 因此要執(zhí)行模運(yùn)算的多項(xiàng)式的最高次冪不大于260. 因此最多只需要5 個(gè)64 位的整數(shù)即可表示要執(zhí)行模運(yùn)算的多項(xiàng)式, 記為e=(E[4],E[3],E[2],E[1],E[0]). 模運(yùn)算運(yùn)行的結(jié)果則使用3 個(gè)64 位整數(shù)表示. 在程序?qū)崿F(xiàn)過(guò)程中, 將兩個(gè)64 位的整數(shù)封裝成一個(gè)128 位的整數(shù), 然后使用Intel SSE 指令進(jìn)行移位、異或等操作, 使程序具有更高的運(yùn)行效率. 具體的模運(yùn)算算法見(jiàn)算法1.
算法1 F2131 模運(yùn)算Input: e = (E[4],E[3],E[2],E[1],E[0]).Output: b = e mod f(z).1 for i = 4 →3 do 2T ←E[i];3E[i-3] ←E[i-3]⊕(T ?61)⊕(T ?62)⊕(T ?63);4E[i-2] ←E[i-2]⊕(T ?10)⊕(T ?1)⊕(T ?2)⊕(T ?3);5E[i-1] ←E[i-1]⊕(T ?54);6 end 7 T ←(E[2]&0xFFFFFFFFFFFFFFFC);8 B[0] ←E[0]⊕(T ?10)⊕(T ?1)⊕(T ?2)⊕(T ?3);9 B[1] ←E[1]⊕(T ?54);10 B[2] ←E[2]&0x7;11 return b ←(B[2],B[1],B[0]);
有限域F2131的快速乘法運(yùn)算可以通過(guò)Karatsuba 算法[15]減少多項(xiàng)式相乘時(shí)所需要執(zhí)行的乘法次數(shù). 在編程實(shí)現(xiàn)上, 通過(guò)使用大部分CPU 都支持的PCLMULQDQ 指令[18], 可以在幾個(gè)CPU 周期內(nèi)完成64 位特征為2 的有限域的乘法運(yùn)算. 具體乘法運(yùn)算如算法2 所示.
算法2 F2131 乘法Input: a = (A[2],A[1],A[0]),b = (B[2],B[1],B[0]).Output: c = ab.1 T0 ←A[1]⊕A[2];2 T1 ←B[1]⊕B[2];3 T2 ←A[0]⊕A[2];4 T3 ←B[0]⊕B[2];5 T4 ←A[0]⊕A[1];6 T5 ←B[0]⊕B[1];7 T0 ←PCLMULQDQ(T 0,T 1);8 T1 ←PCLMULQDQ(T 2,T 3);9 T2 ←PCLMULQDQ(T 4,T 5);10 T3 ←PCLMULQDQ(A[0],B[0]);11 T4 ←PCLMULQDQ(A[1],B[1]);12 T5 ←PCLMULQDQ(A[2],B[2]);13 T6 ←T0 ⊕T4 ⊕T5;14 T7 ←T2 ⊕T4 ⊕T3;15 T8 ←T1 ⊕T3 ⊕T4 ⊕T5;16 c ←(T5,T6,T8,T7,T3);17 return c;
算法3 F2131 開(kāi)平方根Input: e = (e130,e129,··· ,e0).Output: c = ■e.1 eodd ←(e129,e127,··· ,e1);2 eeven ←(e130,e128,··· ,e0);3 c ←■z ×eodd mod f(z);4 c ←c ⊕eeven;5 return c;
有限域F2131的解一元二次方程的具體算法在文獻(xiàn)[19] 中被詳細(xì)地描述. 具體解一元二次方程組如算法4 所示. 令Tr(e) 表示有限域元素的跡(trace), 其中e ∈F2m. 則
令H(e) 表示有限域元素的半跡(half-trace), 其中e ∈F2m. 則跡與半跡的計(jì)算方法在文獻(xiàn)[19] 中有詳細(xì)的介紹.
算法4 F2131 解一元二次方程Input: e = (e130,e129,··· ,e0).Output: z2 +z = e+Tr(e) 的解c, 其中z 為變量.1 預(yù)計(jì)算: ?i ∈{1,2,··· ,65}, 計(jì)算H(z2i-1);2 c ←0;3 for i = 65 →1 do 5 e ←e ⊕zi;4if e2i = 1 then c ←c ⊕zi;7end 8 end 9 c ←c ⊕∑65 i=1 e2i-1H(z2i-1);6 10 return c;
常見(jiàn)的有限域F2131上的求逆算法有基于歐幾里得算法的擴(kuò)展歐幾里得算法和基于費(fèi)馬小定理的Itoh-Tsujii 求逆算法[20]. 雖然廣義歐幾里得算法的平均復(fù)雜度比Itoh-Tsujii 求逆算法的復(fù)雜度要低, 但是通過(guò)大量的實(shí)驗(yàn)發(fā)現(xiàn)Itoh-Tsujii 求逆算法的性能更加穩(wěn)定, 其平均開(kāi)銷(xiāo)會(huì)比廣義歐幾里得算法的要少.因此本文使用Itoh-Tsujii 求逆算法實(shí)現(xiàn)有限域F2131求逆運(yùn)算.
文獻(xiàn)[17] 中給出了有限域F2131的快速模運(yùn)算算法, 算法1 把有限域求模運(yùn)算分為兩個(gè)部分. 算法的輸入是260 位的一個(gè)整數(shù). 算法的第一部分把192 到260 位的數(shù)通過(guò)(2) 式的關(guān)系移位到低位, 并與低位的數(shù)進(jìn)行異或; 算法第二部分把131 到191 位的數(shù)移位到低位并異或. 有限域上的不可約多項(xiàng)式為f(z)=z131+z13+z2+z+1, 則有
第一部分的運(yùn)算會(huì)把[192,260] 位的數(shù)移位異或到[61,142] 位中, 由(2) 式可知算法的第一部分可分為3 次移位異或運(yùn)算. 第一次運(yùn)算把需要異或到[61,63] 位的數(shù)進(jìn)行移位異或; 第二次運(yùn)算把需要異或到[64,127] 位的數(shù)進(jìn)行移位異或; 第三次運(yùn)算把需要異或到[127,142] 位的數(shù)進(jìn)行移位異或.
當(dāng)不可約多項(xiàng)式有z3這一項(xiàng)的時(shí)候, 192 到260 位的數(shù)會(huì)與64 到132 位進(jìn)行異或. 此時(shí)模運(yùn)算的第一部分沒(méi)有第一次運(yùn)算; 第二次運(yùn)算將會(huì)把192 位到255 位的數(shù)異或到64 到127 位中; 第三次運(yùn)算把256 到260 位中的數(shù)異或到128 到132 位中.
本文找出在此模運(yùn)算下運(yùn)行效率最高的其中一個(gè)最優(yōu)不可約多項(xiàng)式:f(z)=z131+z8+z3+z2+1,并編程實(shí)現(xiàn)了其模運(yùn)算.只要選取包含單項(xiàng)的冪次為3 而且除了z131的單項(xiàng)其它單項(xiàng)的冪次不要超過(guò)54 的不可約多項(xiàng)式, 那么該不可約多項(xiàng)式即為最優(yōu)不可約多項(xiàng)式. 不可約多項(xiàng)式f(z)=z131+z8+z3+z2+1的模運(yùn)算的運(yùn)算時(shí)間與ECC2-131 的不可約多項(xiàng)式的模運(yùn)算的運(yùn)算時(shí)間如表1 所示. 測(cè)試環(huán)境為: Intel Core i7-7940X CPU @ 3.10 GHz; Ubuntu 18.04 操作系統(tǒng). 本文對(duì)乘法運(yùn)算和模運(yùn)算進(jìn)行50 次測(cè)試, 每次測(cè)試測(cè)量200 000 次模運(yùn)算和乘法運(yùn)算所使用的CPU 周期, 并去除大于最小值10% 的數(shù)據(jù)取平均, 從而得到模運(yùn)算和乘法運(yùn)算的運(yùn)算效率.
表1 F2131 上非最優(yōu)不可約多項(xiàng)式與最優(yōu)不可約多項(xiàng)式執(zhí)行200 000 次模運(yùn)算和乘法運(yùn)算平均需要的CPU 周期Table 1 Average CPU cycles required to perform 200 000 modular and multiplicative operations between non-optimal irreducible polynomials and optimal irreducible polynomials on F2131
對(duì)于ECC2-131, 通過(guò)Zierler[21]等人的方法找出不可約多項(xiàng)式為f(z)=z131+z8+z3+z2+1 的有限域與ECC2-131 中的有限域的同構(gòu)映射, 可以誘導(dǎo)出ECC2-131 所使用的橢圓曲線在新的有限域中的同構(gòu)映射, 從而找出橢圓曲線加法群的同構(gòu)映射. ECC2-131 中的橢圓曲線加法群的離散對(duì)數(shù)問(wèn)題與同構(gòu)映射后的加法群的離散對(duì)數(shù)問(wèn)題是一致的, 因此本文在實(shí)現(xiàn)分布式Pollard rho 算法時(shí)使用了這種同構(gòu)映射的方法實(shí)現(xiàn)求解ECC2-131 的加速.
本節(jié)將會(huì)對(duì)3 種基于不同迭代函數(shù)的分布式Pollard rho 算法的效率進(jìn)行分析.
分布式Pollard rho 算法的迭代函數(shù)主要涉及橢圓曲線上的三個(gè)運(yùn)算, 分別是點(diǎn)加、倍點(diǎn)和半分運(yùn)算.其中涉及到的有限域上運(yùn)算為加法、乘法、求逆、平方、開(kāi)平方根和解二次方程. 使用多項(xiàng)式基表示的有限域F2m中的元素的加法運(yùn)算在編程實(shí)現(xiàn)上僅需要少數(shù)異或指令, 遠(yuǎn)小于乘法運(yùn)算的開(kāi)銷(xiāo). 因此在加法運(yùn)算與乘法運(yùn)算數(shù)量差不多的情況下, 可以選擇性地忽略加法運(yùn)算的開(kāi)銷(xiāo). 設(shè)乘法運(yùn)算的開(kāi)銷(xiāo)為M, 求逆運(yùn)算的開(kāi)銷(xiāo)為I, 平方運(yùn)算的開(kāi)銷(xiāo)為S, 開(kāi)平方根運(yùn)算的開(kāi)銷(xiāo)為Sqr, 解二次方程的運(yùn)算為Sol. 本文把乘法運(yùn)算的時(shí)間作為基準(zhǔn), 然后測(cè)量出有限域上求逆、平方、開(kāi)平方根和解二次方程這四種運(yùn)算與乘法運(yùn)算的時(shí)間比例, 從而在分析三種分布式Pollard rho 算法效率時(shí)只用乘法運(yùn)算的開(kāi)銷(xiāo)表示. 記:
其中i,s,sqr,sol∈R.
基于迭代函數(shù)為r-加游走、r+q-混合游走與r+h-混合游走的分布式Pollard rho 算法只是迭代函數(shù)不相同, 其初始化過(guò)程都是相同的. 分布式Pollard rho 算法具體流程如下.
已知點(diǎn)P,Q滿足P=sQ. 點(diǎn)P生成的群G的階為n. 首先隨機(jī)生成2r個(gè)正整數(shù)mi,ni ∈{1,2,··· ,n-1}, 計(jì)算出r個(gè)點(diǎn)Mi. 接下來(lái)隨機(jī)生成2 個(gè)序列, 每個(gè)序列中有256 個(gè)元素, 這兩個(gè)序列記為(a0,a1,··· ,a255),(b0,b1,··· ,b255), 其中ai,bi ∈{1,2,··· ,n-1},i=0,1,··· ,255. 用這兩個(gè)序列計(jì)算出256 個(gè)橢圓曲線E上的點(diǎn)Pi=aiP+biQ. 令這256 個(gè)點(diǎn)作為初始點(diǎn), 然后執(zhí)行迭代函數(shù), 直到碰撞發(fā)生. 迭代函數(shù)每生成一個(gè)可區(qū)分點(diǎn), 則判斷是否有碰撞發(fā)生.
r-加游走迭代函數(shù)計(jì)算迭代點(diǎn)Pi需要使用橢圓曲線的點(diǎn)加運(yùn)算.r-加游走迭代函數(shù)還需要計(jì)算滿足Pi=aiP+biQ的兩個(gè)整數(shù)ai,bi.橢圓曲線的點(diǎn)加運(yùn)算需要用到1 個(gè)求逆、2 個(gè)乘法、1 個(gè)平方運(yùn)算.通過(guò)使用Montgomery trick[22]的方法計(jì)算求逆運(yùn)算可以快速求出N個(gè)有限域元素的逆, 把原本需要NI的計(jì)算開(kāi)銷(xiāo)轉(zhuǎn)換為1I+(3N-3)M的計(jì)算開(kāi)銷(xiāo). 本文把快速求出N個(gè)有限域元素的逆Montgomery trick的方法稱為同時(shí)逆算法. 實(shí)現(xiàn)r-加游走迭代函數(shù)的程序中, 定義了N=256 個(gè)起始點(diǎn), 迭代函數(shù)的每一次循環(huán)對(duì)這N個(gè)點(diǎn)分別進(jìn)行一次迭代. 之所以要這樣設(shè)計(jì), 是因?yàn)辄c(diǎn)加運(yùn)算中每個(gè)點(diǎn)都需要進(jìn)行一次求逆運(yùn)算, 而且在仿射坐標(biāo)下求逆運(yùn)算的計(jì)算開(kāi)銷(xiāo)遠(yuǎn)大于乘法運(yùn)算的計(jì)算開(kāi)銷(xiāo). 通過(guò)同時(shí)逆的方法能對(duì)求逆運(yùn)算的運(yùn)算效率有大幅度的提升. 通過(guò)這種設(shè)計(jì), 迭代函數(shù)每一次循環(huán)就生成了N個(gè)迭代點(diǎn). 具體r-加游走迭代函數(shù)的偽代碼如算法5 所示:
算法5 r-加游走迭代函數(shù)Input: N = 256; N 個(gè)初始點(diǎn)Pi,i ∈{0,1,··· ,N -1}; 常數(shù)r; r 個(gè)點(diǎn)Mj,j ∈{0,1,··· ,r-1}.Output: ai,bi,aj,bj 滿足Pi = aiP +biQ,Pj = ajP +bjQ, 其中Pi = Pj,i /= j.1 for i = 0 →(N -1) do 2vi ←Pi.x mod r;3si ←Mvi.x+Pi.x;4 end 5 同時(shí)逆算法求出s-1 i ,i = 0,1,··· ,N -1;6 for i = 0 →(N -1) do 7利用s-1 i 計(jì)算Pi = Pi +Mvi;8ai ←ai +mvi;9bi ←bi +nvi;10if Hamming(Pi.x) <= 28 then 11DList.append(Pi,ai,bi) ;12判斷DList 中是否有(Pj,aj,bj) 滿足Pi = Pj 且i /= j, 若滿足, return ai,bi,aj,bj;13end 14 end 15 goto 1
算法5 中ai,bi的計(jì)算為兩個(gè)數(shù)相加然后模n, 因?yàn)閍i,bi <n, 所以ai+bi <2n. 因此模運(yùn)算可以簡(jiǎn)化為整數(shù)減法運(yùn)算. 則計(jì)算ai,bi最多需要2 個(gè)整數(shù)加法運(yùn)算、2 個(gè)判斷和2 個(gè)整數(shù)減法運(yùn)算, 其計(jì)算量仍然小于有限域上的乘法運(yùn)算的計(jì)算量. 由上述分析以及迭代算法偽代碼可得,r-加游走迭代函數(shù)平均每個(gè)迭代點(diǎn)的計(jì)算開(kāi)銷(xiāo)為:
r+q-混合游走迭代函數(shù)使用橢圓曲線中的點(diǎn)加和倍點(diǎn)運(yùn)算. 由4.2 節(jié)可知點(diǎn)加運(yùn)算的計(jì)算開(kāi)銷(xiāo)為2M+ 1S+ 1I. 橢圓曲線的倍點(diǎn)運(yùn)算需要用到1 個(gè)求逆、2 個(gè)乘法、2 個(gè)平方運(yùn)算, 因此其計(jì)算開(kāi)銷(xiāo)為2M+2S+1I. 與此同時(shí), 倍點(diǎn)運(yùn)算與點(diǎn)加運(yùn)算都需要用到求逆運(yùn)算, 因此同樣可以使用同時(shí)逆的方法進(jìn)行優(yōu)化, 即把求逆運(yùn)算的NI的開(kāi)銷(xiāo)轉(zhuǎn)換為I+(3N-3)M. 在實(shí)現(xiàn)r+q-混合游走迭代函數(shù)的程序中迭代函數(shù)的每一次循環(huán)對(duì)N個(gè)點(diǎn)分別進(jìn)行一次迭代. 此迭代函數(shù)與r-加游走迭代函數(shù)的區(qū)別在于計(jì)算完同時(shí)逆后, 需要使用判斷語(yǔ)句判斷當(dāng)前迭代點(diǎn)執(zhí)行的是倍點(diǎn)還是點(diǎn)加運(yùn)算. 具體r+q-混合游走迭代函數(shù)的偽代碼如算法6 所示:
算法6 r+q-混合游走迭代函數(shù)Input: N = 256; N 個(gè)初始點(diǎn)Pi,i ∈{0,1,··· ,N -1}; 常數(shù)r; r 個(gè)點(diǎn)Mj,j ∈{0,1,··· ,r-1}; 常數(shù)q.Output: ai,bi,aj,bj 滿足Pi = aiP +biQ,Pj = ajP +bjQ, 其中Pi = Pj,i /= j.1 for i = 0 →(N -1) do 2vi ←Pi.x mod (r+q);3si ←Mvi.x+Pi.x;4 end 5 同時(shí)逆算法求出s-1 i ,i = 0,1,··· ,N -1;6 for i = 0 →(N -1) do 7if vi >= r then 8利用s-1 i 計(jì)算Pi = 2Pi;ai = 2ai mod n;10bi = 2bi mod n;11else 9 12利用s-1 i 計(jì)算Pi = Pi +Mvi;13ai ←ai +mvi;14bi ←bi +nvi;15end 16if Hamming(Pi.x) <= 28 then 17DList.append(Pi,ai,bi) ;18判斷DList 中是否有(Pj,aj,bj) 滿足Pi = Pj 且i /= j, 若滿足, return ai,bi,aj,bj;19end 20 end 21 goto 1
根據(jù)文獻(xiàn)[11], Teske 對(duì)自身提出的r+q-混合游走與r-加游走進(jìn)行了詳細(xì)的性能分析. 在r+q-混合游走迭代函數(shù)中,r與q的取值以及兩者的比值不同對(duì)算法的隨機(jī)性與計(jì)算效率都有影響. 對(duì)于r+q-混合游走迭代函數(shù), 需要考慮r,q的取值的對(duì)r+q-混合游走算法平均每個(gè)迭代點(diǎn)的計(jì)算開(kāi)銷(xiāo)的影響, 并忽略有限域加法運(yùn)算以及計(jì)算ai,bi所帶來(lái)的計(jì)算開(kāi)銷(xiāo). 由上述分析以及迭代算法偽代碼可得,r+q-混合游走迭代函數(shù)平均每個(gè)迭代點(diǎn)的計(jì)算開(kāi)銷(xiāo)為:
Zhang[13]等人提出的r+h-混合游走迭代函數(shù)使用了橢圓曲線中的點(diǎn)加和半分運(yùn)算. 橢圓曲線上的半分運(yùn)算需要用到1 個(gè)解二次方程、1 個(gè)開(kāi)平方根和2 個(gè)乘法運(yùn)算. 在計(jì)算ai,bi的時(shí)候, 因?yàn)槌龜?shù)是2,因此可以看成是移位運(yùn)算, 在計(jì)算機(jī)上可以高效實(shí)現(xiàn), 進(jìn)而忽略計(jì)算ai,bi的開(kāi)銷(xiāo). 因此計(jì)算一個(gè)半分運(yùn)算的計(jì)算開(kāi)銷(xiāo)為1Sol+1Sqr+2M. 在實(shí)現(xiàn)r+h-混合游走迭代函數(shù)的程序中, 迭代函數(shù)的每一次循環(huán)同樣是對(duì)N個(gè)迭代點(diǎn)進(jìn)行一次迭代, 但是運(yùn)算與前面兩種迭代函數(shù)有所不同. 程序首先對(duì)每一個(gè)迭代點(diǎn)循環(huán)執(zhí)行半分運(yùn)算, 直到該點(diǎn)的下一次迭代需要執(zhí)行點(diǎn)加運(yùn)算為止. 當(dāng)N個(gè)迭代點(diǎn)都執(zhí)行半分運(yùn)算后, 此時(shí)N個(gè)點(diǎn)都需要執(zhí)行點(diǎn)加運(yùn)算, 可使用同時(shí)逆算法對(duì)點(diǎn)加運(yùn)算進(jìn)行加速. 所有點(diǎn)執(zhí)行完點(diǎn)加運(yùn)算后, 完成迭代函數(shù)的一輪循環(huán). 具體r+h-混合游走迭代函數(shù)的偽代碼如算法7 所示.
可見(jiàn)r+h-混合游走迭代函數(shù)平均每個(gè)迭代點(diǎn)的計(jì)算開(kāi)銷(xiāo)為:
算法7 r+h-混合游走迭代函數(shù)Input: N = 256; N 個(gè)初始點(diǎn)Pi,i ∈{0,1,··· ,N -1}; 常數(shù)r; r 個(gè)點(diǎn)Mj,j ∈{0,1,··· ,r-1}; 常數(shù)h.Output: ai,bi,aj,bj 滿足Pi = aiP +biQ,Pj = ajP +bjQ, 其中Pi = Pj,i /= j.1 for i = 0 →(N -1) do 2vi ←Pi.x mod (r+h);3while vi >= r do 2 Pi;5 4 Pi = 1 ai = (ai&1) ? 1 2(bi +n) : 1 2 ai;6 2 bi;7 2(ai +n) : 1 bi = (bi&1) ? 1 DList.append(Pi,ai,bi) ;8判斷DList 中是否有(Pj,aj,bj) 滿足Pi = Pj 且i /= j, 若滿足, return ai,bi,aj,bj;10end 11vi ←Pi.x mod (r+h);12end 13 end if Hamming(Pi.x) <= 28 then 9 14 同時(shí)逆算法求出s-1 i ,i = 0,1,··· ,N -1;15 for i = 0 →(N -1) do 16利用s-1 i 計(jì)算Pi = Pi +Mvi;17ai ←ai +mvi;18bi ←bi +nvi;19if Hamming(Pi.x) <= 28 then 20DList.append(Pi,ai,bi) ;21判斷DList 中是否有(Pj,aj,bj) 滿足Pi = Pj 且i /= j, 若滿足, return ai,bi,aj,bj;22end 23 end 24 goto 1
4.5.1 實(shí)驗(yàn)環(huán)境
本文對(duì)三種分布式Pollard rho 進(jìn)行了仿真實(shí)現(xiàn). 程序基于Intel 指令集的CPU 進(jìn)行實(shí)現(xiàn), 借助SIMD 指令集對(duì)算法進(jìn)行指令集層面上的優(yōu)化.算法的仿真實(shí)現(xiàn)使用的編程語(yǔ)言是C++,編譯器為GCC 9.3.0, 支持C++14 標(biāo)準(zhǔn); 處理器為: Intel Core i7-7940X CPU @ 3.10 GHz; 操作系統(tǒng)為Ubuntu 18.04.
當(dāng)橢圓曲線上的點(diǎn)的x軸坐標(biāo)的Hamming 重量小于13 時(shí), 定義該點(diǎn)為可區(qū)分點(diǎn). 表2 到表4 給出三種迭代函數(shù)在不同參數(shù)下的隨機(jī)性量化指標(biāo).
表2 r-加游走迭代函數(shù)隨機(jī)性量化指標(biāo)Table 2 Quantification index of randomness of r-adding walk
由表3 和表4 可見(jiàn),當(dāng)q/r和h/r的比例小于1 的時(shí)候,q,h的取值對(duì)迭代函數(shù)的隨機(jī)性影響不大;而當(dāng)比例超過(guò)1 時(shí), 兩種迭代函數(shù)的隨機(jī)性變化相對(duì)明顯. 僅測(cè)量迭代函數(shù)的隨機(jī)性不足以說(shuō)明算法的好壞,還需要測(cè)量出在不同參數(shù)下的分布式Pollard rho 算法的整體效率, 從而得出最優(yōu)的算法參數(shù). 4.5.3 節(jié)將會(huì)給出整體效率的測(cè)量結(jié)果.
表3 r+q-混合游走迭代函數(shù)隨機(jī)性量化指標(biāo)Table 3 Quantification index of randomness of r+q-mixed walk
表4 r+h-混合游走迭代函數(shù)隨機(jī)性量化指標(biāo)Table 4 Quantification index of randomness of r+h-mixed walk
4.5.3 有限域運(yùn)算效率以及最優(yōu)分布式Pollard rho 算法
由于求解ECC2-131 的算法與4.5.2 節(jié)中求解有限域F241上的ECDLP 所使用的算法是相同的, 因此在迭代函數(shù)選擇相同參數(shù)的情況下, 求解ECC2-131 的Pollard rho 算法的隨機(jī)性與4.5.2 節(jié)測(cè)出的相同. 對(duì)于ECC2-131, 考慮到在不同的CPU 上程序的運(yùn)算速度不同, 這可能導(dǎo)致乘法運(yùn)算與求逆、平方等運(yùn)算的運(yùn)行時(shí)間的比值不穩(wěn)定,因此本文在多個(gè)CPU 上進(jìn)行測(cè)試來(lái)證明有限域多種運(yùn)算的比值是穩(wěn)定的.對(duì)有限域F2131的每一個(gè)算術(shù)運(yùn)算進(jìn)行100 次開(kāi)銷(xiāo)測(cè)試, 每一次測(cè)試使該運(yùn)算執(zhí)行1000 000 次, 然后取這100 次測(cè)試的平均值然后除以1000 000 得到每次運(yùn)算所消耗的CPU 周期. 表5 為程序使用多種CPU 測(cè)試所得到的結(jié)果.
表5 不同CPU 上的有限域運(yùn)算的計(jì)算開(kāi)銷(xiāo)Table 5 Computational cost of finite field operations on different CPUs
把乘法運(yùn)算的時(shí)間作為基準(zhǔn), 通過(guò)表5 可以計(jì)算出有限域上求逆、平方、開(kāi)平方根和解二次方程這四種運(yùn)算與乘法運(yùn)算的時(shí)間比例, 具體見(jiàn)表6.
表6 不同CPU 上有限域基本運(yùn)算時(shí)間與乘法運(yùn)算的比值Table 6 Ratio of the basic operation time to multiplication operation in finite field on different CPUs
有限域的解二次方程這個(gè)運(yùn)算需要用到AVX512 指令集, 而Intel Core i7-9700K @3.60 GHz CPU不支持AVX512 指令集而只支持AVX2 指令集, 因此在該CPU 上的解二次方程運(yùn)算所用的時(shí)間比其它兩個(gè)CPU 的要長(zhǎng). 本文對(duì)三種迭代算法的計(jì)算速率的測(cè)試在Intel Core i7-7940X @3.10 GHz CPU 上進(jìn)行測(cè)量. 但是總的來(lái)說(shuō), 由表6 可見(jiàn), 同一套有限域?qū)崿F(xiàn)的代碼不會(huì)因?yàn)镃PU 選擇的不同而導(dǎo)致有限域運(yùn)算的運(yùn)算時(shí)間比例有明顯的差異.
根據(jù)表6 所測(cè)量得到的數(shù)據(jù), 可以對(duì)三種分布式Pollard rho 算法的整體效率進(jìn)行測(cè)量, 從而選出最優(yōu)的情況. 由表6 可得, 式(3) 中的參數(shù)取M= 29.179,i= 97.071,s= 0.737,sqr = 1.564,sol = 2.705. 同時(shí)逆算法同時(shí)計(jì)算的有限域元素的個(gè)數(shù)為N= 256. 測(cè)量的方法為統(tǒng)計(jì)出迭代100 000 000 個(gè)隨機(jī)點(diǎn)所需要的CPU 周期. 為了避免CPU 被其它進(jìn)程影響導(dǎo)致測(cè)量結(jié)果偏大, 本文測(cè)量100 次CPU 周期, 去除大于最小值10% 的值后取平均值作為測(cè)量結(jié)果, 然后把測(cè)量結(jié)果除以100 000 000 得到每一個(gè)迭代點(diǎn)平均計(jì)算開(kāi)銷(xiāo). 由此可得, 三種算法每一個(gè)迭代點(diǎn)平均計(jì)算開(kāi)銷(xiāo)如表7 所示.
表7 三種分布式Pollard rho 算法每一個(gè)迭代點(diǎn)平均計(jì)算開(kāi)銷(xiāo)Table 7 Average computational cost of each iteration point of the three distributed Pollard rho algorithms
由表7 可見(jiàn), 仿真測(cè)試值與理論值相比有一定的誤差, 而且都是仿真測(cè)量值比理論值要大. 這個(gè)誤差的存在是合理的, 因?yàn)殚_(kāi)銷(xiāo)表達(dá)式?jīng)]有包含有限域加法運(yùn)算和計(jì)算ai,bi的開(kāi)銷(xiāo), 而且程序在編程實(shí)現(xiàn)的過(guò)程中會(huì)有額外的開(kāi)銷(xiāo), 比如循環(huán)語(yǔ)句、條件判斷語(yǔ)句和函數(shù)調(diào)用開(kāi)銷(xiāo)等. 與此同時(shí), 仿真測(cè)量值的變化趨勢(shì)沒(méi)有完全按照理論值的變化而變化, 這是由于在不同的r,q,h參數(shù)下, 程序的優(yōu)化情況有所不同. 比如r+q-混合游走迭代函數(shù)在r= 128,q= 1024 的情況下, 相比于r+h-混合游走迭代函數(shù)在h= 64 和h=128 的情況運(yùn)算速率更高, 這是因?yàn)閞=128,h=1024 的參數(shù)更有利于程序的優(yōu)化和執(zhí)行. 從整體效率上看, 基于r+q-混合游走迭代函數(shù)或r+h-混合游走迭代函數(shù)的分布式Pollard rho 算法都是當(dāng)r固定不變時(shí)隨著參數(shù)q和h的增大而整體效率的值增大, 即意味著算法的整體效率更低.
通過(guò)表5 和表6 的有限域運(yùn)算的測(cè)量結(jié)果, 從理論上進(jìn)行分析可知, 基于r-加游走迭代函數(shù)的分布式Pollard rho 算法所需要的每一輪迭代的計(jì)算開(kāi)銷(xiāo)是最小的, 而且在程序?qū)崿F(xiàn)上也能達(dá)到最優(yōu)效果, 其整體效率也是最優(yōu)的. 文獻(xiàn)[13] 中提到在不使用同時(shí)逆算法的情況下, 半分運(yùn)算與倍點(diǎn)運(yùn)算相比有非常大的計(jì)算優(yōu)勢(shì), 因此使用半分運(yùn)算可以獲得非常高的提速效果. 然而, 在使用同時(shí)逆算法對(duì)求逆運(yùn)算進(jìn)行加速的時(shí)候, 半分運(yùn)算的優(yōu)勢(shì)會(huì)被削弱. 理想情況下,r-加游走迭代函數(shù)的迭代時(shí)間為
本節(jié)在實(shí)驗(yàn)仿真層面對(duì)求解ECC2-131 的速度與開(kāi)銷(xiāo)進(jìn)行評(píng)測(cè). 通過(guò)在工作站和天河二號(hào)超級(jí)計(jì)算機(jī)上分別進(jìn)行算法開(kāi)銷(xiāo)的評(píng)測(cè). 自橢圓曲線算法提出以來(lái), ECDLP 的求解就一直在被研究. 與此同時(shí), 計(jì)算機(jī)CPU 指令集的發(fā)展也為ECDLP 的求解帶來(lái)了非常大的幫助. 因此本文選取算法上隨機(jī)性好, 整體效率較優(yōu), 軟件實(shí)現(xiàn)能做到高效實(shí)現(xiàn)的基于r-加游走的分布式Pollard rho 算法求解ECC2-131. 軟件實(shí)現(xiàn)的過(guò)程中r-加游走迭代函數(shù)的參數(shù)根據(jù)表7 中整體效率較優(yōu)的結(jié)果, 并考慮r不能取太大以免導(dǎo)致內(nèi)存訪問(wèn)過(guò)于頻繁從而降低效率, 因此取r=255.
本節(jié)的分布式Pollard rho 算法使用C++ 進(jìn)行實(shí)現(xiàn), 程序支持C++14 標(biāo)準(zhǔn). 分布式Pollard rho算法的實(shí)現(xiàn)借助SIMD 指令集進(jìn)行優(yōu)化. 本節(jié)所使用的工作站的配置為: Intel Core i7-7940X @3.10 GHz CPU; 操作系統(tǒng): Ubuntu 18.04.
超級(jí)計(jì)算機(jī)能執(zhí)行需要大量計(jì)算量的程序. 它主要的優(yōu)勢(shì)是有極大的數(shù)據(jù)存儲(chǔ)容量和非常高的數(shù)據(jù)處理速度, 因此在超級(jí)計(jì)算機(jī)上解決ECDLP 就成為了可能. 本文使用天河二號(hào)超級(jí)計(jì)算機(jī)進(jìn)行求解, 首先計(jì)算出一部分ECC2-131 的可區(qū)分點(diǎn), 測(cè)試在天河二號(hào)中求解ECC2-131 的程序的性能. 然后通過(guò)測(cè)試的結(jié)果分析出使用天河二號(hào)求解ECC2-131 的速度以及可行性.
Daniel 等人[16]的測(cè)量結(jié)果指出,在Core 2 Extreme Q6850 CPU 上進(jìn)行ECDLP 的求解,可以做到單核消耗533 個(gè)CPU 周期運(yùn)行一輪迭代.在他們的軟件程序中,有限域使用多項(xiàng)式基表示,而有限域的數(shù)據(jù)結(jié)構(gòu)使用Bitslice[16]轉(zhuǎn)換得到. 在Bitslice 轉(zhuǎn)換后的數(shù)據(jù)結(jié)構(gòu)下, 有限域的平方、乘法和求逆運(yùn)算與非Bitslice 的形式有很大的不同. 在CPU 層面上, 天河二號(hào)中單個(gè)節(jié)點(diǎn)單核的計(jì)算性能與Core 2 Extreme Q6850 CPU 單核的計(jì)算能力相近. 因此忽略CPU 性能對(duì)兩種算法所帶來(lái)的影響, 本文實(shí)現(xiàn)的程序在天河二號(hào)中能做到單核單線程消耗418.2163 個(gè)CPU 周期完成一輪迭代, 與Daniel 等人的結(jié)果相比有21.5%的效率提速.
考慮到CPU 技術(shù)的不斷發(fā)展, 本文在Intel Core i7-7940X @3.10 GHz CPU 上能做到每個(gè)迭代花費(fèi)229.8463 個(gè)CPU 周期, 與Daniel 等人的結(jié)果相比有了56.9% 的效率提升. 而且本文算法的軟件實(shí)現(xiàn)對(duì)于CPU 的性能要求比較高, 用到了Intel CPU 中的PCLMULQDQ 指令, 近幾年CPU 不斷的發(fā)展, PCLMULQDQ 指令所消耗的時(shí)間越來(lái)越小, 這也是為什么同一套軟件在Intel Core i7-7940X @3.10 GHz CPU 和在天河二號(hào)中性能有明顯的差距.
本文重新分析了基于r-加游走、r+q-混合游走和r+h-混合游走迭代函數(shù)的三種分布式Pollard rho算法的計(jì)算效率, 并對(duì)這三種算法在軟件上進(jìn)行快速實(shí)現(xiàn). 本文發(fā)現(xiàn)基于迭代函數(shù)為r-加游走的分布式Pollard rho 算法在理論分析和仿真測(cè)試中與其它兩種算法相比有著更高的效率.
與此同時(shí), 本文在有限域F2131中找出最優(yōu)不可約多項(xiàng)式, 通過(guò)域同構(gòu)誘導(dǎo)出橢圓曲線的同構(gòu)從而得出橢圓曲線群的同構(gòu). 使用同構(gòu)后的有限域和橢圓曲線進(jìn)行編程實(shí)現(xiàn), 程序的有限域乘法運(yùn)算有11.23%的提速, 模運(yùn)算有11.29% 的提速. 本文在Intel Core i7-7940X @3.10 GHz CPU 和天河二號(hào)超級(jí)計(jì)算機(jī)上對(duì)求解ECC2-131 的軟件程序進(jìn)行效率測(cè)試, 發(fā)現(xiàn)目前為止, 計(jì)算ECC2-131 仍然是一個(gè)困難的問(wèn)題.即使使用天河二號(hào)所有CPU 計(jì)算資源, 也需要消耗12.96 億核時(shí), 花費(fèi)1.296 億元, 運(yùn)算約125.60 天才能計(jì)算出ECC2-131. 這在時(shí)間和經(jīng)濟(jì)上不實(shí)際, 可見(jiàn)ECC2-131 目前仍然無(wú)法被解決.