基于AVX512 的格密碼高速并行實現

雷斗威，何德彪，羅敏，彭聰

（武漢大學國家網絡安全學院，湖北 武漢 430072）

0 引言

隨著量子計算的迅速發展，現有的公鑰密碼體系受到了嚴重威脅。Shor 量子算法［1］利用量子計算機，可在多項式時間內破解大整數分解問題和橢圓曲線上的離散對數問題，從而攻破現有的公鑰密碼算法（RSA）和橢圓加密算法（ECC）。為了應對量子計算帶來的威脅，能夠抵抗量子攻擊的密碼即后量子密碼（PQC）成為國內外的研究熱點。

美國國家標準技術研究院（NIST）于2016 年開啟了后量子密碼標準的征集活動，旨在提出PQC 標準，以應對隨時可能到來的量子計算的威脅。格密碼因性能與安全性的優秀表現，在后量子密碼標準征集活動中占據了重要地位。經過多年的篩選和評估，NIST 于2022 年5 月確定 了4 個PQC 標準算 法，涵蓋公鑰加密（PKE）、密鑰封裝（KEM）和數字簽名，其中3 個算法是基于格的，可以看出格密碼在PQC中的重要地位。

Kyber 是NIST 確定的4 個PQC 標準算法之一，安全性基于模誤差學習（MLWE）問題。簡單而言，判定性MLWE 問題是指：給定上的隨機矩陣A以及上服從特定分布的s和e，區分(A，t=As+e)和隨機選取的(A，t)是困難的。搜索性MLWE 問題是指：給定隨機矩陣A以及t=As+e，其中s和e服從上的某些分布，求得s是困難的?，F有的量子算法并不能有效地解決MLWE 問題，從而保證了Kyber的抗量子安全性。

隨著PQC 標準的確定，為了后續的推廣和落地工作，對標準算法的高速實現的需求也日益增加。Kyber 包括一組PKE 算法以及利用Fujisaki-Okamoto轉換得到的KEM 算法。由于Fujisaki-Okamoto 轉換主要是哈希（Hash）操作，并非實現的重點，因此主要聚焦于對其PKE 算法的優化實現。本文基于512位高級向量擴展（AVX512），針對不同安全級別下的Kyber進行高速實現。根據Kyber 中q為12 bit 的特點，使用16 bit對其進行存儲與運算，將512 bit寄存器分為32組16 bit 進行并行處理，實現32 路并行多項式運算。同時，利用惰性模約減的思想大大減少模操作的次數，充分發揮了16 bit存儲空間的性能。在抽樣過程中，利用AVX512 每次同時產生8 組偽隨機比特串并對其進行合理分配，從而高效完成了對多項式的抽樣。

1 相關工作

