改進(jìn)的蒙哥馬利模乘算法及FPGA實(shí)現(xiàn)

程碧倩，劉光柱，肖 昊

(合肥工業(yè)大學(xué) 微電子學(xué)院，安徽 合肥 230009)

隨著互聯(lián)網(wǎng)和云計(jì)算技術(shù)的飛速發(fā)展，網(wǎng)絡(luò)購(gòu)物、網(wǎng)絡(luò)支付、云存儲(chǔ)等互聯(lián)網(wǎng)新經(jīng)濟(jì)迅速發(fā)展，方便了人們的生產(chǎn)生活。但由于互聯(lián)網(wǎng)存在開(kāi)放性，因此保障云數(shù)據(jù)和云交易的安全成為工業(yè)界和學(xué)術(shù)界關(guān)注的焦點(diǎn)。目前，利用公鑰密碼系統(tǒng)(Public Key Cryptosystem，PKC)對(duì)數(shù)據(jù)信息進(jìn)行加密或簽名保護(hù)是保障網(wǎng)絡(luò)信息安全的主要方法之一[1-3]。大整數(shù)模乘運(yùn)算是PKC的核心運(yùn)算，其計(jì)算效率對(duì)PKC的性能至關(guān)重要。蒙哥馬利模乘算法作為一種高效率的模乘算法，能夠避免模乘中的除法運(yùn)算，減少模乘的計(jì)算量，是優(yōu)化大整數(shù)模乘計(jì)算效率的研究重點(diǎn)[4-5]。

目前，國(guó)內(nèi)外關(guān)于蒙哥馬利模乘算法的優(yōu)化已有不少研究成果。按照其研究目的，主要分為兩類(lèi)：一是提高計(jì)算時(shí)間；二是減少邏輯資源。文獻(xiàn)[6～8]采用Karastuba算法或基于該算法的衍生算法對(duì)輸入數(shù)據(jù)進(jìn)行變形，降低計(jì)算的復(fù)雜度，并進(jìn)行全展開(kāi)的256 bit模乘運(yùn)算，大幅提高了計(jì)算速度，縮短了計(jì)算時(shí)間。但對(duì)于大位寬模乘(如1 024 bit、2048 bit等)，該方法需消耗大量邏輯資源。文獻(xiàn)[9]采用蒙哥馬利模乘算法，減少了計(jì)算的邏輯資源，但其每次運(yùn)算僅能處理2 bit乘法，因此需經(jīng)過(guò)多次迭代完成模乘運(yùn)算。盡管該算法節(jié)省了大量的邏輯資源，但其速度仍有待提升[9]。

針對(duì)大位寬模乘運(yùn)算，為了更好地平衡其計(jì)算時(shí)間和邏輯開(kāi)銷(xiāo)，本文提出了一種多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法，并基于現(xiàn)場(chǎng)可編程門(mén)陣列(Field Programmable Gate Array，F(xiàn)PGA)設(shè)計(jì)實(shí)現(xiàn)了1 024 bit多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法。本文對(duì)所提算法的正確性進(jìn)行了驗(yàn)證，并對(duì)其性能進(jìn)行了對(duì)比分析。

1 算法介紹

假設(shè)有3個(gè)長(zhǎng)度均為nbit的二進(jìn)制大整數(shù)A、B和M，其中A與B范圍在0～M。模乘運(yùn)算表示為S=A×BmodM，通常情況下，先計(jì)算A、B的乘積然后再計(jì)算modM。經(jīng)典蒙哥馬利模乘算法用簡(jiǎn)單的移位運(yùn)算代替耗時(shí)的取模運(yùn)算。本文在此基礎(chǔ)上提出了一種靈活度高、低成本、高效的多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法。

1.1 經(jīng)典蒙哥馬利模乘算法

經(jīng)典蒙哥馬利模乘算法的設(shè)計(jì)思想是利用剩余系性質(zhì)，借助構(gòu)造一個(gè)模M的剩余系，將普通求模運(yùn)算轉(zhuǎn)化成移位和加法運(yùn)算[10-13]。蒙哥馬利算法完成S=A×B×R-1modM。選擇滿(mǎn)足0≤A，BM)，并選取滿(mǎn)足0

R×R-1=1 modM

(1)

M′=-M-1modR

(2)

經(jīng)典蒙哥馬利算法的實(shí)現(xiàn)，主要經(jīng)過(guò)以下4個(gè)步驟：

