传统的公钥密码大多基于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]. 我国学者基于${\mathrm{E}}_{8}$格编码^[7]和NTRU格提出了CTRU和CNTR密钥封装算法^[8]，简记为CTRU，其中，CNTR是CTRU的一个简化版本. CTRU和CNTR本质上还是NTRU机制，其密钥生成和解密过程与传统的NTRU机制保持一致. 其优点是采用了尺度化${\mathrm{E}}_{8}$格编码，提高纠错能力，能够容忍更大的秘密范围和更小的模数，并支持密文压缩，从而在安全性和带宽上相比传统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方案的${\mathbb{Z}}_{q}\left[x\right]/\left({x}^{8}-{\zeta }^{r\left(i\right)}\right)$多项式环求逆，使用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$为正整数. 定义集合${\mathbb{Z}}_{q}=\mathbb{Z}/q\mathbb{Z}=\{0,\;1,\;…, q-1\}$，${\mathbb{R}}_{q}=\mathbb{R}/q$. 定义多项式环$\mathcal{R}=\mathbb{Z}\left[x\right]/({x}^{n}-{x}^{n/2}+1)$，${\mathcal{R}}_{q}={\mathbb{Z}}_{q}\left[x\right]/({x}^{n}-{x}^{n/2}+1)$，其中每个多项式系数均在环${\mathbb{Z}}_{q}$中. 令$f$为${\mathcal{R}}_{q}$中的多项式，可表示为$f={\displaystyle\sum\limits^{n-1}_{i=0}}{f}_{i}{x}^{i}$. 多项式$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

下载: 导出CSV 

| 显示表格

表  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

下载: 导出CSV 

| 显示表格

算法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变换快速计算：

其中${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有如下同构：

其中$\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技巧计算如下：

我们可以将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数

下载: 导出CSV 

| 显示表格

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所示：

图  1  $n=512$多项式同构树形图

Figure  1.  The tree map of $n=512$ polynomial isomorphism

下载: 全尺寸图片 幻灯片

图  2  $n=768$多项式同构树形图

Figure  2.  The tree map of $n=768$ polynomial isomorphism

下载: 全尺寸图片 幻灯片

图  3  $n=1 024$多项式同构树形图

Figure  3.  The tree map of $n=1 024$ polynomial isomorphism

下载: 全尺寸图片 幻灯片

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 CTRU

$n$	多项式环	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)$

下载: 导出CSV 

| 显示表格

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)$

下载: 导出CSV 

| 显示表格

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 455$	$32$
第$5$层	$0\to 7,128\to 135$ $136\to 143$ $256\to 263,384\to 391$ $392\to 399$	$48$
第$6$层	$8\to 15,264\to 271$ $16\to 31,272\to 287$	$48$

下载: 导出CSV 

| 显示表格

表  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 682$	$48$
第$5$层	$0\to 11,384\to 395$ $204\to 215,588\to 599$ $192\to 203,576\to 587$	$72$
第$6$层	$12\to 23,396\to 407$ $24\to 47,408\to 431$	$72$

下载: 导出CSV 

| 显示表格

表  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 911$	$64$
第$5$层	$0\to 15,512\to 527$ $272\to 287,784\to 799$ $256\to 271,768\to 783$	$96$
第$6$层	$16\to 31,528\to 543$ $32\to 63,544\to 575$	$96$

下载: 导出CSV 

| 显示表格

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$阶旋转矩阵：

计算$f$的逆可以转化为计算旋转矩阵的逆. 求旋转矩阵的逆需要求旋转矩阵行列式的逆，若行列式的逆不存在，则逆矩阵也不存在. 令$d$为系数矩阵的行列式，根据Cramer法则，存在唯一的${f}^{-1}，$其中每个系数可以通过如下公式计算：

