面向多橢圓曲線的高速標量乘法器設計與實現

于斌，黃海，劉志偉，趙石磊，那寧

（哈爾濱理工大學軟件與微電子學院，黑龍江 哈爾濱 150080）

1 引言

隨著處理器算力的不斷增強，傳統加密算法的安全性正受到日益嚴峻的考驗，需要更加復雜的加密算法來保證數據的安全。橢圓曲線加密算法是近年來比較活躍的一種非對稱加密算法，與對稱加密算法中的高加密標準（AES,advanced encryption standard）一樣在密碼系統中得到廣泛使用[1]。在最新的傳輸層安全性協議1.3 版（TLS1.3,transport layer security protocol version 1.3）中，使用多種橢圓曲線進行數據加密，無論是在數量上還是在實際使用頻率上，都凸顯了橢圓曲線加密日益重要的位置[2]。

橢圓曲線的運算速度和對多種曲線的兼容一直是研究的重點。Kudithi 等[3]利用可編程邏輯器件（FPGA,field programmable gate array）完成了適用于物聯網的加密算法結構，利用較低面積完成了對P224 和P256 橢圓曲線的支持；Hossain 等[4]設計了同時適用于專用集成電路和FPGA 的標量乘操作步驟，同樣支持了2 條曲線；Liu 等[5]優化了非相鄰形式（NAF,non-adjacent-form）算法并完成了P256 位的標量乘設計；Lee 等[6]采用雙處理單元改進了自右向左的標量乘算法；Chung 等[7]改進了蒙哥馬利模乘的流水結構，提升了標量乘的整體速度。

本文設計了橢圓曲線硬件加密中使用的一種標量乘法器，選取最常用的Curve25519 和secp256r12曲線，使其盡可能復用乘法單元完成標量乘，使用了256 位乘法器得到乘法結果并采用快速模約簡來達到提高模乘速度的目的。此外，考慮標量乘的使用環境，尤其是在簽名和驗簽中的使用情況，本文設置了對于基點和普通點不同的算法，同時支持多標量乘，使設計的標量乘法器能夠滿足橢圓加密的各種應用需求。

2 研究背景

橢圓曲線有很多分類，使用較早的一類曲線是Weierstrass 曲線，滿足式(1)所示的曲線方程。

TLS1.3 所選用由美國國家標準與技術研究所（NIST,National Institute of Standards and Technology）提出的secp256r1、secp384r1 和secp521r1 這3 條曲線，就屬于Weierstrass 曲線。

橢圓曲線上的運算采用模運算，包括模加減、模乘和模逆運算，最核心的運算方式為模乘和模逆。點加和倍點運算以模運算為基礎，并進一步用來實現標量乘。標量乘的性能也是衡量橢圓曲線加密性能的一個重要的指標，為了得到更快的加密速度，需要盡可能快速地完成標量乘。根據標量乘的運算過程，大致可把提高標量乘速度的優化方法分為3 類：優化標量乘算法[8-13]、轉換坐標系[14-16]和優化模運算算法[17-20]。

第一類優化方法是改變標量乘kP的運算過程，降低總體運算量，其中，P是橢圓曲線上的點，k是待乘的倍數。優化的方法包括：使用自左向右和自右向左2 種算法，對于n位的k值需要計算n次倍點和次點加，其中自左向右算法的點加和倍點可以并行運算，但是需要額外的模運算單元[8]；把k做相對簡單的預處理，實現NAF 算法，把點加的次數降低為次[9]；對k做m位的分段處理，實現窗口NAF 法，預計算需要計算窗口內的2m個點運算結果，預計算量較大[10]；使用Comb 算法來快速完成標量乘，但需要的預計算量龐大[11]；采取適用于多標量乘的Shamir 算法，用一次標量乘的時間完成2 個標量乘的相加[12]；采用多基鏈算法，增加三倍點、五倍點形成多基鏈來加速標量乘[13]。