步驟1將大整數(shù)A和B相乘，得到乘積T，稱(chēng)此步驟為模乘法操作；

T=A×B

(3)

步驟2取乘積T的低n位，并與模數(shù)的逆M′相乘，取乘積的低n位，得到商數(shù)q

q=((TmodR)×M′)modR

(4)

步驟3將商數(shù)q乘以模數(shù)M，加上乘積T，取結(jié)果的高n位，稱(chēng)此步驟為模約簡(jiǎn)操作

S=(T+qM)/R

(5)

步驟4判斷(S>M)是否成立，若成立返回(S-M)，否則返回S。

通過(guò)對(duì)經(jīng)典蒙哥馬利算法分析并結(jié)合式(2)得到

q×M=T×M′×M≡-TmodR

(6)

且S×R≡TmodR，即S≡T×R-1modM。由于0≤T+q×M

可以看到，在經(jīng)典蒙哥馬利算法中，邏輯操作包括兩次nbit的大整數(shù)乘法和兩次2nbit的大整數(shù)加法，位寬越大，邏輯操作越復(fù)雜，耗時(shí)越長(zhǎng)。因此，本文主要通過(guò)兩方面實(shí)現(xiàn)對(duì)算法的改進(jìn)：(1)對(duì)輸入數(shù)據(jù)進(jìn)行多項(xiàng)式展開(kāi)確定每次參與運(yùn)算的數(shù)據(jù)位寬；(2)針對(duì)模乘法操作和模約簡(jiǎn)操作的開(kāi)展方式進(jìn)行調(diào)整。

1.2 多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法

本文對(duì)經(jīng)典蒙哥馬利模乘算法做了改進(jìn)，得到多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法，其由兩個(gè)外循環(huán)組成，稱(chēng)為循環(huán)1和循環(huán)2。其中，循環(huán)1求取模約簡(jiǎn)運(yùn)算中的中間參數(shù)qi；循環(huán)2經(jīng)多次迭代完成模乘法運(yùn)算和模約簡(jiǎn)運(yùn)算，并在其內(nèi)嵌的小循環(huán)(循環(huán)3)中完成單次迭代的模乘法運(yùn)算及模約簡(jiǎn)運(yùn)算。

在該算法中，輸入?yún)?shù)：模數(shù)M、乘數(shù)A、被乘數(shù)B，其多項(xiàng)式表達(dá)分別為

k2=n/y

(7)

0≤A≤2M，k2=n/y

(8)

0≤B≤2M，4M<2k1x

(9)

輸入?yún)?shù)：M′0，它為模數(shù)M的逆(即M′)的低xbit，S0= 0。

輸出參數(shù)：S。

該算法步驟如下：

步驟1依次取B的xbit(即bi)乘以A的低xbit(即A[0])，并加上中間結(jié)果Si的低xbit(即Si[0])，完成乘累加運(yùn)算后保留低xbit并乘以M′0，然后取低xbit得到參數(shù)qi：

For(i=0 tok1-1)//循環(huán)1

If(i=0)ThenSi=0//初始化

步驟(1)qi=(((Si[0]+A[0]×bi)modx)×M′0)

modx

End for

步驟2本算法采用由循環(huán)2和循環(huán)3組成的二重循環(huán)。在每次內(nèi)循環(huán)(循環(huán)3)中，每次取A的ybit(即aj)依次乘以B的xbit(即bi，其中i取0～30)，并累加上一次循環(huán)結(jié)果Si的ybit(即Si，j)和進(jìn)位C_ABi+1，j完成乘累加運(yùn)算，得到模乘法結(jié)果。隨后，將qi乘以模數(shù)M的ybit(即mj)，并和上一級(jí)進(jìn)位C_qmi+1，j以及S_ABi1，j+1做連續(xù)累加，完成模約簡(jiǎn)運(yùn)算：

For(i=0 tok1-1)//循環(huán)2

If(i=0)ThenSi=0//初始化

步驟(2)For(j=0 tok2-1)//循環(huán)3

If(j=0)ThenC_ABi+1，j=0；//初始化

步驟(2a)(C_ABi+1，j+1，S_ABi+1，j+1)=bi×

aj+Si，j+C_ABi+1，j

步驟(2b)(C_qmi+1，j+1，S_qmi+1，j+1)=qi×

mj+S_ABi+1，j+1+C_qmi+1，j

End for

Si+1=S_qmi+1/x

