Data Center Heterogeneous Acceleration Software-Hardware System-Level Platform Based on Reconfigurable Architecture
-
摘要:
构建数据中心加速服务的软硬件系统级原型平台,需要考虑高计算能力、扩展性、灵活性和低成本等因素. 为了提高数据中心的能力,从软硬件协同的角度研究数据中心异构计算在云平台架构、硬件实现、高速互连和应用等方面的创新,研究设计并构建了一个可重构组合的软硬件加速原型系统,简化了现有以处理器为中心的系统级计算平台构建方法,实现目标软硬件设计的快速部署与系统级原型验证. 针对以上目标,通过解耦的可重构架构设备虚拟化和远程映射等方法,发掘独立计算单元的潜力,构建了一套ISOF(independent system of FPGA(field programmable gate arrays))软硬件计算平台系统,可使其超越普通服务器设计所能提供的能力,实现计算单元低成本高效扩展,使客户端可灵活使用外设资源,并且为满足系统级通信挑战,设计了一套计算单元之间的通信硬件平台和交互机制. 此外,为提升软硬件系统级平台的敏捷性,ISOF提供了灵活统一的调用接口. 最后,通过对平台目标系统级的分析评估,验证了该平台在满足了当下计算与加速需求下,保证了高速、低延时的通信,以及良好的吞吐率和弹性扩容效率,另外在高速通信的基础上改进的拥塞避免和丢包恢复机制,满足了数据中心规模通信的稳定性需求.
Abstract:Constructing a software and hardware system-level prototype platform for accelerating data center services requires the consideration of factors such as high computing power, scalability, flexibility, and low cost. To enhance data center capabilities, research from the perspective of software-hardware synergy has been conducted on the innovation of heterogeneous computing in cloud platform architecture, hardware implementation, high-speed interconnection, and applications. A reconfigurable and combinable software-hardware acceleration prototype system is designed and built to simplify existing processor-centric system-level computing platform construction methods, enabling rapid deployment and system-level prototype validation of target software-hardware designs. To achieve these objectives, methods such as decoupled reconfigurable architecture device virtualization and remote mapping are utilized to uncover the potential of independent computing units. An ISOF (independent system of FPGA) software-hardware computing platform system is constructed to surpass the capabilities of conventional server designs, enabling low-cost and efficient expansion of computing units while allowing clients to flexibly utilize peripheral resources. To address system-level communication challenges, a communication hardware platform and interaction mechanism between computing units are designed. Additionally, to enhance the agility of the software-hardware system-level platform, ISOF provides a flexible and unified invocation interface. Finally, through the analysis and evaluation of the system-level objectives of the platform, it has been verified that the platform meets the current computing and acceleration requirements, ensuring high-speed, low-latency communication, as well as good throughput and efficient elastic scalability. In addition, improvements have been made in congestion avoidance and packet recovery mechanisms based on high-speed communication, meeting the stability requirements of communication at data center scale.
-
传统的公钥密码大多基于RSA算法[1]以及椭圆曲线密码[2]. Shor算法的提出证明大数分解难题可用量子计算机在多项式时间内解决. 随着量子计算机的不断发展,传统公钥密码面临被攻破的威胁. 2016年,美国国家标准技术研究院(National Institute of Standards and Technology, NIST)在全球范围内启动了针对后量子密码算法的征集[3]. 后量子密码的标准化和迁移是下一代密码技术的重要方向. 在众多的后量子密码算法中,基于格的密码算法由于在安全性、公私钥尺寸、计算速度上达到了更好的平衡,成为最有前景的后量子密码算法之一. 到2022年7月,NIST初选了4个后量子标准算法,其中3个是签名算法,1个是加密算法.4个标准算法中3个都是基于格的密码算法,说明目前格基密码是后量子密码主流技术路线. 伴随着美国NIST后量子密码标准的推出,如今学术与工业界都在加快推进网络安全生态向后量子时代迁移,高效安全的后量子密码实现对后量子迁移与应用尤为重要.
基于格的密钥封装机制主要使用:一般格、理想格和模格以及NTRU格. 目前NIST拟标准化的密钥封装算法Kyber[4]是基于LWE技术路线. 但基于LWE技术路线的后量子密码标准面临较为复杂的专利因素. 相较于LWE技术,NTRU密码体制更为成熟且无专利风险. 传统的基于NTRU格的密钥封装机制通常需要采用相较于基于{R, M}LWE密钥封装机制更大的模数和更小的密钥空间来满足方案正确性[5-6]. 我国学者基于E8格编码[7]和NTRU格提出了CTRU和CNTR密钥封装算法[8],简记为CTRU,其中,CNTR是CTRU的一个简化版本. CTRU和CNTR本质上还是NTRU机制,其密钥生成和解密过程与传统的NTRU机制保持一致. 其优点是采用了尺度化E8格编码,提高纠错能力,能够容忍更大的秘密范围和更小的模数,并支持密文压缩,从而在安全性和带宽上相比传统NTRU机制更有优势. 目前,CTRU算法的标准化工作已经在我国密码标准委员会正式立项[9],其中每个算法提供了3个方案CTRU-512,CTRU-768,CTRU-1024,分别对应不同的安全级别. 因此,对CTRU算法的高效实现研究具有重要的应用价值. 但是,在之前的CTRU密钥封装方案工作中,只提供了CTRU-768方案的C参考实现和AVX2实现,实现的不完整且有较大的性能优化空间,这为CTRU算法在实际中的应用和实验带来很大的不便,迫切需要解决.
典型的格密码系统的核心运算是矩阵算术运算和多项式环算术运算. 多项式乘法是大部分格密码中较为复杂且耗时的运算,一般采用快速数论变换(number theoretic transform, NTT)进行加速. 当前格密码的软件优化实现一般通过SIMD指令实现运算的并行加速,如Intel处理器的AVX/AVX2/AVX512指令集[10-12],ARM处理器上的DSP,NEON指令集[13]. AVX是x86(-64)处理器上的扩展指令集,增加了256b向量寄存器,以及操作向量寄存器的相应指令(AVX2增加了更多指令). AVX2指令集包含在向量寄存器上同时操作多个数据的指令,这种一条指令同时操作多个数据的概念被称为单指令多数据(single instruction multiple data, SIMD). 在众多提交至NIST后量子标准遴选的格密码方案实现中,除了纯C参考实现,AVX2是必要实现之一. 因此对于格密码方案的软件工程实现来说,C参考实现和AVX2向量化并行实现是一个重要的方向.
本文的主要研究目标在于针对CTRU方案提出Intel平台上高效并行的软件实现,不仅对目前已有实现进行较为系统的优化,除此之外还完成了之前没有实现的CTRU-512和CTRU-1024方案的C参考实现和AVX2优化实现,并创造CTRU方案在Intel CPU平台上新的性能记录. 利用AVX2指令的SIMD特性,对CTRU底层较为耗时的多项式乘法、多项式求逆、模约减等运算进行较为系统的优化实现,最大化提升CTRU在Intel平台的性能. 具体而言,本文主要贡献如下:
1)提供了首个CTRU-512和CTRU-1024方案的参考实现. 现有CTRU-768方案的正向NTT实现采用的是6层基-2 NTT结合1层基-3 NTT的混合基NTT.而CTRU-512和CTRU-1024方案不需要进行基3NTT,因此目前CTRU-768方案的代码无法直接扩展到CTRU-512和CTRU-1024方案实现上,而CTRU-512和CTRU-1024方案的正向NTT采用的都是1层分解加6层基-2NTT,因此CTRU-512和CTRU-1024方案的正向NTT参考实现可以使用相同的代码模板. 除此之外,本文在之前延迟约减的基础上进行了进一步优化,利用中心Barrett约减对逆向NTT做出了更精细的针对系数索引的延迟约减.
2)对于CTRU-1024方案的多项式求逆实现,由于多项式矩阵维度大,现有的CTRU-768方案的多项式实现无法直接推广到CTRU-1024方案实现. 本文提出利用Bernstein快速多项式求逆算法实现CTRU-1024方案的Zq[x]/(x8−ζr(i))多项式环求逆,使用Karatsuba技术加速小度数多项式乘法.
3)针对CTRU方案的主要性能瓶颈完成AVX2优化实现:针对正向NTT和逆向NTT提出了层融合实现技巧,进一步降低存取指令,提升多项式乘法速度;利用AVX2并行指令优化耗时的多项式求逆运算实现,对多项式矩阵求逆方法采用循环展开技巧实现模幂运算. 对Bernstein快速多项式求逆方法用AVX2进行向量化实现加速,针对多项式序列化和反序列化进行AVX2优化实现. 之前已有的实现没有对性能瓶颈SHA-3哈希函数进行AVX2优化,本文使用AVX2向量化指令对SHA-3哈希函数进行并行汇编实现.
4)本文提出的CTRU 3个方案的AVX2优化实现相较于本文完成的参考实现在密钥生成,密钥封装和密钥解封装上分别可以达到56%~91%,74%~90%, 70%~83%的性能提速,对于CTRU-768方案,相较于目前已有AVX2实现,提升了8%~11%,是目前CTRU方案在Intel CPU平台上实现的最好性能记录. 这为CTRU方案的实际应用和实验打下重要基础,具有重要的应用价值. 相较于已经标准化的Kyber密钥封装方案的AVX2优化实现,本文的CTRU AVX2优化也在性能上有一定优势.
1. 相关工作
目前有很多针对NIST PQC候选算法的优化实现,其中大部分优化实现主要关注FPGA[14-16]、x86(主要是AVX2指令集优化实现)[17-19] 以及ARM Cortex-M4的实现[20-22],这些是NIST PQC竞赛时推荐的官方实现平台. Chung等人[23]提出了在NTT非友好环上使用NTT加速多项式乘法,他们采用此技术在NIST主要软件优化目标AVX2和Cortex-M4进行实现,创造了Saber方案和NTRU方案的最新性能记录. Hwang等人[24]对NTRU Prime的ntrulpr761/sntrup761这2个参数集提出了向量友好型的多项式乘法实现,与最新的AVX2实现相比达到了一定的提速. Alkim等人[25] 提出了2种针对NTT非友好环的不同NTT技巧,对NTRU Prime方案在ARM Cortex-M4平台上带来了9%~16%的性能提升. Dai等人[26]给出了定制的乘法算法结构以及对应的英特尔AVX2指令,并提出了依赖于参数的优化. Becker等人[27]提出了一种存储高效的Toom-Cook/Karatsuba多项式乘法技巧实现,并且带来了性能上的提升,他们在Cortex-M4上的实现能带来3~5倍的性能提速. Bernstein等人[28]提出在TLS 1.3中对sntrup761实现更高的性能提升,提出了批量密钥生成的AVX2并行优化实现. 在CTRU方案实现中,主要存在的难点是快速数论变换、多项式求逆、和小度数多项式乘法. 对于不同参数,运用的多项式求逆方法不一致. 并且目前CTRU-768方案中采用的多项式矩阵方法的多项式求逆实现由于矩阵维数大,无法在CTRU-1024方案的多项式环上直接实施. 因此本文需要解决如何加速CTRU 3个方案中中均存在的大量的小度数乘法,以及考虑如何实现CTRU-1024环上的多项式求逆,并且进一步优化目前已有的快速数论变换实现.
2. 预备知识
本节给出了本文工作相关的预备知识. 首先定义了变量符号、数学运算与格上困难问题,接着详细描述CTRU方案,并介绍了实现平台与优化方式.
2.1 基本定义及符号说明
令n和q为正整数. 定义集合Zq=Z/qZ={0,1,…,q−1},Rq=R/q. 定义多项式环R=Z[x]/(xn−xn/2+1),Rq=Zq[x]/(xn−xn/2+1),其中每个多项式系数均在环Zq中. 令f为Rq中的多项式,可表示为f=n−1∑i=0fixi. 多项式f的向量形式是 \boldsymbol{f}=({f}_{0},{f}_{1},…,{f}_{n-1}),其中 {f}_{i}\in {\mathbb{Z}}_{q}, i=\mathrm{0,1},…,n-1 . 对于多项式模数 q,{r}'=r{\mathrm{m}\mathrm{o}\mathrm{d}}^{\pm }q 表示 r' 是 r 约减到 [-\dfrac{q}{2},\dfrac{q}{2}) 的结果. {r}'=r{\mathrm{m}\mathrm{o}\mathrm{d}}^{+}q 表示将 r 约减到 [0,q) . 对于一个元素 a\in {\mathbb{Z}}_{q},{\left|\left|a\right|\right|}_{\infty } 表示 \left|a{\mathrm{m}\mathrm{o}\mathrm{d}}^{\pm }q\right| . 对于一个向量 \boldsymbol{b}=({b}_{0},…,{b}_{n-1}), 定义 {\left\|\boldsymbol{b}\right\|}_{\infty }=\underset{i}{\mathrm{m}\mathrm{a}\mathrm{x}}\left\|{b}_{i}\right\|_{\infty } . 对于集合 D,s\leftarrow D 表示 D 中随机选取一个元素 s . 中心二项分布 {B}_{\eta } 与正整数 \eta 有关,中心二项分布定义如下:采样 \left({c}_{1},{c}_{2},… ,{c}_{\eta },{d}_{1},{d}_{2},… ,{d}_{\eta }\right)\leftarrow {\left\{\mathrm{0,1}\right\}}^{2\eta }, 输出 {\displaystyle\sum\limits^{\eta }_{i=1}}({c}_{i}-{d}_{i}) ,采样多项式 f\leftarrow {B}_{\eta } 表示按照分布 {B}_{\eta } 采样多项式 f 中的每个系数.
2.2 CTRU方案介绍
CTRU是基于NTRU格的密钥封装方案,CTRU方案中包括一个IND-CPA安全的公钥加密和一个IND-CCA安全的密钥封装方案. CTRU公钥加密的密钥生成算法、加密算法和解密算法如算法1,2,3所示.CTRU选择的模数 q=3457 . CTRU选取的多项式环 {\mathcal{R}}_{q}={\mathbb{Z}}_{q}\left[x\right]/({x}^{n}-{x}^{\tfrac{n}{2}}+1) . {\varPsi }_{1} 和 {\varPsi }_{2} 都是 \mathcal{R} 上的分布. 密钥生成、加密、解密算法中的多项式采样都是通过中心二项分布采样算法生成的服从 {\varPsi }_{1} 或 {\varPsi }_{2} 分布的多项式. 函数 PolyEncode 和 PolyDecode 是用于对明文进行 {\mathrm{E}}_{8} 格编码和解码的函数,具体细节可见文献[8].
算法1. CTRU.PKE.KeyGen\left(\right) 密钥生成算法.
输出:公钥 pk:=h 、私钥 sk:=f .
① {f}',g\leftarrow {\varPsi }_{1} ;
② f:={pf}'+1 ;
③ 如果 f 在 {\mathcal{R}}_{q} 中不可逆,重启密钥生成;
④ h:=g/f ;
⑤ \mathrm{r}\mathrm{e}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{n}\left(pk:=h,sk:=f\right).
算法2. \boldsymbol{C}\boldsymbol{T}\boldsymbol{R}\boldsymbol{U}. \boldsymbol{P}\boldsymbol{K}\boldsymbol{E}. \boldsymbol{E}\boldsymbol{n}\boldsymbol{c}(\boldsymbol{p}\boldsymbol{k},\boldsymbol{m};\boldsymbol{c}\boldsymbol{o}\boldsymbol{i}\boldsymbol{n}) 加密算法.
输入:公钥CTRU-512和CTRU-768这2个方案、待加密消息 \boldsymbol{m}\in \boldsymbol{M} 、随机数 \boldsymbol{c}\boldsymbol{o}\boldsymbol{i}\boldsymbol{n} ;
输出:密文 \boldsymbol{c}.
① r,e\leftarrow {\varPsi }_{2} ;
② \sigma :=hr+e ;
③ c:=\left\lfloor \left. \dfrac{{q}_{2}}{q}(\sigma +PolyEncode(m\left)\right)\right\rceil \right. \mathrm{m}\mathrm{o}\mathrm{d}\;{q}_{2} ;
④ \mathrm{r}\mathrm{e}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{n}c.
算法3. \text{CTRU}.PKE.Dec(sk,c) 解密算法.
输入: 私钥 sk:=f 、密文 \boldsymbol{c} ;
输出:解密消息 \boldsymbol{m}.
① \text{m}:=PolyDecode\left(cf{\mathrm{m}\mathrm{o}\mathrm{d}}^{\pm }{q}_{2}\right) ;
② \mathrm{r}\mathrm{e}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{n}\;m .
通过前缀哈希FO转换\mathrm{FO}_{ID\left(p{k}\right), m}^{/ \perp} 将CPA安全的公钥加密CTRU. PKE转换为CCA安全的密钥封装. 前缀哈希FO转换 \mathrm{F}_{ID\left(p{k}\right), m}^{/ \perp} 是FO转换[29-30]的变体[31]. CTRU的密钥封装方案包含密钥生成、密钥封装以及密钥解封装,3个算法如算法4,5,6所示. 其中所有的哈希函数都是用SHA-3实例化的. 采样多项式所需要的随机字节流通过SHAKE-256生成. 密钥封装中的哈希函数 \mathcal{H} 用SHA3-512实例化.
算法4. \boldsymbol{C}\boldsymbol{T}\boldsymbol{R}\boldsymbol{U}. \boldsymbol{K}\boldsymbol{E}\boldsymbol{M}. \boldsymbol{K}\boldsymbol{e}\boldsymbol{y}\boldsymbol{G}\boldsymbol{e}\boldsymbol{n} 密钥生成算法.
输出: 公私钥对 \left(\boldsymbol{p}\boldsymbol{k},\boldsymbol{s}\boldsymbol{k}\right).
① \left(\boldsymbol{p}\boldsymbol{k},\boldsymbol{s}\boldsymbol{k}\right)\leftarrow \boldsymbol{C}\boldsymbol{T}\boldsymbol{R}\boldsymbol{U}. \boldsymbol{P}\boldsymbol{K}\boldsymbol{E}. \boldsymbol{K}\boldsymbol{e}\boldsymbol{y}\boldsymbol{G}\boldsymbol{e}\boldsymbol{n}\left(\right) ;
② z\stackrel{\mathrm{\$}}{\leftarrow }{\left\{\mathrm{0,1}\right\}}^{l} ;
③ \mathrm{r}\mathrm{e}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{n}\left(p{k}':=pk,s{k}':=\left(sk,z\right)\right).
算法5. \text{CTRU}.KEM.Encaps 密钥封装算法.
输入: 公钥 pk ;
输出: 协商密钥 K 、封装密钥 c.
① m\stackrel{\mathrm{\$}}{\leftarrow }M ;
② \left(K,coin\right):=\mathcal{H}\left(ID\left(pk\right),m\right) ;
③ c:=CTRU.PKE.Enc\left(pk,m;coin\right) ;
④ \text{return}\left(c,K\right).
CNTR是CTRU方案的变体,算法结构和CTRU基本类似,区别在于CNTR算法在加密过程中去掉了噪音多项式 \boldsymbol{e} . CTRU和CNTR方案的具体参数如表1和表2所示,多项式环维度 \boldsymbol{n} 有512,768,1 024这3种选择,分别对应NIST推荐的3个安全级别:安全级别Ⅰ,Ⅲ和安全级别Ⅴ. \left|\boldsymbol{p}\boldsymbol{k}\right|,\left|\boldsymbol{c}\boldsymbol{t}\right| 为公钥和密文的长度(单位为B).
表 1 CTRU推荐参数集合Table 1. CTRU Recommended Parameter Sets方案 n q {q}_{2} ({\varPsi }_{1},{\varPsi }_{2}) \left|pk\right| \left|ct\right| CTRU-512 512 3 457 {2}^{10} ({B}_{3},{B}_{3}) 768 640 CTRU-768 768 3 457 {2}^{10} ({B}_{2},{B}_{2}) 1 152 960 CTRU- 1024 1 024 3 457 {2}^{11} ({B}_{2},{B}_{2}) 1 536 1 408 表 2 CNTR推荐参数集合Table 2. CNTR Recommended Parameter Sets方案 n q {q}_{2} ({\varPsi }_{1},{\varPsi }_{2}) \left|pk\right| \left|ct\right| CNTR-512 512 3 457 {2}^{10} ({B}_{5},{B}_{5}) 768 640 CNTR-768 768 3 457 {2}^{10} ({B}_{3},{B}_{3}) 1 152 960 CNTR- 1024 1 024 3 457 {2}^{10} ({B}_{2},{B}_{2}) 1 536 1 280 算法6. \boldsymbol{C}\boldsymbol{T}\boldsymbol{R}\boldsymbol{U}. \boldsymbol{K}\boldsymbol{E}\boldsymbol{M}. \boldsymbol{D}\boldsymbol{e}\boldsymbol{c}\boldsymbol{a}\boldsymbol{p}\boldsymbol{s} 解封装算法.
输入:私钥 \boldsymbol{s}{\boldsymbol{k}}':=\left(\boldsymbol{s}\boldsymbol{k},\boldsymbol{z}\right) 、封装密钥 \boldsymbol{c};
输出:封装密钥 {\boldsymbol{K}}' 或 \widetilde{\boldsymbol{K}} .
① {\boldsymbol{m}}':=\boldsymbol{C}\boldsymbol{T}\boldsymbol{R}\boldsymbol{U}. \boldsymbol{P}\boldsymbol{K}\boldsymbol{E}. \boldsymbol{D}\boldsymbol{e}\boldsymbol{c}\left(\boldsymbol{s}\boldsymbol{k},\boldsymbol{c}\right) ;
② \left({\boldsymbol{K}}',\boldsymbol{c}\boldsymbol{o}\boldsymbol{i}{\boldsymbol{n}}'\right):=\mathcal{H}\left(\boldsymbol{I}\boldsymbol{D}\left(\boldsymbol{p}\boldsymbol{k}\right),{\boldsymbol{m}}'\right) ;
③ \stackrel{~}{\boldsymbol{K}}:={\mathcal{H}}_{1}\left(\boldsymbol{I}\boldsymbol{D}\left(\boldsymbol{p}\boldsymbol{k}\right),\boldsymbol{z},\boldsymbol{c}\right) ;
④ \mathbf{i}\mathbf{f} {\boldsymbol{m}}' \ne \perp and \boldsymbol{c}=\boldsymbol{C}\boldsymbol{T}\boldsymbol{R}\boldsymbol{U}. \boldsymbol{P}\boldsymbol{K}\boldsymbol{E}. \boldsymbol{E}\boldsymbol{n}\boldsymbol{c}(\boldsymbol{p}\boldsymbol{k},{\boldsymbol{m}}';\boldsymbol{c}\boldsymbol{o}\boldsymbol{i}\boldsymbol{n}') \mathbf{t}\mathbf{h}\mathbf{e}\mathbf{n}
⑤ \mathbf{r}\mathbf{e}\mathbf{t}\mathbf{u}\mathbf{r}\mathbf{n} {\boldsymbol{K}}' ;
⑥ \mathbf{e}\mathbf{l}\mathbf{s}\mathbf{e}
⑦ \mathbf{r}\mathbf{e}\mathbf{t}\mathbf{u}\mathbf{r}\mathbf{n} \widetilde{\boldsymbol{K}} ;
⑧ end if
2.3 NTT与混合基NTT
快速数论变换是用于加速多项式乘法的一个复杂度为 O\left(n\mathrm{l}\mathrm{o}\mathrm{g}n\right) 的算法[32]. 计算负折叠卷积多项式环乘法 h=f\times g\in {\mathbb{Z}}_{q}\left[x\right]/({x}^{n}+1) , n 为2的幂次,且 2n|q-1 . 可以通过以下NTT变换快速计算:
h={F}_{\mathrm{I}\mathrm{N}\mathrm{T}\mathrm{T}}\left({F}_{\mathrm{N}\mathrm{T}\mathrm{T}}\left(f\right)\circ {F}_{\mathrm{N}\mathrm{T}\mathrm{T}}\left(g\right)\right) \text{,} 其中 {F}_{\mathrm{N}\mathrm{T}\mathrm{T}} 为正向NTT变换, {F}_{\mathrm{I}\mathrm{N}\mathrm{T}\mathrm{T}} 为NTT逆向变换. \circ 为多项式系数之间的逐点相乘,时间复杂度为 O\left(n\right) . 因此负折叠卷积的主要时间转换为NTT变换所需的时间复杂度 O\left(n\mathrm{l}\mathrm{o}\mathrm{g}n\right) . Seiler等人[33]提出了使用FFT trick计算 {\mathbb{Z}}_{q}\left[x\right]/({x}^{n}+1) 环上的NTT变换,正向FFT trick有如下同构:
{\mathbb{Z}}_{q}\left[x\right]/({x}^{n}+1)\cong {\mathbb{Z}}_{q}\left[x\right]/({x}^{\tfrac{n}{2}}-{\zeta }^{\tfrac{n}{2}})\times {\mathbb{Z}}_{q}\left[x\right]/({x}^{\tfrac{n}{2}}+{\zeta }^{\tfrac{n}{2}}) \text{,} 其中 \zeta 为 2n 阶本原单位根,满足 {\zeta }^{2n}\equiv 1\mathrm{m}\mathrm{o}\mathrm{d}\;q .基2 FFT trick每次输入2个系数,输出2个系数,基2 FFT trick的正向变换和逆向变换分别采用Cooley Tukey蝴蝶变换[34]和Gentleman Sande蝴蝶变换[35].Cooley Tukey蝴蝶变换输入 ({f}_{j},{f}_{j+\tfrac{n}{2}}) ,输出 ({f}_{j}+\zeta {f}_{j+\tfrac{n}{2}}, {f}_{j}-\zeta {f}_{j+\tfrac{n}{2}}) ,Gentleman Sande蝴蝶变换输入 ({f}_{j},{f}_{j+\tfrac{n}{2}}) ,输出 ({f}_{j}+{f}_{j+\tfrac{n}{2}},{{\zeta }^{-1}(f}_{j}-{f}_{j+\tfrac{n}{2}})) .
如果多项式环维度并不是2的幂次,且含有因子3的情况,我们可以用基3 FFT trick对多项式环进行NTT变换. 基3 FFT trick变换如下,给定同构 {\mathbb{Z}}_{q}\left[x\right]/({x}^{3m}-{\zeta }^{3})\cong {\mathbb{Z}}_{q}\left[x\right]/({x}^{m}-\zeta )\times {\mathbb{Z}}_{q}\left[x\right]/({x}^{m}-\rho \zeta )\times {\mathbb{Z}}_{q}\left[x\right]/({x}^{m}-{\rho }^{2}\zeta ) ,其中 \zeta 在 {\mathbb{Z}}_{q} 上不可逆,且 \rho 为3阶本原单位根. 对于多项式维度既含有2的幂次因子,又含有3的幂次因子,我们可以采用混合基2加基3的方式对多项式环进行NTT变换. 对于CTRU采用的多项式环,因为不满足NTT多项式环的结构,需要先将其同构到2个可以进行FFT Trick变换的多项式环上,再进行NTT操作,并且对于CTRU-768方案,由于同构后的环维度含有2的幂次因子和3的幂次因子,需要用到混合基NTT,具体方法会在第3节详细介绍.
2.4 Karatsuba算法
Karatsuba算法可以用来加速多项式乘法,对于2个一次多项式 A\left(x\right),B\left(x\right)\in {\mathbb{Z}}_{q}\left[x\right]/\left({x}^{2}-{\zeta }^{r\left(i\right)}\right) ,设 A\left(x\right)= {\mathrm{a}}_{1}x+{a}_{0},B\left(x\right)={b}_{1}x+{b}_{0} , 令{D}_{0}={a}_{0}{b}_{0},{D}_{1}={a}_{1}{b}_{1} , {D}_{\mathrm{0,1}}= ({a}_{0}+{a}_{1})({b}_{0}+{b}_{1}) ,则采用Karatsuba技巧计算如下:
\begin{split} C\left(x\right)=\;&A\left(x\right)\times B\left(x\right)\mathrm{m}\mathrm{o}\mathrm{d}\left({x}^{2}-{\zeta }^{r\left(i\right)}\right)=\\ &{D}_{1}{\zeta }^{r\left(i\right)}+{D}_{0}+\left({D}_{\mathrm{0,1}}-{D}_{0}-{D}_{1}\right)x .\end{split} 我们可以将Karatsuba算法推广到维度为 n 的多项式乘法上,如算法7所示. 具体的正确性分析可见文献[36]. Karatsuba算法的时间复杂度为 {O(n}^{\mathrm{lb}3}) ,其中 n 为多项式的维度. Karatsuba算法一般用于加速维度较小或无法进行NTT的多项式环乘法.
算法7. 度为 n-1 多项式环乘法 basemul\left(a,b\right) .
输入: A\left(x\right) 和 B\left(x\right) 为2个度为 n-1 的多项式 A\left(x\right)={\displaystyle\sum\limits^{n-1}_{i=0}}{a}_{i}{x}^{i},B\left(x\right)={\displaystyle\sum\limits^{n-1}_{i=0}}{b}_{i}{x}^{i} ;
输出:多项式 C\left(x\right)=A\left(x\right)\times B\left(x\right)={\displaystyle\sum\limits^{2n-2}_{i=0}}{c}_{i}{x}^{i}.
① for i from 0 to n do
② {D}_{i}:={a}_{i}{b}_{i} ;
③ end for
④ for i from 1 to 2n−3 do
⑤ {D}_{s,t}:=\left({a}_{s}+{a}_{t}\right)\left({b}_{s}+{b}_{t}\right); /* {s},t 为整数, s+t=i,t > s\ge 0*/
⑥ end for
⑦ {c}_{0}={D}_{0} ;
⑧ {c}_{2n-2}={D}_{n-1} ;
⑨ {c}_{i}= \left\{\begin{aligned}&\sum _{s+t=i;t > s\ge 0}{D}_{s,t}-\sum _{s+t=i;t > s\ge 0}\left({D}_{s}+{D}_{t}\right),奇数i.\\&\sum _{s+t=i;t > s\ge 0}{D}_{s,t}-\sum _{s+t=i;t > s\ge 0}\left({D}_{s}+{D}_{t}\right)+{D}_{\frac{i}{2}},偶数i.\end{aligned}\right.
2.5 Bernstein快速求逆算法
Bernstein等人[37]提出一种高效的常数时间实现求逆算法. 如算法8所示,算法8中需要使用到函数 reverse 和函数 swap . 函数 reverse 为 reverse\left(f,d\right)= {x}^{d}\times f(1/x) . 函数 swap 为 swap(f,g) ,互换多项式 f 和 g 的系数值. 另外,算法8中使用符号 \Delta 表示某些数值,其余符号 \varPhi ,F,V,S 均表示某些多项式. Bernstein快速求逆方法适用于维度较大的多项式环求逆计算.
算法8. 多项式求逆算法.
输入:模多项式 \phi ,待求逆的多项式 f ;
输出: {f}_{\mathrm{i}\mathrm{n}\mathrm{v}}:={f}^{-1}\mathrm{m}\mathrm{o}\mathrm{d}\; \phi .
① {\varPhi }_{0}:=reverse\left(\phi ,n\right) ;
② {F}_{0}:=reverse\left(f,n-1\right) ;
③ {\Delta }_{0}:=1 ;
④ {V}_{0}:=0 ;
⑤ {S}_{0}:=1 ;
⑥ \mathrm{f}\mathrm{o}\mathrm{r} i=\mathrm{0,1},2,… ,2n-2 \mathrm{d}\mathrm{o}
⑦ {V}_{i}:={V}_{i}\times x ;
⑧ \mathrm{i}\mathrm{f} {\Delta }_{i} > 0 \mathrm{a}\mathrm{n}\mathrm{d} {F}_{i}\left(0\right)\ne 0 \mathrm{t}\mathrm{h}\mathrm{e}\mathrm{n}
⑨ {\Delta }_{i}:=-{\Delta }_{i} ;
⑩ swap\left({\Phi }_{i},{F}_{i}\right) ;
⑪ swap\left({V}_{i},{S}_{i}\right) ;
⑫ end if
⑬ {\Delta }_{i+1}:={\Delta }_{i}+1 ;
⑭ {F}_{i+1}:=\dfrac{\left[{\Phi }_{i}\left(0\right)\times {F}_{i}-{F}_{i}\left(0\right)\times {\Phi }_{i}\right]}{x} ;
⑮ {S}_{i+1}:={\varPhi }_{i}\left(0\right)\times {S}_{i}-{F}_{i}\left(0\right)\times {V}_{i} ;
⑯ {\varPhi }_{i+1}:={\varPhi }_{i} ;
⑰ {V}_{i+1}:={V}_{i} ;
⑱ \mathrm{e}\mathrm{n}\mathrm{d}\mathrm{ }\mathrm{f}\mathrm{o}\mathrm{r}
⑲ {f}_{\mathrm{i}\mathrm{n}\mathrm{v}}:=reverse\left(\dfrac{{V}_{2n-1}}{{F}_{2n-1}\left(0\right)},n-1\right) ;
⑳ \mathrm{r}\mathrm{e}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{n} {f}_{\mathrm{i}\mathrm{n}\mathrm{v}} .
2.6 实现平台
本文在Intel x86平台上使用AVX2指令集[38]对CTRU 3个方案进行优化实现. AVX2含有16个256 b向量寄存器(ymm0-ymm15),可以同时对16个16 b数进行并行操作. 特别适用于一些非串行程序的向量化实现,本文使用的主要AVX2指令如表3所示,这些指令对于主要耗时操作多项式乘法实现具有很大帮助.
表 3 CTRU实现主要AVX2指令Table 3. Main AVX2 Instructions Implemented by CTRU汇编指令 Intrinsic函数 描述 vpaddw _mm256_add_epi16 16b并行整数加法 vpsllw _mm256_slli_epi16 16b并行整数减法 vpsraw _mm256_srai_epi16 16b并行整数移位 vpsubw _mm256_sub_epi16 16b并行整数减法 vpmulhw _mm256_mulhi_epi16 16b并行整数乘法取高32b vpblendd _mm_blend_epi32 通过mask对32b数据
顺序进行调整vpunpcklqdq _mm256_unpacklo_epi64 拆包并交错低64b数 vpunpckhqdq mm256_unpackhi_epi64 拆包并交错高64b数 vpbroadcastw _mm256_set1_epi16 广播16b数 3. CTRU方案参考实现设计
在CTRU方案中,密钥生成过程中的多项式求逆和密钥生成、密钥封装、密钥解封装都存在的多项式乘法是核心模块. 我们从这2个运算模块进行展开,给出CTRU方案的多项式算术运算的具体实施方法和采用的对应用于加速该模块的算法.
3.1 多项式乘法
对于CTRU方案的多项式环 f={\mathbb{Z}}_{q}\left[x\right]/({x}^{n}- {x}^{n/2}+1) ,可以将其同构到环 {\mathbb{Z}}_{q}\left[x\right]/({x}^{n/2}-{\zeta }_{1}) 和环 {\mathbb{Z}}_{q}\left[x\right]/ ({x}^{n/2}-{\zeta }_{2}) . 其中 {\mathrm{\zeta }}_{1} 和 {\mathrm{\zeta }}_{2} 满足 {\zeta }_{1}\times {\zeta }_{2}=1,{\zeta }_{1}+{\zeta }_{2}=1 . 由 f\mathrm{m}\mathrm{o}\mathrm{d}\left({x}^{n/2}-{\zeta }_{1}\right)=\left({f}_{0}+{\zeta }_{1}{f}_{n/2}\right)+\dots +\left({f}_{\frac{n}{2}-1}+{\zeta }_{1}{f}_{n-1}\right){x}^{\frac{n}{2}-1} ; f\mathrm{m}\mathrm{o}\mathrm{d}\left({x}^{n/2}-{\zeta }_{2}\right)=\left({f}_{0}+{\zeta }_{2}{f}_{n/2}\right)+\dots +\left({f}_{\frac{n}{2}-1}+{\zeta }_{2}{f}_{n-1}\right){x}^{\frac{n}{2}-1} .
我们将多项式环进行分解,分解后的环 {\mathbb{Z}}_{q}\left[x\right]/({x}^{n/2}-{\zeta }_{1}) 和 {\mathbb{Z}}_{q}\left[x\right]/({x}^{n/2}-{\zeta }_{2}) 可以使用中国剩余定理(Chinese remainder theory, CRT)进行基-2 FFT步骤或基-3 FFT步骤的分解. 每个基-2 FFT trick如下: {\mathbb{Z}}_{q}\left[x\right]/({x}^{2m}-{r}^{2})\cong {\mathbb{Z}}_{q}\left[x\right]/({x}^{m}-r)\times {\mathbb{Z}}_{q}\left[x\right]/({x}^{m}+r) . 因此CTRU方案多项式环上的乘法就转化为正向NTT变换、多项式点乘以及逆向NTT变换3部分.
3.1.1 NTT
本节根据具体参数分析FFT trick的实施流程. 对于 n=512 , q=3 457 ,本文选取384阶本原单位根 \mathrm{\zeta }=55 , {\mathrm{\zeta }}_{1}={\zeta }^{n/8}={\zeta }^{64}\mathrm{m}\mathrm{o}\mathrm{d}\;q,{\zeta }_{2}={\zeta }^{5n/8}{=\zeta }^{320}\mathrm{m}\mathrm{o}\mathrm{d}\;q . 本文对 {\mathbb{Z}}_{q}[x]/({x}^{256}-{\zeta }^{64}) 和 {\mathbb{Z}}_{q}[x]/({x}^{256}-{\zeta }^{320}) 这2个环进行基-2FFT trick分解,这2个环的基-2 FFT trick步骤可以进行6层,直到分解到 64 个模 4 次多项式的项,例如 {\mathbb{Z}}_{q}[x]/({x}^{4}-{\zeta }^{5}) . 对于 n=768, 本文选取
1152 阶本原单位根 \zeta =5 , {\mathrm{\zeta }}_{1}={\zeta }^{n/4}={\zeta }^{192}\mathrm{m}\mathrm{o}\mathrm{d}\;q,{\zeta }_{2}={\zeta }^{5n/4}{=\zeta }^{960}\mathrm{m}\mathrm{o}\mathrm{d}\;q . 对 {\mathbb{Z}}_{q}[x]/({x}^{384}-{\zeta }^{192}) 和 {\mathbb{Z}}_{q}[x]/({x}^{384}-{\zeta }^{960}) 这2个环同样进行6层基-2 FFT trick分解,分解到环 {\mathbb{Z}}_{q}[x]/ ({x}^{6}-{\zeta }^{3}) ,此时还能进一步用基-3 FFT trick进行分解,分解到 {\mathbb{Z}}_{q}[x]/({x}^{2}-{\zeta }^{r(i)}) . 具体地说,对于环 {\mathbb{Z}}_{q}[x]/({x}^{6}-{\zeta }^{3}) 存在以下同构 {\mathbb{Z}}_{q}[x]/({x}^{6}-{\zeta }^{3})\cong {\mathbb{Z}}_{q}[x]/ ({x}^{2}-\zeta )\times {\mathbb{Z}}_{q}[x]/({x}^{2}-\rho \zeta )\times {\mathbb{Z}}_{q}[x]/({x}^{2}-{\rho }^{2}\zeta ) . 其中 \rho ={\zeta }^{n/3}\mathrm{m}\mathrm{o}\mathrm{d}\;q 为3次单位根,满足 {\rho }^{2}+\rho +1\equiv 0\;\mathrm{m}\mathrm{o}\mathrm{d}\;q . 例如 {\mathbb{Z}}_{q}[x]/({x}^{6}-{\zeta }^{3}) 中的多项式表示为 f={f}_{x}+{f}_{y}{x}^{2}+{f}_{z}{x}^{4}, 其中 {f}_{x},{f}_{y},{f}_{z} 为一次多项式,令基-3 NTT分解后的多项式为 {f}_{a},{f}_{b},{f}_{c} . 可得分解后的3个一次多项式{f}_{a}=f\mathrm{m}\mathrm{o}\mathrm{d}({x}^{2}-\zeta )={f}_{x}+\zeta {f}_{y}+{\zeta }^{2}{f}_{z},{f}_{b}= f\mathrm{m}\mathrm{o}\mathrm{d}({x}^{2}-\rho \zeta )={f}_{x}+\rho \zeta {f}_{y}+{(\rho \zeta )}^{2}{f}_{z}, {f}_{c}={f}''\mathrm{m}\mathrm{o}\mathrm{d}({x}^{2}-{\rho }^{2}\zeta )= {f}_{x}+{\rho }^{2}\zeta {f}_{y}+{({\rho }^{2}\zeta )}^{2}{f}_{z} ,对于 n=1 024 ,本文同样选取384阶本原单位根, {\zeta }^{n/16}={\zeta }^{64}\mathrm{m}\mathrm{o}\mathrm{d}\;q,{\zeta }_{2}={\zeta }^{5n/16}{=\zeta }^{320}\mathrm{m}\mathrm{o}\mathrm{d}\;q ,对环 {\mathbb{Z}}_{q}[x]/({x}^{512}-{\zeta }^{64}) 和环 {\mathbb{Z}}_{q}[x]/({x}^{512}-{\zeta }^{320}) 进行6层基-2 FFT分解到 {\mathbb{Z}}_{q}[x]/({x}^{8}-{\zeta }^{r(i)}) . 根据以上分析,可以画出3个多项式环对应的正向FFT trick分解图,如图1,2,3所示:3.1.2 多项式点乘
在进行正向之后,CTRU-512,CTRU-768,CTRU-1024方案的多项式环分解到 {\mathbb{Z}}_{q}[x]/({x}^{,4}-{\zeta }^{r(i)}) 、 {\mathbb{Z}}_{q}[x]/ ({x}^{2}-{\zeta }^{r(i)}) 和 {\mathbb{Z}}_{q}[x]/({x}^{8}-{\zeta }^{r(i)}) . 因此需要在这3个环上进行多项式乘法操作,这里本文采用Karatsuba算法加速多项式乘法. 与Schoolbook乘法相比,Karatsuba技巧减少了一个乘法操作,增加了加法操作;在Intel CPU上的C参考实现和AVX2实现中,乘法一般是较为耗时的,因此Karatsuba可以有效加速C实现和AVX2实现. {n}=512 多项式环需要进行 n/4 个 {\mathbb{Z}}_{q}[x]/ ({x}^{4}-{\zeta }^{r(i)}) 多项式环点乘, n=768 需要进行 n/2 个 {\mathbb{Z}}_{q}[x]/ ({x}^{2}-{\zeta }^{r(i)}) 多项式环点乘, n=1 024 需要进行 n/8 个 {\mathbb{Z}}_{q}[x]/ ({x}^{8}-{\zeta }^{r(i)}) 多项式环点乘. 而这3个环采用Schoolbook算法进行多项式乘法运算时的时间复杂度分别为 O(16) , O(4) , O(64) . (在这里我们分析时间复杂度时只考虑乘法的个数,因为相较于乘法运算,加减法运算在Intel CPU软件实现中运算基本不耗时). Karatsuba算法进行多项式乘法的时间复杂度为 O(10) , O(3) , O(22) . 因此我们可以得出3个方案进行多项式点乘的总时间复杂度如表4所示,可以看出采用Karatsuba算法可以进一步降低多项式点乘的时间复杂度,从而加速多项式点乘运算.
表 4 CTRU 3个方案多项式逐点相乘时间复杂度分析Table 4. Time Complexity Analysis of Polynomial Pointwise Multiplication in Three Schemes of CTRUn 多项式环 Karatsuba算法 Schoolbook算法 512 {\mathbb{Z}}_{q}\left[x\right]/\left({x}^{4}-{\zeta }^{r\left(i\right)}\right) O\left(\dfrac{5}{2}n\right) O\left(4n\right) 768 {\mathbb{Z}}_{q}\left[x\right]/\left({x}^{2}-{\zeta }^{r\left(i\right)}\right) O\left(\dfrac{3}{2}n\right) O\left(2n\right) 1 024 {\mathbb{Z}}_{q}\left[x\right]/\left({x}^{8}-{\zeta }^{r\left(i\right)}\right) O\left(\dfrac{11}{4}n\right) O\left(8n\right) 3.2 多项式乘法时间复杂度分析
根据之前提出的CTRU方案多项式环不同参数下多项式环的NTT变换操作. 我们可以给出CTRU方案不同参数下的多项式乘法时间复杂度,总体来说,采用NTT加速的多项式乘法包含2个正向NTT操作一个多项式点乘操作和一个多项式求逆,由于CTRU-512和CTRU-1024方案都是采用6层基2 NTT方法,每层乘法个数为 n/2 ,因此CTRU-512和CTRU-1024方案的NTT时间复杂度为 O\left(3n\right) ,算上正向、逆向NTT,NTT变换的总时间复杂度为 O\left(9n\right) . 而CTRU-768方案采用的是6层基2加1层基3混合基NTT,需要算上一层基-3 NTT的时间复杂度,一层基-3 NTT的乘法个数为 O(3\times \dfrac{n}{3}=n) . 因此CTRU-768方案NTT部分的时间复杂度为 O(3n+n=4n) ,总时间复杂度为 O(3\times 4n=12n) . 再加上3.1.2节分析的多项式点乘的时间复杂度,我们可以得到CTRU 3个方案的多项式乘法时间复杂度,如表5所示,相较于Schoolbook算法,使用NTT算法可以大幅度的降低多项式乘法的时间复杂度.
表 5 CTRU 3个方案多项式乘法时间复杂度分析Table 5. Time Complexity Analysis of Polynomial Multiplication in three Schemes of CTRU方案 NTT算法 Schoolbook算法 CTRU-512 O\left(9n+\dfrac{5}{2}n\right) {O(n}^{2}) CTRU-768 O\left(12n+\dfrac{3}{2}n\right) O\left({n}^{2}\right) CTRU- 1024 O\left(9n+\dfrac{11}{4}n\right) O\left({n}^{2}\right) 3.3 约减算法
考虑到CTRU方案的模数,我们在参考实现中选取16b有符号数存储CTRU方案的系数,在多项式算术运算过程中,存在对系数的加减和乘法操作,使得运算结果可能会超出16b有符号数范围,为了防止系数的溢出,需要对系数进行约减,Barrett约减[39]和Montgomery约减[40]是2种高效且常数时间的约减算法.Seiler给出了这2种算法在有符号整数下的优化版本和其参考实现[33].Barrett约减的输出范围是 (0,q) ,而Montgomery约减的输出范围是 \left(-q,q\right). 由于Montgomery约减算法是对MONT域范围的整数进行约减,因此在约减结束后还需要将结果转换到正常域. 在正向NTT和逆向NTT过程中, \beta 值可以提前预计算在旋转因子表中,不需要额外进行域的转换. 由于CTRU方案的模数在C语言16b有符号数范围内,在本文实现中选取 \beta ={2}^{16} .
为了更好地实施延迟约减,减少约减次数. 本文采用文献[41]中提到的输出范围在 [-\dfrac{q}{2},\dfrac{q}{2}) 的中心Barrett约减算法对系数进行约减.
在CTRU的解密过程中,需要用到多模数NTT方法优化模数 {q}_{2}=1 024 的非素数模数多项式乘法,根据文献[8],选取 q=3 457 和 {q}'=7 681. 7 681可以用2的幂次项表示: {q}'=7 681={2}^{13}-{2}^{9}+1. 我们根据模数 q' 的2的幂次项之和表示形式给出针对其的特殊约减算法,只需要移位和加减操作就能完成约减计算,减少约减时间.
算法9. 对于模数 q=7681 的特殊约减.
输入: -{2}^{15}\le a < {2}^{15} ;
输出: r\equiv a\left(\mathrm{m}\mathrm{o}\mathrm{d}\;q\right), 其中 -{2}^{15}+4q\le r < {2}^{15}-3q.
① t\leftarrow \left\lfloor {a/{2}^{13}} \right\rfloor ;
② u\leftarrow a\;\mathrm{m}\mathrm{o}\mathrm{d}\;{2}^{13} ;
③ u\leftarrow u-t ;
④ \mathrm{t}\leftarrow t\times {2}^{9} ;
⑤ r\leftarrow u+t .
3.4 系数范围分析与延迟约减
为了最大程度减少约减的次数,我们在NTT过程中对系数进行延迟约减,通过分析系数范围尽可能推迟约减的层数.CTRU的NTT使用16b有符号整数,可存储最大范围为 -32768\le \boldsymbol{k}\boldsymbol{q}\le 32767 . 其中模数 \boldsymbol{q}=3457 ,因此 \boldsymbol{k}=9. 由于CTRU环的特殊性,进行正向基-2 NTT之前还要进行分解到环 {\mathbb{Z}}_{\boldsymbol{q}}\left[\boldsymbol{x}\right]/ ({\boldsymbol{x}}^{\boldsymbol{n}/2}-{\boldsymbol{\zeta }}_{1}) 和环 {\mathbb{Z}}_{\boldsymbol{q}}\left[\boldsymbol{x}\right]/({\boldsymbol{x}}^{\boldsymbol{n}/2}-{\boldsymbol{\zeta }}_{2}) 的操作,此次分解后系数增长到 \left[-2\boldsymbol{q},2\boldsymbol{q}\right] . 正向NTT的蝴蝶操作为Cooley-Tukey蝴蝶操作: (\boldsymbol{a},\boldsymbol{b})\to (\boldsymbol{a}+\boldsymbol{\zeta }\boldsymbol{b},\boldsymbol{a}-\boldsymbol{\zeta }\boldsymbol{b}) ,6层基-2 NTT操作结束,系数范围为 [-8\boldsymbol{q},8\boldsymbol{q}], 不会溢出. 如果只进行6层基-2 NTT,则正向NTT不需要进行Barrett约减. CTRU-768的多项式环正向NTT时还需要进行1层基-3 NTT,因此在正向基-2 NTT结束后需要对所有系数进行Barrett约减.
我们在逆向NTT第1层之前对一半系数做中心Barrett约减,因此在进行逆向基-2 NTT第1层之后的系数范围为 [-q,q) . 逆向NTT进行的是Gentleman-Sande蝴蝶操作: (a,b)\to (a+b,\zeta \left(a-b\right)) ,后半部分系数每次都会进行Montgomery约减到 \left(-q,q\right). 因此,并不是每一层中所有的系数都需要进行Barrett约减,此时,可以对不同层固定的下标索引进行分析. 在纯基-2 NTT中,需要约减的层数、系数索引和对应Barrett约减个数如表6,7,8所示. 在基-2、基-3混合基NTT中,本文采取基-3 NTT后对所有系数做Barrett约减,再按刚才分析的基-2 NTT约减情况进行约减. 在 n=512 时Barrett约总次数为256+128=384次,在 n=768 时Barrett约减次数为384+192=576次,在 n=1 024 时Barrett约减次数为512+640=768次.
表 6 n=512 逆NTT Barrett约减分析Table 6. \boldsymbol{n}=512 Inverse NTT Barrett Reduction Analysis层数 索引 约减次数 第 4 层 64\to 71,192\to 199
320\to 327,448\to 45532 第 5 层 0\to 7,128\to 135
136\to 143
256\to 263,384\to 391
392\to 39948 第 6 层 8\to 15,264\to 271
16\to 31,272\to 28748 表 7 n=768 逆NTT Barrett约减分析Table 7. \boldsymbol{n}=768 Inverse NTT Barrett Reduction Analysis层数 索引 约减次数 第 4 层 96\to 107,480\to 491
288\to 299,672\to 68248 第 5 层 0\to 11,384\to 395
204\to 215,588\to 599
192\to 203,576\to 58772 第 6 层 12\to 23,396\to 407
24\to 47,408\to 43172 表 8 n=1 024 逆NTT Barrett约减分析Table 8. \boldsymbol{n}=1024 Inverse NTT Barrett Reduction Analysis层数 索引 约减次数 第 4 层 128\to 143,640\to 655
384\to 399,896\to 91164 第 5 层 0\to 15,512\to 527
272\to 287,784\to 799
256\to 271,768\to 78396 第 6 层 16\to 31,528\to 543
32\to 63,544\to 57596 3.5 多项式求逆
多项式求逆是CTRU方案中的核心多项式算术运算之一,出现在CTRU方案的密钥生成过程中,需要先计算多项式 f 在 {R}_{q} 中的逆 {f}^{-1} . 但由于多项式 f 先进行了NTT变换,我们可以转向对NTT变换后分解的子多项式进行求逆操作,对于 n=\{512,768,1\mathrm{ }024\} ,正向NTT后分别分解到环 {\mathbb{Z}}_{q}[x]/({x}^{4}-{\zeta }^{r(i)}) ,环 {\mathbb{Z}}_{q}[x]/ ({x}^{2}-{\zeta }^{r(i)}) ,和环 {\mathbb{Z}}_{q}[x]/({x}^{8}-{\zeta }^{r(i)}) . 因此多项式环 {\mathbb{Z}}_{q}[x]/ ({x}^{n}-{x}^{n/2}+1) 上的求逆可以转化为 n/4 个 {\mathbb{Z}}_{q}[x]/ ({x}^{4}-{\zeta }^{r(i)}) , n/2 个 {\mathbb{Z}}_{q}[x]/({x}^{2}-{\zeta }^{r(i)}) 和 n/8 个 {\mathbb{Z}}_{q}[x]/({x}^{8}-{\zeta }^{r(i)}) 多项式环上的求逆. 对于多项式环 {\mathbb{Z}}_{q}[x]/({x}^{2}-{\zeta }^{r(i)}) 和 {\mathbb{Z}}_{q}[x]/({x}^{4}-{\zeta }^{r(i)}) 的求逆计算可以用Cramer’s 方法[42]计算. 以 {\mathbb{Z}}_{q}[x]/({x}^{4}-\zeta ) 为例. 令 f\in {\mathbb{Z}}_{q}[x]/({x}^{4}-\zeta ), 其多项式的逆定义为 {f}^{-1} ,有 f\times {f}^{-1}=1\mathrm{m}\mathrm{o}\mathrm{d}({x}^{4}-\zeta ). 本文可以将多项式乘法 f\times {f}^{-1} 写成 {\mathbb{Z}}_{q} 上的矩阵向量乘法,其中包含多项式 f 的 4\times 4 阶旋转矩阵:
\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right]=\left[\begin{array}{cccc}{f}_{0}& \zeta {f}_{3}& \zeta {f}_{2}& \zeta {f}_{1}\\ {f}_{1}& {f}_{0}& \zeta {f}_{3}& \zeta {f}_{2}\\ {f}_{2}& {f}_{1}& {f}_{0}& \zeta {f}_{3}\\ {f}_{3}& {f}_{2}& {f}_{1}& {f}_{0}\end{array}\right]\cdot \left[\begin{array}{c}{f}_{0}^{-1}\\ {f}_{1}^{-1}\\ {f}_{2}^{-1}\\ {f}_{3}^{-1}\end{array}\right]. 计算 f 的逆可以转化为计算旋转矩阵的逆. 求旋转矩阵的逆需要求旋转矩阵行列式的逆,若行列式的逆不存在,则逆矩阵也不存在. 令 d 为系数矩阵的行列式,根据Cramer法则,存在唯一的 {f}^{-1}, 其中每个系数可以通过如下公式计算:
{f}_{i}^{-1}=\frac{{d}_{i}}{d},i=0,1,2,3. {d}_{i} 是替换系数矩阵的第 i+1 列为 {(\mathrm{1,0},\mathrm{0,0})}^{\mathrm{T}} 后得到的矩阵的行列式. {d}^{-1} 可以通过费马小定理计算: {d}^{-1}\equiv {d}^{q-2}={d}^{3 455}\mathrm{m}\mathrm{o}\mathrm{d} 3 457. 同样 {\mathbb{Z}}_{q}[x]/({x}^{2}-\zeta ) 的多项式的逆也可以按照如上方法计算. 但是在 {\mathbb{Z}}_{q}[x]/ ({x}^{8}-\zeta ) 下,矩阵为 8\times 8 阶矩阵,此时矩阵行列式的项数为 8!=40\;320 . 由于项数过多,给实现带来了困难. 因此对于 {\mathbb{Z}}_{q}[x]/({x}^{8}-\zeta ) 的多项式求逆,本文采用Bernstein快速求逆方法,Bernstein比传统的扩展欧几里得算法要快,适用于加速 {\mathbb{Z}}_{q}[x]/({x}^{8}-\zeta ) 下无法采用多项式矩阵求逆方法的多项式求逆实现. 为了分析Bernstein快速求逆和多项式矩阵求逆在CTRU-512和CTRU-768这2个方案下哪个方法更好. 我们对CTRU-512和CTRU-768这2个方案的多项式求逆分别用Bernstein快速求逆方法和多项式矩阵求逆方法实现了多项式求逆,结果如表9所示.
表 9 多项式求逆方性能比较Table 9. Performance Comparison of Polynomial Inversion Techniques方案 Bernstein快速求逆/时钟周期 多项式矩阵求逆/时钟周期 CTRU-512 129 352 91 860 CTRU-768 195 300 72 234 从表9中可以看出,对于CTRU-512和CTRU-768这2个方案来说,采用Bernstein快速求逆方法还是要远慢于多项式矩阵求逆方法,因此对于CTRU-512和CTRU-768这2个方案我们采取多项式矩阵求逆方法进行加速.
4. CTRU方案AVX2优化实现设计
4.1 性能分析
密码算法的AVX2优化实现相较于参考实现,主要目的在于提高算法的性能与效率. 其核心思想在于针对参考实现中消耗CPU时钟周期较多的函数,即性能瓶颈部分,进行优化改进. 因此本节给出了perf工具下CTRU-512,CTRU-768,CTRU-1024这3个方案参考实现的软件分析测试结果,该性能分析数据通过运行CTRU-768代码1 000次并计算平均执行时间;结果如表10所示.
表 10 CTRU-512/768/1024 方案时间占比分析Table 10. Analysis of the Time Proportion of CTRU-512/768/1024 Schemes函数 CTRU-512 CTRU-768 CTRU- 1024 Montgomery reduce 31.60 24.31 13.29 Poly Decode 7.30 15.50 10.08 Base Inversion 3.84 13.53 30.08 {F}_{\mathrm{N}\mathrm{T}\mathrm{T}} 10.55 10.92 7.80 Base multiplication 2.16 2.45 3.42 {F}_{\mathrm{I}\mathrm{N}\mathrm{T}\mathrm{T}} 5.10 6.34 3.26 KeccakF1600_StatePermutel 6.95 8.09 4.49 Barrett reduce 15.86 4.37 16.91 其中Poly Decode是CTRU算法中的解码函数,Base Inversion是多项式求逆函数,KeccakF1600_StatePermute是SHA-3哈希函数的轮函数,正向NTT和逆向NTT分别表示正向数论变换和逆向数论变换,Base multiplication表示多项式点乘. 从表10中可以看出,在3个方案的参考实现中,Montgomery约减,正向NTT和逆向NTT以及多项式点乘、多项式求逆都是较为耗时的操作,特别地,因为CTRU-
1024 方案只能采取Bernstein快速求逆方法,多项式求逆成为了耗时最多的运算.根据以上分析,本文对CTRU算法主要耗时模块进行AVX2优化实现,主要优化点如下:
1)给出了CTRU-512,CTRU-768,CTRU-1024 3个方案下多项式环乘法的实现,之前的CTRU-768方案实现并没有针对哈希部分进行AVX2优化实现,并且也没有基于索引的延迟约减,针对已有参数实现,本文进行优化,采用基于索引的延迟约减,并对耗时模块哈希部分进行优化. 对于CTRU-512和CTRU-1024方案,本文给出首个AVX2实现,采用前面提到的基于索引的延迟约减,并采用层融合、系数置乱等方法加速正向NTT和逆向NTT.
2)给出了CTRU-512,CTRU-768,CTRU-
1024 这3个方案下多项式环求逆的AVX2优化实现,提出了首个Bernstein快速多项式求逆的AVX2向量化实现.3)给出了多项式采样、多项式算术运算(多项式相加、多项式相减)等函数的AVX2优化实现.
4)给出了多项式序列化和反序列化的AVX2优化实现.
5)哈希函数SHA-3部分进行AVX2优化实现.
4.2 NTT AVX2优化实现
在NTT的AVX2优化实现中,关键技术有NTT变换时的层融合、系数置乱操作,蝴蝶操作的实现. 本节也将从这几个方面展开讨论CTRU NTT AVX2优化实现细节. 由于CTRU-768方案在文献[8]中已经给出AVX2的实现细节. 本节中主要给出其余2组的AVX2优化实现细节.
4.2.1 多项式存储格式
本文将每个含有 n 个系数的多项式表示为长度为 n 的16b有符号数整型数组. 为了降低访存时间,我们将每个系数数组进行32B的对齐. 一个256b ymm向量寄存器可以存储16个16b数,达到16并行效果.
4.2.2 层融合和系数置乱
由于AVX2向量寄存器有限,不能一次性将所有的多项式系数载入向量寄存器,在正向NTT和逆向NTT的AVX2汇编实现中,本文需要采用层融合技术. 层融合技术是指对加载到向量寄存器的系数完成多层计算再存回内存. n=512 和 n=1 024 都是6层基2-NTT,因此这里本文以 n=512 为例讲述在NTT中如何应用层融合技术. 本文使用vmovdqa指令从内存中载入连续的16个系数. 对于基-2 NTT前的分解操作,需要对间隔为256的系数进行蝴蝶操作. 本文需要将间隔为256的系数载入到向量寄存器. 如图4所示.1组向量寄存器对中的系数完成相应运算后写回内存.
按图4顺序载入系数,共需要8个向量寄存器;并且环分解操作和基-2 NTT的第1层到第2层的系数对都在寄存器中,因此可以进行融合操作,以减少访存次数,本文采用调用宏加传入偏移地址的方法实现多项式所有系数的运算.
环分解操作及基-2 NTT 第1层到第2层. 对于正向基-2 NTT的前2层,本文一次性载入128个系数到8个向量寄存器,载入系数顺序如图5所示.
前2层系数对蝴蝶操作间隔分别为128和64,按照图5所示顺序载入系数刚好可以完成前2层系数对间隔的蝴蝶操作. 基-2 NTT第3层到第6层. NTT后4层载入系数顺序如图6所示. 本文使用8个向量寄存器一次性载入128个系数.
对于系数置乱,后4层的间隔分别为32,16,8,4;由于只能载入连续的16个系数到向量寄存器,因此上述存放顺序无法完成相应系数间隔的蝴蝶操作,需要调整向量寄存器中系数顺序,这里本文称调整向量寄存器中系数顺序的操作为Shuffle.对目前已经载入向量寄存器的系数做Shuffle操作,本文需要对每8,4,2,1个系数进行Shuffle,这里本文称为Shuffle8,Shuffle4,Shuffle2,Shuffle1.后三者的shuffle操作流程如图7所示.
Shuffle8的主要操作如图8所示,取0~15的前8个系数和32~47的前8个系数拼起来放到1个向量寄存器中,后8个系数8~15和40~47放到1个寄存器中.Shuffle8后寄存器中存放系数顺序如图9所示,两两寄存器之间的连线和蝴蝶表示这2对寄存器之间完成1对蝴蝶操作,对应第5层需要进行间隔8的蝴蝶操作.
对于第6层,因为需要对间隔为4的系数进行蝴蝶操作,以上寄存器中系数顺序无法满足蝴蝶操作需求,需要进一步对寄存器中系数进行Shuffle4. Shuffle4操作如图8所示.
为了后续多项式点乘方便,本文在正向NTT结束后又对系数进行了Shuffle2和Shuffle1的调整,具体原因在4.3节会讲解,Shuffle2,Shuffle1的操作如图10和图11所示. 本文使用vperm2i128指令,通过设置mask值为0x20和0x31达到Shuffle8的效果. 对于Shuffle4操作,AVX2的vpunpcklqdq和vpunpckhqdq指令可以完成取出2个寄存器中每128位中的低64b和高64b. Shuffle2本文使用移位指令vpsllq和vpsrlq指令调整每64b中32b数据顺序,最后通过vpblendd指令设置mask为0xAA,将2个寄存器中的每低32b进行拼接. Shuffle1的实现思路和Shuffle2类似,采用移位指令和vpblendw指令完成. 对于 n=1 024 多项式环NTT, AVX2实现与 n=512 的实现思路类似.
4.2.3 蝴蝶操作
蝴蝶变换中系数和旋转因子进行乘法运算后做Montgomery约减,为了对多项式系数乘积 a 进行模约减操作,本文首先需要计算 m={aq}^{-1}\mathrm{m}\mathrm{o}\mathrm{d}\beta ,这里 \beta ={2}^{16}, 因为旋转因子是常数,本文可以将Montgomery约减中第2步乘 {q}^{-1} 操作预计算进旋转因子表中,这样可以省去1条乘法指令. 本文使用vpmullw指令和vpmulhw指令计算2个16 b整数乘积的高16 b和低16 b,这样避免了先计算完整32 b乘法再取高低位操作.
4.2.4 延迟约减
本文采取第3节所提到的基于索引的延迟约减策略,但是因为1次只能对一整个向量寄存器进行处理,会存在额外对部分系数进行约减的情况. 但总体上比对延迟约减层所有系数进行约减的次数更少.
4.3 多项式点乘AVX2实现
前面提到 n=512 时多项式点乘是2个一次多项式之间的乘法,具体而言,是在 {\mathbb{Z}}_{q}\left[x\right]/({x}^{4}-{\zeta }^{t\left(i\right)}) 中计算2个多项式的模乘. 在AVX2向量化实现中,本文只能对2个向量寄存器进行运算操作,同一向量寄存器中位于不同位置的数据无法进行运算. 具体地说,如果本文按照0,1,2,…的顺序存储并加载到向量寄存器,那么无法对第0个系数和第1个系数直接做运算. 因此需要对系数顺序进行额外的置乱操作. 另一种思路是在做完正向NTT后,保持最后Shuffle4结束的顺序而不调整回自然顺序,则输入多项式点乘是乱序,但只要保证对应位置进行相乘,则多项式乘法结果是正确的. 因此在之前顺序的基础上再做1次Shuffle2和Shuffle1,使得间隔为1的系数不在同一个向量寄存器中.
这样做可以减少逆向NTT变换时候的系数置乱操作,方便了点乘和多项式求逆的运算. 在 n=1 024 进行点乘时也是一样思路.
4.4 多项式求逆AVX2实现
多项式矩阵求逆 AVX2实现在之前的文献中已经给出,本文主要给出Bernstein快速实现的AVX2优化实现思路. 由于Bernstein快速多项式求逆的特性,并不是很适合直接的向量化. 但多项式求逆又是整个算法十分耗时的一个部分,因此本文采取一些技术对Bernstein快速多项式求逆进行了首个AVX2向量化实现,本文提出的Bernstein快速多项式求逆AVX2实现并行处理16个多项式的求逆运算. 由于寄存器数量有限,因此开辟栈空间保存循环中会使用到的变量 P,F,V,S,Delta,Swap ,中间值 uq,ur 和返回值 r . 这些变量对应的作用和在栈中的地址如表11所示.
表 11 多项式求逆实现变量作用及地址Table 11. Polynomial inversion realizes variable function and address变量 作用 地址 P 保存旋转因子 0~287 F 保存原始多项式的系数 288~575 V 逆多项式系数的中间值 576~863 S 逆多项式系数的中间值 864~1 151 Delta 计算返回值 r 的参数 1 152~1 215 Swap 交换 Phi 与 F , V 与 S 1 216~1 311 Ur 计算乘积之差的中间值 1 312~1 343 Uq 计算乘积之差的中间值 1 344~1 407 r 返回值 1 408~1 409 将堆栈中数据按照图12进行初始化, {P}_{i} 列表示第 i 个多项式的变量, {z}_{i} 表示第 i 个多项式对应的旋转因子, {f}_{ij} 表示第 j 个多项式第 i 个系数. D{a}_{i} , S{p}_{i} , U{r}_{i} , U{q}_{i} 分别表示 Delta , Swap , Ur , Uq 的第 i 个值, i=0, 1,…,15 .
多项式求逆函数的执行步骤如下,首先初始化 Phi,F,V,S ,执行15次循环,对应算法15的6~19行. 计算常数项 {\Phi }_{2n-1}\left(0\right) 在 {\mathbb{Z}}_{\mathrm{q}} 中的逆元. 执行8次循环,对应算法15的第20行. 这里本文需要计算32b数乘积,本文采用vpmovzxdq指令将每64b的高32b和低32b扩展成64b分别放置在2个向量寄存器中进行计算. 最后需要返回是否可逆的真值,我们将向量寄存器中每个数据单元中 Delta 的值相加后进行返回.
4.5 多项式序列化与反序列化AVX2实现
本文将多项式系数序列化为字节流和字节流反序列回多项式系数这2个过程进行AVX2向量化实现.CTRU的公钥系数实际占12b,因此2个多项式系数可以生成3个字节(24b). C实现的思路是将第1个系数的前8b作为第1个字节,再将第1个系数剩下的4位和第1个系数的前4位拼接作为第2个字节,最后将第2个系数剩下的8b作为第3个字节进行拼接. 如图13所示. 由于私钥 sk 对应的多项式 f 系数存在负数的情况,这里本文先对所有系数通过加上界的处理转换到正数域,再进行拼接,同时反序列化的时候需要从正数域再恢复到原值. 密文序列化与反序列化的思路和公钥序列化的思路类似.
在AVX2优化实现中思路有所不同,因为向量寄存器中是按相同数据单元进行对应操作,为了方便向量化操作,本文将系数拼接到对应16b数据单元中. 最后再按8b数据单元读回内存中,因此需要4个系数才能完整拼接3个16b数据单元,如图14所示,本文使用vpsllw,vpsrlw,vpor,vpand这些16b数据单元移位指令对数据进行移位拼接操作.
4.6 SHA-3哈希函数AVX2优化实现
本文对文献[8]中CTRU AVX2实现的SHA-3模块进行了优化. 根据SHA-3标准文档FIPS202[43]可知,SHA-3模块是由一个海绵结构组成,海绵结构包含吸收和挤压2个部分,吸收和挤压通过调用Keccak轮函数与输入的字节流进行分组异或得到最终输出结果. Keccak轮函数采用一个1 600b长的状态,在实现中将这个1 600b状态存储为一个64b的长度为25的1维数组. 一个AVX2向量寄存器可以存储4个64b状态,因此7个向量寄存器就能存放完25个64b状态,其中1个向量寄存器通过广播指令vpbroadcastq存放第1个状态. 通过这样的方式,本文可以一次性对4个状态进行置乱和异或操作,相较于C参考实现可以达到4并行.
5. 实验结果与性能比较
5.1 测试环境
本文给出了CTRU密钥封装方案3个方案的AVX2优化实现,本节介绍了3个方案优化实现的实现性能和测试分析. 本文的AVX2测试硬件环境为3.7 GHz的Intel(R) Core(TM) i9-10900X CPU和16GB内存,编译器为gcc-9.4.0,本文使用的编译参数为-Wall -Wextra -Wpedantic -Wmissing-prototypes -Wredundant-decls -Wshadow -Wpointer-arith -mavx2 -mbmi2 -mpopcnt -maes -march=native -mtune=native -O3 -fomit-frame-pointer.本文关闭了处理器的睿频加速技术(turbo boost)和超线程技术(hyper-threading).
本文采用CPU周期数作为衡量算法效率的标准,AVX2代码测试的CPU周期数为对应的函数执行
100000 次后所得结果的中位数.5.2 多项式乘法实现性能比较
表12展示了参考实现和AVX2实现的正向NTT、逆向NTT以及多项式点乘的速度对比,从结果可以看出,3个方案下的多项式环正向NTT、逆向NTT以及点乘的AVX2优化实现在性能上有显著提升,3个方案的多项式乘法优化率都达到了95%以上,我们还另外测试了Kyber多项式乘法的性能,从表12中数据可以看出我们AVX2实现的优化率与Kyber接近.NTT实现AVX2优化性能达到显著提升主要原因是NTT比较向量化友好,并且本文通过层融合技术减少了内存和向量寄存器之间的访存次数和合理安排寄存器达到寄存器利用率最大化.
表 12 多项式乘法实现性能比较Table 12. Performance Comparisons of Polynomial Multiplication Implementation方案 测试函数 参考实现/
时钟周期AVX2实现/
时钟周期优化率/% CTRU-512 正向NTT 11 120 502 95 逆向NTT 16 736 530 97 点乘 9 754 462 95 CTRU-768 正向NTT 24 520 882 96 逆向NTT 31 592 950 97 点乘 9 998 278 97 CTRU- 1024 正向NTT 58 442 1 010 98 逆向NTT 55 862 1 070 98 点乘 43 658 1 852 96 Kyber[44] 正向NTT 7 042 216 97 逆向NTT 11 002 232 98 点乘 3 708 84 98 5.3 多项式求逆实现性能比较
多项式求逆是CTRU中最为耗时的运算之一,其需要返回是否可以求逆的真值以及存在的一些模幂运算给向量化实现带来了困难. 如表13所示,因为本文并没与对CTRU-768方案进行多项式求逆实现上的优化,因此我们这里只给出CTRU-512和CTRU-
1024 方案的优化率. 可以看到CTRU-512方案的优化率达到了96.9%. 而CTRU-1024 方案的优化率只有46.8%,因为此时采用的Bernstein多项式快速求逆方法在向量化实现中存在一些困难,需要不断申请栈空间和访问内存,因此Bernstein多项式快速求逆AVX2实现的加速比最低.表 13 多项式求逆实现性能比较Table 13. Performance Comparisons of Polynomial inversion Implementation方案 参考实现/时钟周期 AVX2实现/时钟周期 优化率/% CTRU-512 92 026 2 898 96.9 CTRU- 1024 398 542 212 048 46.8 5.4 多项式序列化与反序列化实现性能比较
在3个方案下的多项式序列化和反序列化AVX2实现性能比较如表11所示,从表14中可以看出,AVX2实现对多项式序列化与反序列化带来的性能提升较为显著,总体带来了56%~90%以上的性能提升. 性能提升主要来自于我们将比特拼接操作利用AVX2指令进行向量化.
表 14 多项式序列化与反序列化实现性能比较Table 14. Performance Comparison of Polynomial Serialization and Deserialization Implementation测试函数 参考实现/时钟周期 AVX2实现/时钟周期 优化率/% 序列化公钥 396 32 91.9 反序列化公钥 506 38 92.5 序列化私钥 82 36 56.1 反序列化私钥 590 38 93.6 序列化密文 796 34 95.7 反序列化密文 934 42 95.5 5.5 CTRU方案整体实现性能比较
5.5.1 CTRU-768 AVX2优化比较
本文对已有的CTRU-768方案AVX2实现进行了一系列优化,从表15可以看出,本文的实现相较于目前已有的最新实现有着8%~11%的性能提升,性能提升主要来自于哈希SHA-3模块的AVX2优化实现.
表 15 CTRU方案实现性能比较Table 15. Performance Comparisons of CTRU Scheme Implementation测试函数 文献[8]/时钟周期 本文实现/时钟周期 优化率/% 密钥生成 17 376 18 958 8 密钥封装 12 146 13 650 11 密钥解封装 48 292 53 792 10 5.5.2 AVX2实现与参考实现比较
为了测试本文的AVX2优化实现技术对CTRU整体方案性能带来的提升,本文测试了3个方案的AVX2实现时钟周期和参考实现时钟周期. 从表16和表17中可以看出,CTRU-512和CTRU-768方案的密钥生成提速较为显著,均有90%的提速. 但CTRU-
1024 方案的密钥生成仅有56%的提速,主要原因还是由于CTRU-1024 方案无法使用多项式矩阵快速计算多项式求逆,而Bernstein快速求逆算法又不是向量化友好的,因此并不能很好的利用AVX2的并行度加速. 除此之外密钥封装和解封装都有70%~90%以上的提速. 除此之外由于CTRU方案基于RLWE问题,CNTR方案基于RLWR问题,在密钥封装过程中需要多生成1个噪音多项式,因此CTRU方案相较于CNTR方案的封装过程会更加耗时. 但在CNTR-512方案下由于CNTR的多项式噪音向量分布为 ({B}_{5},{B}_{5}) ,相较于CTRU的噪音多项式分布 ({B}_{3},{B}_{3}) 需要更多的字节流用于采样生成多项式,哈希的时间也会增加,因此CNTR-512方案的密钥生成、密钥封装、密钥解封装要慢于CTRU-512方案.表 16 CTRU方案实现性能比较Table 16. Performance Comparisons of CTRU Scheme Implementation方案 测试函数 参考实现/
时钟周期AVX2实现/
时钟周期优化率/% CTRU-512 密钥生成 143 244 13 860 90 密钥封装 57 962 12 230 79 密钥解封装 131 964 37 010 72 CTRU-768 密钥生成 154 156 13 940 91 密钥封装 91 636 13 340 85 密钥解封装 218 654 49 182 78 CTRU- 1024 密钥生成 589 650 259 782 56 密钥封装 190 894 20 018 90 密钥解封装 396 076 70 948 82 表 17 CNTR方案实现性能比较Table 17. Performance Comparisons of CNTR Scheme Implementation方案 测试函数 参考实现/
时钟周期AVX2实现/
时钟周期优化率/% CTRU-512 密钥生成 150 294 18 350 88 密钥封装 63 210 16 278 74 密钥解封装 137 492 41 114 70 CTRU-768 密钥生成 157 326 17 376 89 密钥封装 89 372 12 146 86 密钥解封装 221 926 48 292 78 CTRU- 1024 密钥生成 589 054 263 832 55 密钥封装 189 334 17 312 91 密钥解封装 393 648 68 378 83 5.5.3 AVX2实现与Kyber方案比较
除了与参考实现进行比较,本文还测试了NIST标准化的Kyber密钥封装方案的AVX2实现,与我们的CTRU/CNTR方案进行比较. 因为文献[31]中没有提供前缀哈希FO转换的Kyber开源代码,因此我们采用的Kyber官方代码[44]在其上面进行修改. 测试数据如表18所示,在CTRU-512和CTRU-
1024 方案下,CTRU方案的密钥生成和密钥封装性能都要远优于Kyber的性能. 但在CTRU-1024 方案下,CTRU密钥封装方案比Kyber慢很多,主要原因还是在多项式求逆部分AVX2优化提速有限,导致整个密钥生成的提速不明显,而Kyber由于解封装不存在较为耗时且难以向量化的多项式求逆. 总体来说,CTRU方案相较于Kyber方案还是有一定性能优势.表 18 CTRU/CNTR/Kyber方案AVX2实现性能比较Table 18. Performance Comparisons of CNTR Scheme Implementation方案 测试函数 CTRU方案/时钟周期 CNTR方案/时钟周期 Kyber方案/时钟周期 CTRU-512 密钥生成 13 860 18 350 28 936 密钥封装 12 230 16 278 23 524 密钥解封装 37 010 41 114 31 688 CTRU-768 密钥生成 13 940 17 376 42 554 密钥封装 13 340 12 146 35 204 密钥解封装 49 182 48 292 49 456 CTRU- 1024 密钥生成 259 782 263 832 54 150 密钥封装 20 018 17 312 48 190 密钥解封装 70 948 68 378 69 918 6. 总 结
伴随着美国NIST后量子密码标准的推出,学术与工业界都在加快推进网络安全生态向后量子时代迁移.CTRU作为我国自主可控的抗量子密码算法,其标准化工作目前已经在我国密码标准委员会正式立项. 自主可控的后量子密码的高效安全实现对我国后量子迁移与应用尤为重要.
本文首次完成了CTRU-512和CTRU-
1024 方案的完整参考实现和AVX2优化实现,并对已有的CTRU-768方案的AVX2实现进行了优化,较为系统地针对不同实现给出性能加速策略和优化技巧,提出了Bernstein快速求逆AVX2优化实现方案. 实验结果表明,本文给出的高效实现创造了CTRU方案在Intel CPU平台上新的性能记录. 性能提升的主要来源于AVX2向量化实现技巧如NTT向量化实现的层融合和系数置乱,指令调度,多项式求逆耗时模块的优化实现. 这为CTRU方案的实际应用打下重要基础,具有重要的参考和应用价值. 本文提出的基于索引的延迟约减以及对CTRU-1024 方案的多项式求逆采用Bernstein快速求逆方法是在算法层面的优化,不依赖于特定指令集,可以应用在ARM处理器等嵌入式平台.作者贡献声明:郑婕妤负责实现方案设计和论文撰写和代码测试;宋振宇负责多项式求逆代码实现;朱浩亮负责论文检查和文献调研;赵运磊、林璟锵、范晶负责论文检查和修改.
-
表 1 iRDMA协议字段
Table 1 iRDMA Protocol Fields
字段 字段含义 SE 重发请求指示 M 路径迁移指示 Pad 补充字节数,载荷4B对齐发送 TVer BTH的版本信息,默认为0 A 需要返回应答包指示 R_key 与存储的地址空间对应 WPSN 写(读)操作包序列号 WMSN 写(读)消息序列号 Syndrome 操作成功指示 表 2 初始化接口
Table 2 Initialize Interface
初始化API 功能 Icloud_init 实现FPGA云编程接口,初始化FPGA
云平台软件内部资源Icloud_close 云编程接口关闭,释放FPGA云平台软件内部资源 icfReadReg_n 读取FPGA多个寄存器内容 icWriteReg_n 向FPGA寄存器填写内容 表 3 配置接口
Table 3 Configuration Interface
配置 kernel API 功能 ifcConfigKernel 完成将每个输入数据的加1操作后将数据
放到定义参数result addr中ifeStartKernel 启动对应设备ID的FPGA的kernel功能 icfStartLoadPR 通过配置相关寄存器,进行PR(partial reconfiguration)加载 icfFlashPcieMem 烧写指定的格的文件到FPGA设备 icfLoadApplicationlmage 加载FPGA applicaton文件 表 4 数据接口
Table 4 Data Interface
数据传输API 功能 ifcHostAddrData 实现FPGA云编程接口,初始化FPGA
云平台软件内部资源ifcFpgaSwFpgaDatauint8 通过RDMA网络接口在2块 FPGA
设备之间传输数据ifcConfiglkl 将FPGA的kerel计算加1的数据传输到目的FPGA的数据存放地 表 5 流完成时间
Table 5 Time for Flow Completion
μs 测试内容 FCT write_4B 10 read_4B 4 write_1KB 10 read_1KB 4 write_1MB 103 read_1MB 94 write_10MB 939 read_10MB 905 write_100MB 9126 read_100MB 8962 write_1GB 91121 read_1GB 89198 表 6 硬件配置
Table 6 Hardware Configuration
硬件 配置 10 Gbps光纤交换机 Ruijie S6220 CPU型号 Intel Xeon E5-2650 v4@2.4GHz FPGA加速卡 型号:Inspur F10A FPGA芯片:10AX115H3F34E2SG 连接方式:SFP+ 10Gb Ethernet /PCIe Gen3 x8 -
[1] Zhu Zongwei, Zhang Junneng, Zhao Jinjin, et al. A hardware and software task-scheduling framework based on CPU+FPGA heterogeneous architecture in edge computing[J]. IEEE Access, 2019, 7: 148975−148988 doi: 10.1109/ACCESS.2019.2943179
[2] Choi Y, Cong J, Fang Zhenman, et al. A quantitative analysis on microarchitectures of modern CPU-FPGA platforms[C/OL]//Proc of the 53rd Annual Design Automation Conf. New York: ACM, 2016[2024-07-09]. https://dl.acm.org/doi/abs/10.1145/2897937.2897972
[3] Man Xingchen, Zhu Jianfeng, Song Guihuan, et al. CaSMap: Agile mapper for reconfigurable spatial architectures by automatically clustering intermediate representations and scattering mapping process[C]//Proc of the 49th Annual Int Symp on Computer Architecture. New York: ACM, 2022: 259−273
[4] 齐乐,常轶松,陈欲晓,等. 基于SoC-FPGA的RISC-V处理器软硬件系统级平台[J]. 计算机研究与发展,2023,60(6):1204−1215 doi: 10.7544/issn1000-1239.202330060 Qi Le, Chang Yisong, Chen Yuxiao, et al. A system-level platform with SoC-FPGA for RISC-V hardware-software integration[J]. Journal of Computer Research and Development, 2023, 60(6): 1204−1215 (in Chinese) doi: 10.7544/issn1000-1239.202330060
[5] Zha Yue, Li Jing. Virtualizing FPGAs in the cloud[C]//Proc of the 25th Int Conf on Architectural Support for Programming Languages and Operating Systems. New York: ACM, 2020: 845−858
[6] Chung E, Fowers J, Ovtcharov K, et al. Serving DNNs in real time at datacenter scale with project brainwave[J]. IEEE Micro, 2018, 38(2): 8−20 doi: 10.1109/MM.2018.022071131
[7] Suda N, Chandra V, Dasika G, et al. Throughput-optimized OpenCL-based FPGA accelerator for large-scale convolutional neural networks[C]//Proc of the 2016 ACM/SIGDA Int Symp Field-Programmable Gate Arrays. New York: ACM, 2016: 16−25
[8] Zhang Jialiang, Li Jing. Improving the performance of OpenCL-based FPGA accelerator for convolutional neural network[C]//Proc of the 2017 ACM/SIGDA Int Symp on Field-Programmable Gate Arrays. New York: ACM, 2017: 25−34
[9] Tine B, Yalamarthy K P, Elsabbagh F, et al. Vortex: Extending the RISC-V ISA for GPGPU and 3D-graphics[C]//Proc of the 54th Annual IEEE/ACM Int Symp on Microarchitecture. Piscataway, NJ: IEEE, 2021: 754−766
[10] Caulfield A M, Chung E S, Putnam A, et al. A cloud-scale acceleration architecture[C/OL]//Proc of the 49th Annual IEEE/ACM Int Symp on microarchitecture (MICRO). Piscataway, NJ: IEEE, 2016[2024-07-09]. https://ieeexplore.ieee.org/abstract/document/7783710
[11] Amazon Web Services EC2. FPGA hardware and software development kit[EB/OL]. [2023-01-28]. https://github.com/aws/aws-fpga
[12] Tarafdar N, Thomas L, Fukuda E, et al. Enabling flexible network FPGA clusters in a heterogeneous cloud data center[C]//Proc of the 2017 ACM/SIGDA Int Symp on Field-Programmable Gate Arrays. New York: ACM, 2017: 237−246
[13] Shu Ran, Cheng Peng, Chen Guo, et al. Direct universal access: Making data center resources available to FPGA[C]//Proc of the 16th USENIX Symp on Networked Systems Design and Implementation (NSDI 19). Berkeley, CA: USENIX Association, 2019: 127−140
[14] Yu Xiaoyu, Wang Yuwei, Miao Jie, et al. A data-center FPGA acceleration platform for convolutional neural networks[C]//Proc of the 29th Int Conf on Field Programmable Logic and Applications (FPL). Piscataway, NJ: IEEE, 2019: 151−158
[15] Choi Y K, Jason C, Fang Zheman, et al. In-depth analysis on microarchitectures of modern heterogeneous CPU-FPGA platforms[J]. ACM Transactions on Reconfigurable Technology and Systems, 2019, 12(1): 1−20
[16] Fleming K, Adler M. The LEAP FPGA Operating System[M]//FPGAs for Software Programmers. Berlin: Springer, 2016: 245−258
[17] Khawaja A, Landgraf J, Prakash R, et al. Sharing, protection, and compatibility for reconfigurable fabric with AmorphOS[C]//Proc of the 13th USENIX Symp on Operating Systems Design and Implementation (OSDI'18). Berkeley, CA: USENIX Association, 2018: 107−127
[18] Baxter R, Booth S, Bull M, et al. Maxwell-a 64 FPGA supercomputer[C]//Proc of the 2nd NASA/ESA Conf on Adaptive Hardware and Systems (AHS 2007). Piscataway, NJ: IEEE, 2007: 287−294
[19] Jeremy F, Kalin O, Michael P, et al. A configurable cloud-scale DNN processor for real-time AI[C/OL]//Proc of the 45th Annual Int Symp on Computer Architecture. Piscataway, NJ: IEEE, 2018[2024-07-09]. https://ieeexplore.ieee.org/abstract/document/8416814
[20] Ouyang J, Shiding L, Qi Wei, et al. SDA: Software-defined accelerator for large-scale DNN systems[C]//Proc of the 26th IEEE Hot Chips Symp (HCS). Piscataway, NJ: IEEE, 2014: 10–12
[21] Vesper M, Koch D, Vipin K, et al. JetStream: An open-source high-performance PCI express 3 streaming library for FPGA-to-Host and FPGA-to-FPGA communication[C/OL]//Proc of the 26th Int Conf on Field Programmable Logic and Applications (FPL). Piscataway, NJ: IEEE, 2016[2024-07-09]. https://ieeexplore.ieee.org/abstract/document/7577334
[22] Jacobsen, M, Richmond, D, Hogains, M, et al. RIFFA 2.1: A reusable integration framework for FPGA accelerators[J]. ACM Transactions on Reconfigurable Technology and Systems, 2015, 8(4): 1−23
[23] Zeke W, Zhang Shuhao, He Bingsheng, et al. Melia: A MapReduce framework on OpenCL-Based FPGAs[J]. IEEE Transactions on Parallel and Distributed Systems, 2016, 27(12): 3547−3560 doi: 10.1109/TPDS.2016.2537805
[24] Sharma D D, Blankenship R, Berger D S. An introduction to the compute express link (CXL) interconnect[J]. arXiv preprint, arXiv: 2306.11227, 2023
[25] Wang Fu, Yan Fulong, Xue Xuwei, et al. Traffic load balancing based on probabilistic routing in data center networks[C/OL]//Proc of the Int Conf on Optical Network Design and Modeling (ONDM). Piscataway, NJ: IEEE, 2020[2024-07-09]. https://ieeexplore.ieee.org/abstract/document/9133002
[26] Mittal R, Shpiner A, Panda A, et al. Revisiting network support for RDMA[C]//Proc of the 2018 Conf of the ACM Special Interest Group on Data Communication. New York: ACM, 2018: 313−326
[27] Biookaghazadeh S, Zhao Ming, Ren Fengbo. Are FPGAs suitable for edge computing[C]//Proc of the USENIX Workshop on Hot Topics in Edge Computing (HotEdge'18). Berkeley, CA: USENIX Association, 2018. https://www.usenix.org/conference/hotedge18/presentation/biookaghazadeh
[28] Belabed T, Coutinho M G F, Fernandes M A C, et al. User driven FPGA-based design automated framework of deep neural networks for low-power low-cost edge computing[J]. IEEE Access, 2021, 9: 89162−89180 doi: 10.1109/ACCESS.2021.3090196
[29] Ross S M. Introduction to Probability Models[M]. Amsterdam, Netherlands: Elsevier, 2014
[30] 段田田,郭仪,李博,等. PieBridge:一种按需可扩展的跨链架构[J]. 计算机研究与发展,2023,60(11):2520−2533 doi: 10.7544/issn1000-1239.202230284 Duan Tiantian, Guo Y, Li Bo, et al. PieBridge: An on-demand scalable cross-chain architecture[J]. Journal of Computer Research and Development, 2023, 60(11): 2520−2533 (in Chinese) doi: 10.7544/issn1000-1239.202230284
[31] 张帆,胡成臣. 一种适用突发流量的数据中心网络流调度策略[J]. 软件学报,2018,28(s2):81−89 Zhang Fan, Hu Chengchen. Flow scheduling policy for burst traffic in data center networks[J]. Journal of Software, 2018, 28(s2): 81−89 (in Chinese)