第二類方法采用轉換坐標系，目的是優化標量乘中的點加和倍點運算。原曲線方程使用二元坐標系，點加運算需要一次模逆和3 次模乘（模平方記為模乘，忽略模加減運算，下文相同），倍點運算需要一次模逆和4 次模乘，而模逆運算需要消耗大量時鐘周期，因此可以通過坐標系的轉換改變點加和倍點運算計算式，消除模逆運算，例如射影坐標系，把坐標(x,y)按轉換式轉換為(X,Y,Z)，點加運算需要14 次模乘，倍點運算需要12 次模乘[14]。雅可比坐標系使用比較普遍，轉換式為這需要16次模乘來完成點加運算，10次模乘完成倍點運算[15]。此外還有修正的雅可比坐標系[16]、混合坐標系[16]等變換，其目的都是在開始時做一次變換，隨后在計算點加和倍點過程中消除模逆運算，最后的計算結果再轉換為二元坐標，整個標量乘只使用一次模逆運算。

第三類優化方法是從基本的模運算入手，提高模運算的速度，主要是模乘和模逆運算。模乘運算使用較多的是蒙哥馬利模乘算法，通過將多位的模乘分解為2m位的模乘，利用二進制的特點，進行多次迭代運算得到最后的模乘結果[17]。IEEE 給出對于NIST 中所使用的模p快速約簡公式，把2 個數據的乘法結果直接做快速的模約簡，得到模乘結果[18]。模逆運算的計算式一般會采用費馬小定理，如式(2)所示。

通過簡單整理即可得到模逆運算的計算式為

這樣可以不需要模逆運算，只需要使用模乘即可完成模逆運算，也可以采用專門的模逆模塊來完成模逆運算。蒙哥馬利模逆算法與蒙哥馬利模乘配合使用，不需要經過轉換也可以完成模逆運算[19]，也可以使用二進制右移算法，通過多次簡單比較和加減來計算模逆運算，資源開銷較小[20]，獨立設計模逆運算主要是為了釋放模乘單元，達到并行的目的，同時減少計算周期數。

還有一類橢圓曲線是蒙哥馬利曲線，其方程為

Curve25519 曲線屬于蒙哥馬利曲線，由Bernstein[21]于 2005 年提出，并在 2016 年與Curve448 一起成為國際標準[22]，這2 種曲線也被TLS1.3 采用。此類曲線提出較晚，在近年開始得到廣泛使用，其標量乘本身具有抗SPA（simple power analysis）特性，所以一般直接使用蒙哥馬利階梯算法而不會做修改[23]，只是在模乘部分做優化，Düll等[24]計算給出了標量乘所需的模乘和模平方數。

上述算法雖各有優點，但都是僅考慮單一使用而優化的。在橢圓曲線加密中，標量乘有很多使用場合，涉及基點G或者普通點P的標量乘，在驗簽時需要計算多標量乘λP+μG，同時服務器端簽名驗簽頻繁，對標量乘速度的要求極高，需要一個能在各種使用場合下都具備較高性能的標量乘法器。

本文全面考慮標量乘法器在橢圓加密算法中的使用情況，完成了如下工作。

1) 選擇適用于不同情況的標量乘算法，使標量乘計算普通點P、基點G和多標量乘λP+μG時都達到較快的運算速度。

2) 推導了Curve25519 的快速模約簡公式，結合已有的secp256r1 的快速模約簡公式，用快速模約簡實現模乘。

3) 結合利用模約簡計算模乘的運算特點，排布了secp256r1 和Curve25519 這2 條曲線的點加倍點操作步驟，根據Comb 算法特點簡化點加操作，并排布了雅可比坐標系下特殊點加（有一個點的Z值為1 時）的操作步驟。

4) 在保證模乘速度的前提下，設計充分考慮復用以降低面積，完成適用于secp256r1和Curve25519曲線的標量乘法器。

3 標量乘法器設計