為了更好地構建格密碼，REGEV［2］提出容錯學習（LWE）問題，使得格密碼方案的構造更加便捷。為了得到更優的計算效率，LANGLOIS 等［3］將LWE問題進行了擴展，得到了其模格形式MLWE，極大提高了格 密碼的 計算效 率。BOS 等［4］基 于MLWE 問題，構建了選擇明文不可區分（IND-CPA）安全的公鑰加密方案，并使 用Fujisaki-Okamoto 轉換［5］，得到了選擇密文不可區分（IND-CCA）下安全的密鑰封裝方案。SEILER［6］利用高級向量擴展2（AVX2）指令集，加速實現了快速數論變換（NTT），并將其運用到格密碼方案NewHope 上，獲得了較好的性能提升。LONGA 等［7］針對特殊的素數q，提出一種新的模約減技巧，提高了模約減效率，并利用該技巧加速NTT的實現。GUERON 等［8］利用冗余抽樣技術，提高了拒絕抽樣的成功率，降低了對偽隨機比特串的需求，同時結合AVX2 指令集，實現了高效的并行抽樣。ROY［9］利用AVX2，并行實現了不同安全級別下的格密碼算 法Saber。CABRAL 等［10］利 用AVX512 指 令集，以2 種不同的方式并行實現了安全哈希算法3（SHA-3）系列函數：第1 種方式是通過將函數內部并行化，提升了單個函數的執行效率；第2 種方式是同時執行8 組不同輸入的SHA-3 函數，從而可以同時輸 出8 組不同 的Hash 運算結果。ALKIM 等［11］在ARM Cortex-M4 處理器上實現了Kyber 和NewHope兩種格密鑰方案，填補了該方面實現的空白。ALKIM 等［12］對RISC-V 指令集進行擴展，并基于該擴展設計一種非常緊湊的PQC 軟硬件協同實現方式，獲得了較大的性能提升。WONG 等［13］設計一種平衡速度和資源的架構，并用該架構在現場可編程門陣列（FPGA）上實現了SHA-3 系列函數，實驗結果顯示：最終吞吐量達到較為理想的16.51 Gb/s。顧麗紅等［14］通過分析AES 加解密算法，結合龍芯平臺體系結構特征，提出基于多媒體指令擴展（SIMD）優化AES 性能的方法，獲得了較好的性能提升。KARMAKAR 等［15］研究了Saber 在資源受限的ARM Cortex-M 系列微控制器中的實現方式，并針對速度和內存資源的平衡進行優化，為未來物聯網設備中使用格密碼方案做好了前期準備工作。FRITZMANN 等［16］針對格密碼方案NewHope、Kyber和Saber 提出一系列硬件加速器，并將其整合進RISA-V 的流水線，同時增加了29 條新的RISA-V 指令集，并在FPGA 上進行了模擬，最終提升了格密碼方案的實現性能。ZHOU 等［17］針對Dilithium 方案提出一種基于FPGA 的軟硬件協同處理方式，提高了實現效率。郭渝洛等［18］提出一種基于SIMD 的并行傅里葉空間圖像相似度算法，經該算法優化后的程序可獲得平均5.132 倍的加速，并且具有較強的魯棒性。

2 相關知識

2.1 基本符號和定義

2.2 快速數論變換

對于環Rq=Zq[X]/(Xn+1)上的多項式乘法，如果直接進行相乘，則計算復雜度是O(n2)。在格密碼方案頻繁使用多項式相乘的情況下，效率低下。通過快速數論變化可以將多項式乘法的計算復雜度從O(n2)降到O(nlbn)，從而有效提高計算效率。

NTT 是有限域上的離散傅里葉變換，詳細內容可參見文獻［19］，下面對其核心思路進行概括描述。在Kyber 所選參數中，n取256，滿足n|q-1 且n=2d（d=8）的條件，因此有限域Zq上存在n原根，記作ζn，滿足條 件且那么環Rq=Zq[X]/(Xn+1) 可以寫作另一種形式Rq=通過利用環上的中國剩余定理［20］可以得 到，其中×表示直積。利用同構性質，可以將環上的多項式計算投射到環上。繼續對環使用環上的中國剩余定理并重復迭代該過程，最終可以得到同構

2.3 Kyber 和AVX512 指令集

Kyber 詳細算法可參見文獻［21］，不再進行重復說明。

2008 年，英特爾發布了高級矢量擴展（AVX）指令集，引入256 bit 位寬矢量指令。2013 年，英特爾將其位寬擴展至512 bit 并命名為AVX512 指令集。利用AVX512 指令集，計算機既可直接對512 bit 數據進行處理，如進行比特異或非等操作，又可將512 bit數據進行劃分，并行處理多組數據，如32 組16 bit 數據、16 組32 bit 數據、8 組64 bit 數據等。AVX512 指令集的誕生擴展了計算機處理位寬，提高了計算機并行效率。

為了能夠使編程人員不用進行匯編編程就可以直接編寫匯編語句，英特爾設計了內聯指令并將SIMD 指令集（包括AVX512 指令集）與內聯指令進行一一對應。采用內聯指令避免了繁瑣的匯編編程，增強了程序的易讀性和可移植性。在后續使用AVX512 時，本文將使用內聯指令描述算法，使表達過程更加通俗易懂。內聯指令_mm512_loadu_epi16可通過給定地址將512 bit 數據從內存裝入寄存器中，_mm512_storeu_epi16 可將寄存器內512 bit 的數據存到給定地址的內存中，在使用AVX512 指令集進行實現的過程中，會頻繁使用這兩條指令進行數據的裝載和存儲。在后續部分，將忽略重復性的數據裝載和存儲操作，并專注于算法本身，從而使對算法的描述更加簡潔。對于一些算法計算過程中所要用到的參數，利用_mm512_set 和_mm512_set1 系 列指令，通過立即數裝載進寄存器，這在后續部分將多次重復使用且不再贅述。