End for

步驟3本算法通過(guò)擴(kuò)大輸入數(shù)據(jù)范圍，省去經(jīng)典蒙哥馬利模乘算法的減法操作[14-15]。計(jì)算結(jié)束，輸出結(jié)果S。

本算法的改進(jìn)主要集中在兩方面：(1)多項(xiàng)式展開(kāi)執(zhí)行邏輯運(yùn)算的數(shù)據(jù)位寬；(2)模乘中涉及的乘法和模約簡(jiǎn)操作開(kāi)展方式。

1.2.1 多項(xiàng)式展開(kāi)

假設(shè)nbit的大整數(shù)B以rB為基表示，B=(bk1-1，bk1-2，…，b1，b0)rB，取rB=2x，每部分用bi表示，bi范圍在0～x，通過(guò)式(7)計(jì)算得到。其中，i取0～k1-1，整數(shù)x和整數(shù)k1滿(mǎn)足4M< 2^(k1x)，得到k1bit的大整數(shù)B(式(8))。

bi=(B>>(x×i))modrB

(10)

(11)

另外，k1bit的大整數(shù)B也可按照(k1-1)次的多項(xiàng)式展開(kāi)，如式(9)所示，其中，dB=rB。

(12)

同理，假設(shè)nbit的大整數(shù)A以rA為基表示，A=(ak2-1，ak2-2，…，a1，a0)rA，取rA=2y，每部分用aj表示，aj范圍在0～y，通過(guò)式(10)計(jì)算得到。其中j取0～k2-1，k2=[n/y]，則得到A為

aj=(A?(dA×j))modrA

(13)

(14)

此外，k2bit的大整數(shù)A也可按照(k2-1)次的多項(xiàng)式展開(kāi)，如式(12)，其中dA=rA。

(15)

T=A×B=A(dA)×B(dB)=

(16)

若輸入的數(shù)據(jù)為nbit的乘數(shù)A和被乘數(shù)B，計(jì)算過(guò)程中最多會(huì)產(chǎn)生2nbit的中間結(jié)果；若對(duì)輸入數(shù)據(jù)以多項(xiàng)式形式展開(kāi)，執(zhí)行多精度計(jì)算，每次僅取部分A、B(即aj與bi)執(zhí)行乘法運(yùn)算，產(chǎn)生的中間結(jié)果縮小至(x+y)bit。在硬件實(shí)現(xiàn)中，該乘法器可被復(fù)用，其大小將影響電路的使用面積，執(zhí)行本操作有利于節(jié)省面積開(kāi)銷(xiāo)。

1.2.2 模乘法和模約簡(jiǎn)交叉執(zhí)行

本文利用上述多項(xiàng)式展開(kāi)模乘算法中涉及的乘法和模約簡(jiǎn)運(yùn)算，并結(jié)合多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法里的步驟2進(jìn)行舉例說(shuō)明。步驟2主要由循環(huán)2和循環(huán)3組成，循環(huán)2內(nèi)嵌循環(huán)3，循環(huán)3完成式(14)和式(15)的運(yùn)算，分別實(shí)現(xiàn)模乘法(步驟(2a))和模約簡(jiǎn)操作(步驟(2b))。循環(huán)2的計(jì)算如式(16)所示。

=bi×A

(17)

=qi×M+S_ABi+1

(18)

=q×M+B×A

(19)

經(jīng)典的蒙哥馬利模乘算法先進(jìn)行模乘法操作，再完成模約簡(jiǎn)運(yùn)算。本文提出數(shù)據(jù)以多項(xiàng)式展開(kāi)來(lái)交叉執(zhí)行模乘法和模約簡(jiǎn)運(yùn)算，相比經(jīng)典的蒙哥馬利模乘算法，本文所提算法每次僅執(zhí)行小位寬的模乘法運(yùn)算，便可開(kāi)始執(zhí)行模約簡(jiǎn)運(yùn)算，縮短了模約簡(jiǎn)開(kāi)始計(jì)算的時(shí)間。若采用兩組乘法單元分別完成模乘法和模約簡(jiǎn)功能，可實(shí)現(xiàn)多組數(shù)據(jù)并行執(zhí)行，進(jìn)一步提高了大整數(shù)模乘計(jì)算效率。