標量乘法器的設計也同優化一樣，需要完成標量乘算法的設計、點加倍點的設計和模運算的設計。在實際使用情況中，secp256r1 曲線和Curve25519 曲線使用頻繁且具有相近的位數，故將這2 條曲線的標量乘集成為一個單元進行設計。由于采用了3 種標量乘算法來適應不同需求，整個設計采用了3 個獨立的子控制模塊來實現不同算法的控制功能，在工程中可以根據不同的使用場合來切換算法，調用運算單元，達到最快的速度。標量乘法器的整體結構如圖1 所示。

圖1 標量乘法器的整體結構

運算單元完成所有點加倍點運算，除必需的控制邏輯之外，主要運算部分是一個256 位的乘法器，通過快速模約簡單元完成模乘功能，所有預計算數值和中間計算結果都存儲至臨時寄存器堆中，最大限度地復用了硬件資源。由于標量乘計算結果需使用模逆運算轉換至二元坐標系，考慮到橢圓加密的運算過程，將模逆端口一并引出，可供外部調用，使整個設計可以滿足橢圓加密中所有關鍵步驟的使用需求。

3.1 標量乘算法

使用secp256r1 曲線時，加密系統會在簽名和公鑰生成時使用基點G的標量乘，在密鑰交換時使用普通點P的標量乘，在驗簽時使用λP+μG的多標量乘，所以選擇了Comb 算法和Shamir 算法來應對不同使用場合。算法1 為編碼長度為4 的Comb算法[11]。

算法1編碼長度為4 的Comb 算法

輸入256 位二進制數λ={λ255λ254λ253…λ1λ0}和橢圓曲線secp256r1 的基點G。

輸出標量乘Q=λG

Comb 算法需要的預計算量非常大，例如{1111}G就需要計算(2192+2128+264+20)G，在計算普通點時該算法并不占優勢。但是由于基點G已知，Comb 算法的預計算表可以提前完成并存儲，這樣在計算時直接使用該表，再使用64 次倍點和最多64 次點加即可完成基點G的標量乘，需要額外付出的代價是15 個點坐標的存儲陣列，其中，{0000}G不需要存儲。

Shamir 算法可以計算λP+μG的多標量乘，經過預計算后可以用一次標量乘的時間完成λP+μG的計算。算法2 為窗口為2 的Shamir 算法[12]。

算法2窗口為2 的Shamir 算法

輸入256位二進制數λ={λ255λ254λ253…λ1λ0}，μ={μ255μ254μ253…μ1μ0}，橢圓曲線secp256r1上的普通點P和基點G。

輸出Q=λP+μG

當λμ≠0 時，對于輸入的P、G點坐標，需要3 次倍點和4 次點加來計算預計算表，預計算結束后進入循環，4 倍點不設置專用的計算式，復用倍點計算式2 次來完成4 倍點功能，這樣可與Comb 算法復用倍點運算單元，所以在窗口為2 的Shamir 算法中，依然會計算256 次倍點和最多128 次點加，同時存儲除了(00)P+(00)G之外的15 個點坐標。如果只算普通點乘的λP，則令μ=0 即可，此時進行一次點加和倍點即可完成預計算，標量乘循環計算時所需次數與λμ≠0 時相同。

對于Curve25519 曲線采用蒙哥馬利階梯算法來完成標量乘，維持抗SPA 特性，其算法步驟如算法3 所示[23]。

算法3Curve25519 標量乘算法

輸入255 位二進制數k，點坐標P1=(x1,y1)的橫坐標x1

輸出標量乘結果kP1=(x2,y2)的橫坐標x2

其中，cswap 操作是根據swap 的值來交換2 個輸入的坐標值，算法不單獨列出；階梯算法Ladder step在3.2 節中給出詳細步驟。

3.2 高速點加和倍點運算設計

運算單元的功能是執行點加和倍點運算，按2 條曲線劃分，執行的計算式不同。對于橢圓曲線secp256r1，采用的標量乘算法中點加操作相對較少，故采用倍點較快的雅可比坐標系，倍點計算式[16]為

由于采用快速模約簡進行模乘，模乘結果需要經過一個周期的乘法和一個周期的模約簡才可以被繼續使用，為減少模乘次數，并結合式(5)，將倍點運算的操作步驟排布如表1 所示，共需9 個周期完成倍點運算(X3,Y3,Z3)=2(X1,Y1,Z1)。