3 多項式計算優化

Kyber 上的計算包括多項式向量相加減、多項式向量相乘、多項式矩陣乘以多項式向量以及多項式向量內積。由于后兩者的計算可由前兩者得到，因此本文的重點在于前兩者的優化，即多項式向量相加減和多項式向量相乘的優化。

3.1 Zq上的計算優化

多項式Rq=Zq[X]/(Xn+1)上的計算在具體實現過程中最終取決于Zq上的模運算。為了提高多項式的計算效率，需要對Zq上的模運算進行優化。

3.1.1 惰性模約減

惰性模約減技術即不再每一次計算過程中進行模約減，而是等到數據快要溢出存儲空間時再進行統一處理。在Kyber 中，q取3 329，占12 個比特長，但本文采用16 bit 存儲Zq上的數據。因此，在計算過程中，如果數據未超過16 bit 而是僅超出絕對最小剩余系Sq，便可正常存儲與計算，將不再對其進行模約減操作。利用該方法，通過合理利用計算機存儲資源，減少了大量的模約減操作，提高了Zq上的計算效率。對于模加減而言，當結果未超出16 bit 時，直接進行加減忽略模約減操作，將這種直接加減法記作DAdd 和DSub。

3.1.2 優化的蒙哥馬利模約減

當Zq上的數據達到16 bit 長度時需要對其進行模約減操作。在進行Zq上的模乘運算時會得到一個無窮范數較大的相乘結果，因其比特長度已經超過16 bit，需要對其進行模約減操作。一般而言，對無窮范數較大的數據進行模約減操作需要除法運算，對于計算機而言，這將耗費大量的計算資源。1985 年，MONTGOMERY［22］提出一種快速的模約減方法，即蒙哥馬利模約減。蒙哥馬利模約減僅利用少量的加法和乘法操作，便可完成對較大數據的模約減操作，解決了一般模約減過程中除法效率過低的問題。傳統的蒙哥馬利模約減會在最后一步進行一次判斷操作，并依據判斷結果對數據進行修正，從而將數據徹底約減到Sq。仿照惰性模約減的思路，取消最后一次判斷以及修正操作，在數據不溢出16 bit 存儲空間的情況下，不會影響最終的計算結果，具體見算法1 所示，在后文中，將算法1 簡記作MontRed。需要注意到的是：對于最終得到的結果r1≡a·b·β-1(modq)中額外引入的β-1，將在后續統一處理。

3.1.3 模乘運算

由于相乘后的數據一般超過16 bit 的存儲范圍，因此需要對結果進行模約減，直接采用MontRed 進行模約減。在Kyber 的具體實現過程中，所用到的模乘可分為2 種：1）一般情況下的模乘，兩乘數均未知，這常用于多項式相乘中，本文將其記作MMul；2）一乘數提前確定的模乘，這常用于NTT 和INTT轉換的過程中，將其記作FMMul。對于MMul，對相乘后的結果直接使用MontRed 進行模約減，如算法2所示；對于FMMul，對提前確定的乘數進行預處理再進行后續操作，如算法3 所示。在后續部分，將采用合適的方法消除MontRed 引入的β-1。

3.2 Rq 上的計算優化

3.2.1 NTT 優化方案

在Kyber 中進行NTT 轉換的多項式系數范圍為(-q，q)。在算法4的第6行中，DAdd 和DSub不會進行模操作，因此每一輪迭代均會使數據的范圍發生改變。在經過7 輪迭代后，多項式系數的范圍從(-q，q)增加到(-8q，8q)。由于(-8q，8q)內的數據仍能被16 bit 正常存儲，因此不對其進行額外處理。

3.2.2 多項式乘法優化方案