綜上，基于上述兩方面改進(jìn)得到的模乘算法更適合大整數(shù)模乘運(yùn)算的硬件實(shí)現(xiàn)。其每次僅以小位寬交替執(zhí)行模乘法操作和模約簡(jiǎn)操作，提高了計(jì)算效率。新算法通過(guò)減少乘法器的面積節(jié)省了硬件資源，實(shí)現(xiàn)低資源消耗的高效模乘。

2 算法實(shí)現(xiàn)與性能分析

為了驗(yàn)證所提算法的正確性及性能，本文對(duì)其進(jìn)行FPGA實(shí)現(xiàn)。本節(jié)內(nèi)容主要由兩部分組成：(1)介紹基于多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法的計(jì)算流程及其硬件結(jié)構(gòu)；(2)先對(duì)本文提出算法的正確性進(jìn)行硬件仿真，再將本算法和經(jīng)典蒙哥馬利模乘算法性能進(jìn)行對(duì)比，最后與其它文獻(xiàn)基于FPGA實(shí)現(xiàn)的模乘算法進(jìn)行性能對(duì)比。

2.1 算法功能實(shí)現(xiàn)

圖1為多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法計(jì)算流程圖。當(dāng)輸入數(shù)據(jù)給定后，初始化Si(即S0=0)和C_ABi+1，j(即C_ABi+1，0=0)，然后進(jìn)入k1次迭代運(yùn)算，迭代完成后輸出結(jié)果，即結(jié)束運(yùn)算。

圖1 多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘計(jì)算流程圖

圖2為改進(jìn)算法的硬件結(jié)構(gòu)圖，由多個(gè)乘累加單元、乘法器構(gòu)成。計(jì)算過(guò)程主要分為3個(gè)階段進(jìn)行：(1)第1階段。將B的xbit(即bi)與A的低xbit(即A[0])以及部分結(jié)果Si的低xbit(即Si[0])作為輸入，經(jīng)過(guò)k1次迭代，在乘累加單元1里完成乘累加運(yùn)算，并保留運(yùn)算結(jié)果的低xbit，和乘累加單元1的另一輸入M′0相乘，保留低xbit得到參數(shù)qi；(2)第2階段。以ybit為單位從A的低位到高位得到aj，乘以bi并疊加上輪迭代結(jié)果Si，j(其中Si，j為Si的ybit，j取0～k2-1)和上一級(jí)進(jìn)位C_ABi+1，j，經(jīng)過(guò)k1次迭代，在乘累加單元2內(nèi)完成對(duì)應(yīng)的模乘法運(yùn)算；(3)第3階段。將第1階段得到的商參數(shù)qi依次與模數(shù)M的ybit(即mj)相乘，并加上模乘法結(jié)果S_ABi+1，j+1和上一級(jí)進(jìn)位C_qmi+1，j，在乘累加單元3內(nèi)完成模約簡(jiǎn)操作。經(jīng)過(guò)這3個(gè)階段的計(jì)算，實(shí)現(xiàn)1 024 bit的多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法。

圖2 多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘的硬件結(jié)構(gòu)

2.2 性能與分析

本文使用Verilog HDL對(duì)所提多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法進(jìn)行RTL級(jí)實(shí)現(xiàn)，并使用Vivado2017.4進(jìn)行功能仿真和邏輯綜合。圖3為其中一組隨機(jī)產(chǎn)生的測(cè)試數(shù)據(jù)(十六進(jìn)制)及其對(duì)應(yīng)的計(jì)算結(jié)果。該結(jié)果與圖4的仿真結(jié)果一致，證明了本文算法的正確性。其中A和B為輸入數(shù)據(jù)，M為模數(shù)，M′0為模數(shù)M的逆的低xbit。圖4的start信號(hào)拉高，啟動(dòng)模乘計(jì)算，done信號(hào)拉高，計(jì)算結(jié)果為S。

圖3 測(cè)試數(shù)據(jù)及結(jié)果

圖4 測(cè)試仿真圖