表1 倍點運算步驟

表1 的步驟5 中，t5得到2S，t1得到M，步驟7的X3即為T值，其余各值依次代入即可驗證是否與式(5)一致。為了保證數據的有效性，模乘結果需間隔一周期才能作為輸入，例如步驟2 的乘法結果t2需在步驟3 進行模約簡，在步驟4 得到模乘結果并繼續參與計算，模約簡過程并未列在算法中。模加減單元作了特殊設計，增加了3 倍模加、被減數為1.5 倍的模減，保證乘法器的連續使用并減少模乘的次數。同時考慮256 位乘法的時間時延較大，而模約簡的時間時延相對較小，所以在模約簡之后可以安排一次模加減，進一步壓縮操作步驟，如步驟2 中t1的使用就是如此。表1、表2 和表3 均采用同樣的設計思路排布，模運算具體內容見3.3 節。

正常情況下表1 的倍點運算需要8 個周期模乘加一個周期的約簡等待，但倍點運算后立即會進行點加或者倍點運算，所以通過調整(X,Y,Z)的輸出順序，可以在下次乘法操作同時進行模約簡，實際工作中需要8 個周期可完成一次倍點運算。

對于點加公式，需要做分類處理以達到最優化，在Shamir 算法中需構造預計算表，但是表內數據都是在計算開始后得到的，所以都是普通的三元坐標(X,Y,Z)，按式(6)排布操作步驟如表2 所示，共需 17 個周期完成普通點加運算(X3,Y3,Z3)=(X1,Y1,Z1)+(X2,Y2,Z2)。

表2 中步驟3 計算U2，步驟4 計算U1，步驟5的t4計算H，步驟8 計算S2，步驟9 計算S1，步驟11的t3計算r，依次代入后，X3可得

表2 普通點加運算步驟

步驟16 中的Y3可得

步驟17 中對Y3做0.5 倍模加，Z3代入即可，所得結果與式(6)一致。

在Comb 算法中，由于每次點加的αiG都是預先計算并存儲的，其點坐標格式均為(X,Y,1)，Z=1不需要存儲，且可以將式(6)做進一步優化和整理，并按照時序排布操作步驟如表4 所示，共需13 個周期可完成一次特殊點加運算(X3,Y3,Z3)=(X1,Y1,Z1)+(X2,Y2,1)。

表4 與表2 相似，依次代入后即可驗證是否與式(6)一致。點加運算結束后會立刻進行倍點運算，所以在2 種點加運算的操作步驟上都對最后輸出的坐標做了調整，以適應點加倍點的連續運算，最后一步的模加減均與下一運算的乘法同時運行，讓乘法器使用率達到最高。

Curve25519 的點加和倍點是采用階梯算法一起完成的，按RFC7748 中的參考算法，將操作步驟重新排布如表3 所示，共需12 個周期完成操作，可根據輸入的(X1,X2,Z2,X3,Z3)計算得到(X2,Z2,X3,Z3)。

表3 Ladder step 階梯運算步驟

表4 特殊點加運算步驟

這樣，根據曲線和算法來調用不同的點加控制和倍點控制，使用乘法器和模約簡以及其他必需的模運算單元來完成對應操作。由于表1～表4 的運算步驟僅僅是控制數據的選擇輸入和輸出，所有運算都由模運算完成，中間計算結果也存儲在寄存器堆，所以運算本身占用硬件資源較少。

3.3 模運算設計

模運算部分主要是模加減、模乘和模逆的設計。在點加倍點的操作周期中需要正常模加和模減單元各一個，還出現了特殊模運算單元，例如對于素數域中2 個數值M和N，需要計算3M、8M、1.5M?N、M?2N，以及模逆需要的0.5M，此類運算算法簡單，均設置獨立單元用于運算。此外，由于乘法消耗時間較長，為配合點加和倍點的周期操作，在模約簡的輸出端直接接入一個模加單元和一個模減單元，使模約簡的結果可以直接進行模加減操作，從而能在一個周期內安排更多的操作，整個模運算部分的結構如圖2 所示。