通過NTT將n-1 階多項式上的乘法轉換為個一階多項式的乘法。不妨設(a1x+a0)，(b1x+b0)∈Zq[X]/(X2-σi)，相乘的過程可以寫作(a1x+a0)·(b1x+b0)=(a1·b1)x2+(a0·b1+a1·b0)x+a0·b0=(a0·b1+a1·b0)x+(a1·b1·σi+a0·b0)。對于a0·b1、a1·b0、a1·b1、a0·b0，使用MMul進行乘操作；對于(a1·b1) ·σi，由于σi是提前確定的乘數，預計算σi·βmod±q及σi·β·q-1mod±β，并使用FMMul 進行乘操作。對于加法，使用DAdd。由于省略了部分模約減操作，最終得到的結果范圍為(-2q，2q)。由于使用MontRed，最終計算的結果乘上了額外的β-1。這些將在INTT中統一處理。

3.2.3 INTT 優化方案

在INTT 的過程中，需要處理由MontRed 引入的β-1，同時需要防止計算過程中數據溢出16 bit 長度。在INTT 的每一輪迭代過程中系數的無窮范數均會翻倍，而16 bit 最多可以存儲的數據范圍為(-9q，9q)。因此，在數據即將產生溢出前，使用MontRed 對其進行模約減。設輸入數據的系數范圍為(-4q，4q)，整個INTT 迭代過程的數據變化如下：

第2 輪迭代：(-q，q)→(-2q，2q)；

第3 輪迭代：(-2q，2q)→(-4q，4q)；

第5 輪迭代：(-q，q)→(-2q，2q)；

第6 輪迭代：(-2q，2q)→(-4q，4q)；

在第1 和4 輪迭代中使用MontRed 進行模約減，同時不會將預乘上β。這樣在每一次模約減后，數據范圍將變成(-q，q)且會引入β-1，如算法5 的第11 和12 行所示。在INTT 的最后一輪迭代過程中，通過乘上2-(m-1)來消除每一輪迭代過程中引入的倍數2，通過乘上β-2來消除第1 和4 輪迭代中進行模約減引入的β-2，通過乘上β來消除最后一輪迭代過程中使用FMMul 進行模乘引入的β-1，如算法5 的第21 和22 行所示。利用改進的INTT 算法，消除了之前運算引入的β-1并得到了無窮范數較小的結果。

4 基于AVX512 的并行實現

本節詳細描述如何利用AVX512 進行具體實現。采用16 bit 存儲Zq上的元素，每512 bit 可存儲32 組16 bit 數 據。

4.1 Rq 上計算的并行實現

4.1.1 32 路并行加減

對于加減法而言，使用惰性模約減技術，利用_mm512_add_epi16 和_mm512_sub_epi16 實 現32 路并行加減法。

4.1.2 32 路并行模約減

4.1.3 32 路并行模乘

MMul 的實現思路是將兩輸入乘數相乘后進行MontRed 處理。本文同樣利用_mm512_mulhi_epi16和_mm512_mullo_epi16，32 路并行獲得兩相乘結果的高16 bit 和低16 bit 并將得到的結果輸入32 路并行MontRed 進行處理，如算法7 所示。FFMul 的32 路并行計算方法和上述類似，但需要注意的是由于提前確定了1/2 的乘數并存入了寄存器y中，可以預計算得到yqinv=y·q-1mod±β，因此，與MMul 相比，FFMul 少了一 次_mm512_mullo_epi16 運 算，提高了計算效率，如算法8 所示。

4.1.4 32 路并行NTT 和INTT

將第4.1.1～4.1.3 節中的32 路并行計算方法運用到NTT 和INTT 過程中，即可得到32 路并行NTT 和INTT，但有2 點需要注意：1）將預計算的參數進行冗余存儲使其達到512 bit，通過該方法，當向512 bit 寄存器裝載數據時可以直接讀入512 bit 冗余存儲的數據，而不再需要通過復雜的廣播指令進行數據的擴展；2）在NTT 和INTT 的第1 次迭代過程中，將多項式的系數裝載進寄存器，256 個系數共需要8 個寄存器，這對于AVX512 來說是可實現的。在所有迭代完成后，再將8 個寄存器中的數據存入內存，通過該方法，省去了大量裝載數據與存儲數據的操作，從而提升了計算效率。