為了更好地量化多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法的性能，將本文提出的算法和經(jīng)典蒙哥馬利算法進(jìn)行性能對(duì)比，如表1所示。其中，存儲(chǔ)空間由輸入數(shù)據(jù)和運(yùn)算單元的資源開(kāi)銷(xiāo)兩部分組成。對(duì)于經(jīng)典蒙哥馬利模乘算法，設(shè)輸入數(shù)據(jù)均為nbit，故A、B、M、M′消耗4n存儲(chǔ)空間，乘法器位寬為2n(代表A×B或q×M)，加法器位寬為2n+1，存儲(chǔ)器可被復(fù)用，故運(yùn)算單元的存儲(chǔ)空間占用(2n+1)bit，因此共消耗(4n+2n+1)bit存儲(chǔ)空間。對(duì)于本文提出的算法，輸入數(shù)據(jù)拆分成xbit和ybit(其中x=[(n+2)/k1]，y=[n/k2])，故A、B、M、M′0共消耗(3n+x)bit存儲(chǔ)空間，兩個(gè)乘法器位寬分別為x+y(代表aj×bi)與2x(代表A[0]×bi)，兩個(gè)加法器位寬分別為x+y+2(代表步驟(2a)或(2b))與2x+1(代表步驟(1))，因此存儲(chǔ)空間共消耗(3n+x+2×(x+y+2)+2x+1)bit(其中2×(x+y+2)為乘累加單元2和3的存儲(chǔ)開(kāi)銷(xiāo)之和，(2x+1)為乘累加單元1的結(jié)果存儲(chǔ)空間)。由于k1和k2的取值通常大于10，因此x和y遠(yuǎn)小于n，故表1中本算法的3個(gè)性能參數(shù)均遠(yuǎn)小于經(jīng)典的蒙哥馬利模乘算法。綜上所述，本文提出的算法相比經(jīng)典蒙哥馬利算法，節(jié)省了大量的硬件資源。

表1 兩種算法性能對(duì)比

本文實(shí)驗(yàn)在Xilinx Virtex7 FPGA上實(shí)現(xiàn)，芯片型號(hào)是XC7VX485T，驗(yàn)證選取參數(shù)x=48，y=34，dA=248，dA=234，k1=22，k2=31，提出的多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法計(jì)算1次1 024 bit的模乘運(yùn)算，經(jīng)過(guò)353個(gè)周期，在237 MHz工作頻率下耗時(shí)1.62 μs，消耗2 317個(gè)LUT和131個(gè)DSP。

基于FPGA實(shí)現(xiàn)的不同蒙哥馬利模乘算法性能的比較如表2所示，其中LUT和DSP的使用數(shù)量代表了算法的邏輯開(kāi)銷(xiāo)情況，反映了相應(yīng)算法的復(fù)雜度；面積時(shí)間積AT1(AT1=LUT的使用數(shù)量×?xí)r間×10-3)和AT2(AT2=DSP的使用數(shù)量×?xí)r間×10-3)用以評(píng)估算法的計(jì)算效率和邏輯開(kāi)銷(xiāo)的折衷性能[12]。文獻(xiàn)[16]使用了27個(gè)DSP和13 811個(gè)DSP，完成1 024 bit模乘運(yùn)算耗時(shí)3.60 μs。本文相比文獻(xiàn)[16]提高55%的計(jì)算速度，節(jié)省了90% LUT數(shù)量，減少了93.2%的AT1，僅增加1倍的面積時(shí)間積AT2。文獻(xiàn)[17]使用22 625個(gè)LUT，耗時(shí)1.01 μs完成1 024 bit模乘運(yùn)算。本設(shè)計(jì)相比于文獻(xiàn)[17]節(jié)省了89.8%的LUT資源，減少了85.2%的面積時(shí)間積AT1，僅增加131個(gè)DSP。本文相比文獻(xiàn)[18]減少96.5%的AT1，減少69%的AT2。經(jīng)上述分析可見(jiàn)，本文提出的算法相比其它參考文獻(xiàn)，可更好地平衡1 024 bit蒙哥馬利模乘運(yùn)算的計(jì)算時(shí)間和面積開(kāi)銷(xiāo)。

表2 不同蒙哥馬利模乘算法的性能比較

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

本文分析了經(jīng)典的蒙哥馬利模乘算法，提出了一種改進(jìn)的蒙哥馬利算法，即多項(xiàng)式展開(kāi)的交叉蒙哥馬利模乘算法。該算法通過(guò)對(duì)數(shù)據(jù)以多項(xiàng)式展開(kāi)來(lái)交叉執(zhí)行模乘法和模約簡(jiǎn)運(yùn)算，可實(shí)現(xiàn)低成本、高效的模乘。本文基于FPGA對(duì)所提算法進(jìn)行了驗(yàn)證，結(jié)果表明本文算法具有良好的靈活性和通用性，不僅可以在FPGA平臺(tái)上實(shí)現(xiàn)，也可應(yīng)用到其他平臺(tái)。