圖2 模運算單元結構

模乘采用模約簡結構，即先進行乘法運算，再對乘法結果做模約簡。由于secp256r1 的數據位寬為256 位，Curve25519 曲線的數據位寬為255位，故直接采用單元庫中已有的256 位乘法器，滿足2 條曲線的位寬需求，并在一個周期內完成乘法。模約簡需單獨設計，secp256r1 曲線的p=2256?2224+2192+296?1，乘法操作后得到一個512位的乘法結果A，易知A＜p2，用二進制表示后分為16 個數據段，設A=A152224+A142192+…+A1232+A0，其中Ai的寬度為32 位，則約簡公式如算法4所示[18]。

算法4secp256r1 的快速模約簡公式

輸入A=(A15||A14||A13||…||A2||A1||A0)

輸出B=Amodp

快速約簡通過多個數據的加減操作來實現，僅在最后一步取模運算。需要注意，由于S2的高位為0，因此最后得到的待約簡數據至多有5 個進位或5 個借位，比較常見的方式是采用高位為0 的形式化簡以上計算式，進行再次約簡，但即使再次模約簡，也會出現超過模p的情況，依然需要進行判斷并計算輸出。由于設計包含2 條曲線的模約簡，需要考慮單元間的復用關系。

對于Curve25519，p=2255?19=2255?24?22+1，位寬為255 位，一次乘法操作的結果A在二進制下位寬最高位510 位，根據p的形式特點，設A=A2542508+A2532506+A2522504+…+A122+A0，Ai的寬度為2 位，利用此形式的A對p做模約簡。

模約簡過程分為兩輪，第一輪處理高254 位，首先要對A的最高項A2542508做約簡，需利用A2542253乘以p值，得到多項式A2542508?A2542257?A2542255+A2542253，消去最高項后得到剩余項A2542257+A2542255?A2542253，以此類推，直至用A12821乘以p得到多項式A1282256?A12825?A12823+A12821，消去A1282256項后得到剩余項A2542257+A2542255?A2542253。第一輪運算得到的約簡結果分為兩部分，第一部分是A的低256 位，即A1272254+…+A122+A0，此部分未操作，直接保留設為為T；第二部分是所有剩余項的和，整理后為A2542257+(A253+A254)2255+(A252+A253?A254)2253+(A251+A252?A253)2251+…+(A128+A129?A130)25+(A128–A129)23+(?A128)21，對第二部分的高兩項進行第二輪約簡，消去A2542257+(A253+A254)2255，得到新增剩余項為A25426+(A253+2A254)24+A25322?(A253+A254)。至此，兩輪約簡整理完畢，將各位置的數值按規律重新整理為和項和差項后，可得算法5 所示的快速約簡算法。

算法5Curve25519 曲線快速約簡算法

輸入A=(A254||A253||A252||…||A2||A1||A0)

輸出B=Amodp

考慮到高位為0 的各部分，以及T所取的位數為256 位，最后一步待約簡數據會有至多5 個進位或2 個借位，與secp256r1 曲線的情況基本類似，故采用加減法陣列計算出2 種模約簡所有最后可能的輸出結果，通過判斷加減法的標志位來選擇輸出正確結果，其結構如圖3 所示。

圖3 模約簡結構

加減法陣列中需參數p～5p參與運算，需要在寄存器堆中存放2 條曲線的5p參數，其中，參數p在標量乘中多次使用，設置為參數形式存儲；2p在運算開始時通過p相加得到；3p通過3 倍模加單元計算得到；4p通過模減單元5p?p運算得到，這些值在首周期結束時固定在加減法陣列的輸入端，并在次周期開始參與運算，與乘法前后銜接。

模逆運算使用二進制右移算法，其基本算法如算法6 所示[20]。

算法6二進制右移算法求模逆

輸入素數p和a∈[1,p?1]