4.1.5 32 路并行多項式乘法

利用NTT 將n-1 階多項式上的乘法轉換為個一階多項式的乘法。因為輸入的多項式是不確定 的，所以采 用MMul 進行相 乘。將32 組a1、a0、b1、b0分別裝載 到4 個不同的512 bit 寄 存器，再使用相應的32 路并行MMul 和32 路并行加法完成32 組一階多項式的相乘，從而實現32 路并行多項式乘法。

4.2 8 路并行SHAKE 函數

SHA-3 是NIST 于2015 年8 月5 日 發布的安全 哈希算法系列標準。SHAKE128、SHAKE256 是SHA-3中的兩個函數［23-24］，最大的特點是可以輸出任何比特長度。Kyber 通過SHAKE128 和SHAKE256 擴展隨機種子，獲得所需長度的偽隨機比特串，用于執行后續的抽樣操作。使用AVX512 可以同時執行8 路SHAKE128（SHAKE256）操作［10］，從而同時拋出8 組偽隨機比特串，具體如圖1 所示。

圖1 8 路并行SHAKE 函數Fig.1 8-way parallel SHAKE functions

4.3 并行抽樣

在Kyber 中利用拒絕抽樣技術［25］和SHAKE128，通過給定的種子生成上的隨機矩陣，具體如算法9 所示。由于是上的隨機矩陣，因此可以將其作為已經經過NTT 運算的矩陣使用，從而節省NTT 運算所耗費時間。

4.3.1 并行拒絕抽樣

對于隨機的12 bit，當數值在[0，q)時將其取出，否則將不使用該12 bit，繼續使用下一串隨機的12 bit 進行上述處理。根據條件概率，當成功抽取一個數值時，滿足[0，q)的均勻分布。這種抽樣方式被稱為拒絕抽樣［25］。設是第i行、第j列的元素的每一 個系數 范圍均 為[0，q)。為了使滿足上的均勻分布，需要使的每一個系數滿足[0，q)上的均勻分布。在具體抽樣過程中，將重復進行拒絕抽樣，直至抽出整個。

一個512 bit 寄存器可以裝載64 Byte 數據，每3 Byte 可以組成2 組12 bit 數據。利用_mm512_permutexvar_epi8，將3 Byte的中間部分重復一次，接著使用_mm512_mask_srli_epi16、_mm512_and_si512以及預先計算好的mask 消除多余部分，如圖2 所示。利用上述方法可以將寄存器中前48 Byte 數據同時轉換為32 組12 bit 數據存放在一個512 bit 寄存器中，其余數據在下一次抽樣中使用。使用_mm512_cmplt_epi8_mask 判斷哪些數據在[0，q)內并生成標志位。利用標志位，將滿足要求的數據作為的系數存入內存。

圖2 數據排列與處理Fig.2 Arrangement and processing of data

4.3.2 偽隨機比特串分配

由于直接生成隨機比特串較為復雜，Kyber 中利用SHAKE128 擴展輸入的種子得到偽隨機比特串，并利用該偽隨機比特串替換隨機比特串完成對A的抽樣，從而提高效率且不會降低安全性。如4.2 節所述，通過AVX512 能夠進行8 路SHAK128 操作，同時生成8 組偽隨機比特串，充分利用該8 組偽隨機比特串對抽樣性能的提升非常重要。本文的核心方法是通過對多項式進行拆分抽樣，使其能夠完全利用8 組偽隨機比特串。

當k取2 時中共有4 個多項式。對于中第i行、第j列的元素αi，j，利 用SHAKE128(seed‖ ‖i‖j0)生成其低128 bit 系數，SHAKE128(seed‖i‖j‖1)生成其高128 bit 系數。通過該方法，生成4 個多項式需要用到8 組偽隨機比特串，從而充分利用了8 路并行SHAKE128 函數且降低了SHAK128 的循環次數。

當k取3 時，中共有9 個多項式。對于中第3 行、3 列的元素，將其拆分成8 個部分，每個部分有32 項系數。利 用SHAKE128(seed‖ 3 ‖3 ‖index)為每一部分生成偽隨機比特串。對于的其他元素，直接利用SHAKE128(seed‖i‖j)生成對應的偽隨機比特串，不再進行拆分。通過該方法，充分利用了8 路并行SHAKE128 函 數。