${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

下载: 导出CSV 

| 显示表格

从表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

下载: 导出CSV 

| 显示表格

其中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  正向NTT载入系数顺序

Figure  4.  Load coefficient order of forward NTT

下载: 全尺寸图片 幻灯片

按图4顺序载入系数，共需要8个向量寄存器；并且环分解操作和基-2 NTT的第1层到第2层的系数对都在寄存器中，因此可以进行融合操作，以减少访存次数，本文采用调用宏加传入偏移地址的方法实现多项式所有系数的运算.

环分解操作及基-2 NTT 第1层到第2层. 对于正向基-2 NTT的前2层，本文一次性载入128个系数到8个向量寄存器，载入系数顺序如图5所示.

图  5  正向NTT第1层载入系数顺序

Figure  5.  Load coefficient order of forward NTT level 1

下载: 全尺寸图片 幻灯片

前2层系数对蝴蝶操作间隔分别为128和64，按照图5所示顺序载入系数刚好可以完成前2层系数对间隔的蝴蝶操作. 基-2 NTT第3层到第6层. NTT后4层载入系数顺序如图6所示. 本文使用8个向量寄存器一次性载入128个系数.

图  6  正向NTT第4层载入系数顺序

Figure  6.  Load coefficient order of forward NTT level 4

下载: 全尺寸图片 幻灯片

对于系数置乱，后4层的间隔分别为32，16，8，4；由于只能载入连续的16个系数到向量寄存器，因此上述存放顺序无法完成相应系数间隔的蝴蝶操作，需要调整向量寄存器中系数顺序，这里本文称调整向量寄存器中系数顺序的操作为Shuffle.对目前已经载入向量寄存器的系数做Shuffle操作，本文需要对每8，4，2，1个系数进行Shuffle，这里本文称为Shuffle8，Shuffle4，Shuffle2，Shuffle1.后三者的shuffle操作流程如图7所示.

图  7  基-2 NTT后3层系数置乱流程图

Figure  7.  Flow chart of 3-layer coefficient scrambling after Radix-2NTT

下载: 全尺寸图片 幻灯片

Shuffle8的主要操作如图8所示，取0～15的前8个系数和32～47的前8个系数拼起来放到1个向量寄存器中，后8个系数8～15和40～47放到1个寄存器中.Shuffle8后寄存器中存放系数顺序如图9所示，两两寄存器之间的连线和蝴蝶表示这2对寄存器之间完成1对蝴蝶操作，对应第5层需要进行间隔8的蝴蝶操作.

图  8  调整8个系数顺序和调整4个系数顺序

Figure  8.  Adjusting the order of 8 coefficients and adjusting the order of 4 coefficients

下载: 全尺寸图片 幻灯片

图  9  Shuffle8后向量寄存器存放顺序

Figure  9.  Storage order of vector registers after Shuffle8

下载: 全尺寸图片 幻灯片

对于第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$的实现思路类似.

图  10  调整2个系数顺序

Figure  10.  Adjusting the order of the two coefficients

下载: 全尺寸图片 幻灯片

图  11  调整一个系数顺序

Figure  11.  Adjusting the order of a coefficient

下载: 全尺寸图片 幻灯片

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

下载: 导出CSV 

| 显示表格

将堆栈中数据按照图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$.

图  12  变量存储格式

Figure  12.  Variants storage format

下载: 全尺寸图片 幻灯片

多项式求逆函数的执行步骤如下，首先初始化$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$系数存在负数的情况，这里本文先对所有系数通过加上界的处理转换到正数域，再进行拼接，同时反序列化的时候需要从正数域再恢复到原值. 密文序列化与反序列化的思路和公钥序列化的思路类似.

图  13  公钥序列化

Figure  13.  Public key serialization

下载: 全尺寸图片 幻灯片

在AVX2优化实现中思路有所不同，因为向量寄存器中是按相同数据单元进行对应操作，为了方便向量化操作，本文将系数拼接到对应16b数据单元中. 最后再按8b数据单元读回内存中，因此需要4个系数才能完整拼接3个16b数据单元，如图14所示，本文使用vpsllw，vpsrlw，vpor，vpand这些16b数据单元移位指令对数据进行移位拼接操作.

图  14  公钥序列化AVX2实现

Figure  14.  Public key serialization AVX2 implementation

下载: 全尺寸图片 幻灯片

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

下载: 导出CSV 

| 显示表格

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

下载: 导出CSV 

| 显示表格

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

下载: 导出CSV 

| 显示表格

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

下载: 导出CSV 

| 显示表格

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

下载: 导出CSV 

| 显示表格

表  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

下载: 导出CSV 

| 显示表格

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

下载: 导出CSV 

| 显示表格

6.   总　　结

伴随着美国NIST后量子密码标准的推出，学术与工业界都在加快推进网络安全生态向后量子时代迁移.CTRU作为我国自主可控的抗量子密码算法，其标准化工作目前已经在我国密码标准委员会正式立项. 自主可控的后量子密码的高效安全实现对我国后量子迁移与应用尤为重要.

本文首次完成了CTRU-512和CTRU-1024方案的完整参考实现和AVX2优化实现，并对已有的CTRU-768方案的AVX2实现进行了优化，较为系统地针对不同实现给出性能加速策略和优化技巧，提出了Bernstein快速求逆AVX2优化实现方案. 实验结果表明，本文给出的高效实现创造了CTRU方案在Intel CPU平台上新的性能记录. 性能提升的主要来源于AVX2向量化实现技巧如NTT向量化实现的层融合和系数置乱，指令调度，多项式求逆耗时模块的优化实现. 这为CTRU方案的实际应用打下重要基础，具有重要的参考和应用价值. 本文提出的基于索引的延迟约减以及对CTRU-1024方案的多项式求逆采用Bernstein快速求逆方法是在算法层面的优化，不依赖于特定指令集，可以应用在ARM处理器等嵌入式平台.

作者贡献声明：郑婕妤负责实现方案设计和论文撰写和代码测试；宋振宇负责多项式求逆代码实现；朱浩亮负责论文检查和文献调研；赵运磊、林璟锵、范晶负责论文检查和修改.