輸出a?1modp算，算法中每次向右移位一次，并需要做除2 操作，但對于二進制計算很容易實現。如需更快的速度，則可以考慮每次右移2 位，但需要額外做除4 操作并需要更多的判斷分支，故本文沒有采用?？紤]到使用標量乘進行橢圓曲線加密時也需要使用模逆運算，且與標量乘操作無時間上的沖突，所以將模逆模塊單獨引出，設立標志位以供標量乘內部使用或外部調用。

4 驗證與實現

完整的標量乘工作流程如圖4 所示。圖4 僅表示了重要的操作步驟，點加Z和點加1 分別使用表2和表3 中的點加流程，臨時寄存器在預計算和迭代計算部分都有使用，但圖4 中未全部標出，模約簡也未在圖中給出。運算時各環節所需要消耗的時鐘周期如表5 所示，其中包含了各工作狀態之間轉換消耗的周期數。

經綜合后所得面積速度與文獻[4-7,25-26]進行對比，得到表6 所示的結果，表6 中的AT是門數與單次運行時間的乘積。

表5 不同運算消耗周期情況

最終設計采用55 nmCMOS 工藝庫綜合，主頻最高可到625 MHz，此主頻下的門數為1 022×103個，按表5 中平均消耗周期計算，對于secp256r1曲線，計算基點的標量乘單次只需2.60 μs，每秒可計算38.5 萬次；計算普通點的標量乘需要6.56 μs，每秒可計算15.3 萬次；對于Curve25519，計算一次標量乘需要6.33 μs，每秒可計算15.8 萬次。

圖4 完整標量乘工作流程

表6 標量乘對比結果

本文設計消耗的門數是文獻[4]的2.3 倍，最高主頻提高了70 MHz，由于文獻[4]采用分段多次迭代的方式實現模乘，而本文設計采用了256 位的乘法器配合模約簡來實現模乘，頻率受乘法器限制，故主頻并未得到顯著提升，但配合3.2 節中的點加倍點操作，可以實現乘法和模約簡的二級流水，標量乘計算過程中每個周期都可以得到模乘結果，所以實際每秒標量乘的運算次數遠超文獻[4]。其余文獻也類似，差異在于模乘分段位數是32 位、64 位還是128 位，如文獻[5]采用的是128 位的分段模乘，每秒可計算標量乘8 萬次，是所有文獻中每秒運算次數最高的，但消耗的門數是本文設計的3.4 倍，同時本文設計的每秒運算量可達文獻[5]的1.9 倍。文獻[6-7,25-26]的門數消耗小于本文設計，但本文設計的主頻和門數均約為文獻[6]的3 倍，實際每秒運算次數遠超文獻[6]；文獻[7]的門數為本文設計的52%，但主頻只能到本文設計的30%，單次標量乘運行時間較短，但也是本文設計的18 倍；文獻[25]采用較低門數完成設計，且能達到相對較高的主頻，AT是所有文獻中最小的，但也是本文設計的3.8 倍；文獻[26]也僅在門數一項上低于本文設計，運算速度和AT與本文設計相比均無優勢。以上數據比較發生在普通點的計算上，本文設計在簽名時還能提供基點G的標量乘38.5 萬次，適用于服務器端的高速運算。

5 結束語

本文面向橢圓加密的運算需求設計了標量乘單元，采用快速模約簡的方式實現模乘運算，支持secp256r1 和Curve25519 曲線的標量乘，并能夠在輸入基點時得到更快的運算速度以應對橢圓曲線簽名和公鑰的生成需求，設計同時支持多標量乘以應對橢圓曲線驗簽使用，將模逆端口單獨引出以應對簽名和驗簽操作中的模逆需求，只需調用本單元即可完成加密算法中絕大部分操作。相較于其他設計，本文設計的運算速度有較大優勢，能夠提供secp256r1 曲線上每秒38.5 萬次的基點標量乘和每秒15.3 萬次的普通點標量乘，可以實現Curve25519 每秒15.8 萬次標量乘。