當k取4 時，中共含16 個多項 式，正好是8 的倍 數。因 此，對 于中 第i行、第j列的元 素，直 接利用SHAKE128(seed‖i‖j)生成對應的偽隨機比特串，不再進行拆分，從而充分利用了8 路并行SHAKE128 函數。

4.4 并行抽樣Sη

在Kyber 中，利用中心二項分布（CBD）抽樣抽取Sη，如算法10 所示。CBD 抽樣利用固定長度的偽隨機比特串抽取Sη上的元素，未使用較為繁瑣的拒絕抽樣技術，運行效率更高。

由于直接生成隨機比特串較為復雜，Kyber 中利用SHAKE256 擴展輸入的種子得到偽隨機比特串，利用偽隨機比特串替換隨機比特串完成對Sη的抽樣，η取值為2 或者3。當η為2 時，每個長度為4 的偽隨機比特串可通過CBD 抽樣得到Sη的一個系數，整個Sη的CBD 抽樣需要消耗長度為1 024 bit 的偽隨機比特串，而SHAKE256 一次循環能產生長度為1 088 的偽隨機比特串，因此只需一次Squeeze 操作，便可完成一次Sη的CBD 抽 樣。當η為3 時，每個長度為6 的偽隨機比特串可通過CBD 抽樣得到Sη的一個系數，整個Sη的CBD 抽樣需要消耗長度為1 536 bit 的偽隨機比特串，需要SHAKE256 循環2 次。當進行實現時，較難處理3 bit 的并行操作，因此利用冗余比特技術，對每4 bit 進行并行處理，并忽略最高比特。通過上述方法，整個Sη的CBD 抽樣需要消耗長度為2 048 bit 的偽隨機比特串，SHAKE256依然只需循環2 次，未造成額外的開銷，如圖3所示。

圖3 CBD 抽樣的比特處理Fig.3 Bit processing of CBD sampling

5 性能分析

利用AVX512 實現了不同安全參數下的Kyber算法，即Kyber512、Kyber768 和Kyber1 024，并與其C 語言版本的實現進行了性能對比與分析［21］。使用的處理器是11th Gen Intel?CoreTMi5-1135G7@2.40 GHz 2.42 GHz，內存大小為16 GB。操作系統選用Windows 10，編譯器選用GCC，優化等級為-O3。將各算法運行20 000 次，并計算其平均時鐘周期作為性能指標。

通過表1～表3 可以看出：在多項式計算方面，由于對算法進行優化且利用AVX512 實現了32 路并行，從而獲得了55～95 倍的加速；在多項式抽樣方面，利用8 路SHAKE 函數以及優化的抽樣方法獲得了4～7 倍的加速；從整體來看，對于包含抽樣操作的密鑰生成算法和加密算法獲得了10～16 倍的加速，對于主要由多項式計算構成的解密算法獲得了約56 倍的加速。

表1 Kyber512 的性能測試與對比 Table 1 Performance testing and comparison of Kyber512

表2 Kyber768 的性能測試與對比 Table 2 Performance testing and comparison of Kyber768

表3 Kyber1 024 的性能測試與對比 Table 3 Performance testing and comparison of Kyber1 024

6 結束語

本文針對Kyber 算法的公鑰加密部分進行高速并行實現。在多項式計算方面，利用惰性模約減、優化的蒙哥馬利模約減、一般模乘和固定乘數模乘等技術，提高了模運算的效率，采用NTT 和INTT 變換提高了多項式計算的效率，將優化后的模運算融入NTT 和INTT 過程，同時利用AVX512 進行高速并行實現，獲得了55～95 倍的加速。在多項式抽樣方面，通過拆分抽樣、冗余抽樣等技術優化抽樣過程，同時利用AVX512 進行高速并行實現，獲得了4～7 倍的加速。在整體上，獲得了10～56 倍的性能提升。在后續工作中將繼續研究不同后量子密碼標準在不同硬件平臺上的高速實現。