基于FPGA的SM9数字签名算法的快速实现

Fast Implementation of SM9 Digital Signature Algorithm Based on FPGA

摘要: SM9算法是由我国推出的基于双线性对的标识密码算法，其中使用的R-ate双线性对在计算效率上优于Weil对和Tate对，但在实际应用中仍有不足。针对上述问题，为了提升SM9中双线性对的计算效率，提出了一种在FPGA平台上使用滑动窗口改进R-ate对运算的方法。该方法以蒙哥马利模乘为基础改进扩域模逆运算并通过滑动窗口算法优化R-ate双线性对运算中的Miller循环部分。同时根据蒙哥马利模乘特性减少椭圆曲线上点运算计算量提升签名算法运算效率。仿真结果表明，采用该方法可以将SM9中R-ate对的运算效率提升约18.46%，对SM9签名算法效率提升约为13.55%。

Abstract: SM9 algorithm is an identity-based cryptographic algorithm based on bilinear pairing introduced by our country, in which the R-ate pairing used is better than Weil pairing and Tate pairing in terms of computational efficiency, but it still has deficiencies in practical applications. To improve the computational efficiency of the bilinear pairing in SM9, a method was proposed to improve the R-ate pairing operation by using a sliding window on the FPGA platform. The method was based on the Montgomery modular multiplication to enhance the modular inverse operation in the twelfth extension field and optimize the Miller loop in the R-ate pairing operation through the sliding window algorithm. At the same time, the computational efficiency of the signature algorithm was improved by reducing the number of point operations on elliptic curves according to the Montgomery modular multiplication property. Simulation results showed that this method could improve the computational efficiency of R-ate pairing in SM9 by approximately 18.46%, and the efficiency of the SM9 signature algorithm by approximately 13.55%.

1. 引言

在现代密码学中，双线性对(Bilinear Pairing)技术因其广泛的应用场景和独特的数学特性，成为众多密码学协议的核心组成部分。双线性对的应用涵盖了身份基加密[1] (Identity-based cryptography, IBC)、签名方案[2]-[6]和各种基于配对的密码协议[7] [8]。其中，R-ate对[9]作为一种优化的双线性配对，逐渐引起了研究者们的广泛关注。R-ate对相较于传统的Weil对[10]和Tate对[11]，具有更高的计算效率和更优的性能表现。这使得R-ate对在高性能需求的密码学应用中显示出巨大潜力。然而，尽管R-ate对已经在理论和实现方面取得了显著进展，其在实际应用中的计算效率仍不能完全满足实际需求。因此，如何进一步优化双线性对的计算过程，成为当前研究的重要方向。

研究者们为了提升双线性对的计算效率，探索并提出了多种创新策略。其中包括：通过将最终模幂的指数进行分解[12] [13]，利用模乘及Frobenius映射进行计算，降低了计算复杂度。通过采用不同坐标系[14] (如仿射坐标系和雅克比坐标系)，改进椭圆曲线上点的运算以优化配对运算。提出一种利用分母约化技术[15]来提升运算效率的方法。从算法角度考虑，将双线性对运算中参数采用非相邻表示型[16] (non-adjacent form, NAF)算法展开，降低参数汉明重量以减少计算量。此外，有学者突破传统计算框架，设计出一种椭圆网算法[17]替代Miller算法[18]用于双线性对的计算。也有学者通过将双基数链[19]应用到双线性对的计算中，有效降低了Tate对运算的计算复杂度。

自2016年国家密码管理局发布《SM9标识密码算法》[20]以来，为了促进SM9算法的普及与实际应用，国内学者提出了多种针对SM9中R-ate对的运算优化方案。在SM9算法流程中，算法的安全性依赖于R-ate双线性对的计算，相应的R-ate对的计算量占据整个算法中很大一部分。因此许多学者都在探索如何能够在不损失双线性对安全性的同时提升其计算的效率。为此有学者提出在算法底层通过改变R-ate对计算中同构映射的操作顺序[21]，将大部分逆元运算从大特征域转换到小特征域，减少了逆元求解的计算消耗。对于双线性对映射中群的运算提出了一系列有效的优化方法[22]以加速双线性算法的执行，包括NAF加速、蒙哥马利域下运算和十二次扩域稀疏乘法等。而文献[23]设计了一种十二次扩域的等价转换方法，提出了一种针对椭圆曲线点乘的稀疏乘法，从而减少了有限域乘法的数量。文献[24]则对运算中Miller循环、最终模幂等步骤进行优化并在硬件上进行实现。尽管已有多种针对SM9算法的优化方案[21]-[25]，但距离算法的普及使用仍有一定差距，仍需要更多工作以实现SM9算法的效率提升。

本文以Jacobian加重射影坐标系为基础，底层运算采用Montgomery模乘以及辗转相除法，利用滑动窗口算法改进R-ate对运算中Miller循环部分计算，有效降低R-ate对运算过程中线函数运算以及椭圆曲线上点加运算的数量。在此基础上利用Montgomery特性优化椭圆曲线上点乘计算，减少计算开销，最终得到一种在FPGA平台上快速计算R-ate对以及SM9数字签名的方法，并通过仿真实验证明其正确性。最后，给出算法的计算量分析以及效率评估。

2. 基础知识

2.1. Miller算法

目前计算配对的算法主要有Miller算法以及椭圆网算法。椭圆网算法优势在于计算过程中无需求逆操作，其计算同样适用于Weil对、Tate对，但是运算效率相对于Miller算法仍然偏低。并且Miller算法的使用范围要比椭圆网算法更加广泛。Miller函数f_T,Q满足的迭代关系如下：

$f_{m+n,Q}=f_{m,Q}{f}_{n,Q}\frac{g_{\left[m\right]Q,\left[n\right]Q}}{g_{\left[m+n\right]Q}}$ (1)

其中$g_{\left[m\right]Q,\left[n\right]Q}$ 表示过点[m]Q和[n]Q的直线方程，$g_{\left[m+n\right]Q}$ 表示过点[m+n]Q的垂线方程。Miller算法则是采用迭代关系中n = m以及n = 1的两种特殊情况，将Miller函数按照二进制展开的形式求解。

2.2. R-ate对

R-ate对作为SM9中最常用的双线性对，可以看作是在Tate对基础上优化幂次的配对。在计算过程中大多采用Miller算法进行运算，R-ate对具体计算如算法1所示。

算法1：R-ate对的计算

输入：$P\in E\left(F_q\right)\left[r\right],Q\in E^{\prime }\left(F_{q^2}\right)\left[r\right],a=6t+2$

输出：$R_a\left(Q,P\right)$

1) 设$a={\displaystyle {\sum }_{i=0}^{L-1}{a}_i}{2}^i,a_{L-1}=1$

2) 置f = 1，T = Q

3) For i from L-2 to 0 do：

4) 计算$f=f^2\cdot g_{T,T}\left(P\right)$ ，$T=\left[2\right]T$

5) if $a_i=1$

6) 计算$f=f\cdot g_{T,Q}\left(P\right)$ ，$T=T+Q$

7) end if

8) end for

9) 计算$Q_1={\pi }_q\left(Q\right),Q_2={\pi }_{q^2}\left(Q\right)$

10) 计算$f=f\cdot g_{T,Q_1}\left(P\right),T=T+Q_1$

11) 计算$f=f\cdot g_{T,-Q_2}\left(P\right),T=T-Q_2$

12) 计算$f=f^{\left(q^{12}-1\right)/r}$

13) 输出f

2.3. SM9数字签名算法

根据GM/T 0044.2-2016中SM9签名算法如算法2所示，其中待签名的消息为比特串M，消息M的数字签名为(h, S)，签名主公钥为P_pub-s，签名密钥为ds_A，群G₁生成元为P₁，密码杂凑函数一般使用SM3算法。

算法2：SM9数字签名算法

输入：M, P_pub-s, ds_A

输出：(h, S)

1) 计算群G_T中的元素g = e (P₁, P_pub-s)

2) while l = 0 do

3) 产生随机数r ∈ [1, N − 1]

4) 计算群G_T中的元素w = g^r

5) 计算整数h = H₂ (M || w, N)

6) 计算整数l = (r − h) mod N

7) end while

8) 计算S = [l] ds_A

9) 确定数字签名(h, S)

3. 二次扩域的优化模逆算法

SM9算法作为基于椭圆曲线的密码算法，其底层运算主要为有限域模运算，在此基础上按照1-2-4-12的塔式扩张规则构成扩域运算。有限域模运算主要包括模乘运算、模加减运算以及模逆运算。其中模乘运算可以采用Montgomery模乘算法进行运算；模逆运算可以采用辗转相除法进行运算。扩域上的运算可以基于Karatsuba思想，将运算中的乘积项转换为运算中间项的加减组合，从而减少模乘运算的使用次数，变相转换为模加减运算，从而减少运算量。

在有限域运算的基础上构建扩域运算可以根据有限域算法特性以改进运算过程，以二次扩域模逆运算为例，假设存在二次扩域元素a = a₀ + a₁u，其中a₀, a₁∈F_q，其对应的逆为c = c₀ + c₁u，其中c₀, c₁∈F_q。根据塔式扩张运算规则c₀ = a₀/(a₀² + 2a₁²)，c₁ = −a₁/(a₀² + 2a₁²)。传统计算方式一般先计算分母，对分母求逆后乘上各元素对应分子得到结果。但由于Montgomery模乘算法相对比传统模乘运算计算过程可以分为模乘和形式转换两个部分，即对于待乘数a和b，第一步运算得到结果abR⁻¹，第二步将结果与R²再运算得到模乘结果。因此可以将二次扩域模逆运算中涉及到有限域模乘的运算步骤进行拆分以减少计算量。具体运算流程如算法3所示。

算法3：二次扩域模逆算法

输入：a = a₀ + a₁u，其中a₀, a₁∈F_q

输出：c = c₀ + c₁u

1) m = a2 0R⁻¹

2) n = a2 1R⁻¹

3) n = n + n

4) m = m + n

5) m = m⁻¹

6) c₀ = a₀∙m

7) c₁ = a₁∙m

8) c₁ = −c₁

通过算法3在计算分母过程中使用素数域元素直接进行Montgomery运算得到素数域数值带R⁻¹的结果，然后在分母求逆的过程将数值转换为Montgomery形式之后再与对应分子相乘得到素数域结果，这样利用Montgomery模乘的特性，可以减少常规运算流程中进行Montgomery形式转换的步骤。假设以A, M, I分别代表有限域模加减运算、模乘运算以及模逆运算所需计算量，其中M为两次Montgomery运算，则有限域及扩域各运算所需计算量如表1所示。其中二次扩域模逆运算的计算量仅需2M + 3A + I，相比于文献[23]中的2M + 2S + I + 3A以及文献[24]中的4M + 2I + 3A在二次扩域模逆上的运算效率都要更高。

Table 1. Computation analysis

表1. 计算量分析

	模加减	模乘	模平方	模逆
有限域	A	M	M	I
二次扩域	2A	3M + 3A	2M + 4A	2M + 3A + I
四次扩域	4A	9M + 21A	60M + 20A	14M + 19A + I
十二次扩域	~	54M + 190A	36M + 172A	113M + 196 + I

4. SM9数字签名算法快速实现

对于椭圆曲线上倍点运算，在传统仿射坐标系下运算简便但涉及模逆运算，因此计算效率低，而在雅克比坐标系下则可以节省模逆运算，将坐标点的求解转换为模加减运算以及模乘运算，计算效率更高。对于椭圆曲线上点加运算，在传统仿射坐标与雅克比坐标混合下计算效率要高于仅使用其中一种坐标系，这是由于传统仿射系仍旧需要进行模逆运算，而仅使用雅克比坐标系要比使用混合坐标系所需计算的参数更多。

除此之外，在双线性对的运算过程中椭圆曲线上点的运算总是与直线函数的运算成对出现，因此可以将椭圆曲线上点的运算与直线函数运算结合为一个整体进行计算。表2中总结了椭圆曲线上点的运算以及直线函数运算在优化前后所需运算数量。

Table 2. Computational analysis of points on elliptic curves

表2. 椭圆曲线上点计算量分析

	$F_{P^2}$ 平方	$F_{P^2}$ 乘法	$F_{P^2}$ 加减	$F_P$ 乘法
T = [2]T	2	6	4	10
T = T + Q	4	15	6	0
gT,T (P)	2	5	5	8
gT,Q (P)	2	12	5	4
优化后倍点	5	6	15	2
优化后点加	3	11	8	2

4.1. 扩域元素乘法设计

R-ate对的计算在Miller循环部分涉及大量直线函数与十二次扩域上元素乘法运算，但是直线函数运算结果在六个维度上的具体数值均为0。因此根据直线函数的稀疏特点，采用稀疏乘法算法以降低十二次扩域元素与直线函数之间计算的复杂度，减少不必要计算，可以提升计算效率。

根据稀疏乘法，假设存在十二次扩域元素f = f₀ + f₁∙w + f₂∙w₂以及一个直线函数运算结果g = l₀ + l₁∙v + l₂∙w₂其中f_i∈F_q4，l_i∈F_q2，元素的乘法结果如下：

$\begin{array}{c}f\cdot g=\left(f_0+f_1\cdot w+f_2\cdot w^2\right)\cdot \left(l_0+l_1\cdot v+l_2\cdot w^2\right)\\ =\left[f_0{l}_0+\left(f_0{l}_1+f_1{l}_2\right)v\right]+\left[f_1{l}_0+\left(f_1{l}_1+f_2{l}_2\right)v\right]w+\left(f_0{l}_2+f_2{l}_0+f_2{l}_{1}v\right)w^2\end{array}$ (2)

根据上式计算十二次扩域元素以及直线函数的乘法运算要相对原有十二次扩域乘法运算涉及维度更低，拆解到二次扩域上所需运算量更少，运算后结果同样为十二次扩域元素，不影响后续运算。相比于传统十二次扩域乘法运算需要55M + 190A，$f\cdot g$ 稀疏乘法运算可以将计算量节省至39M + 99A。

4.2. 基于滑动窗口的R-ate对快速实现

传统R-ate对计算采用Miller算法，即通过算法1中步骤4以及步骤6迭代计算实现，进而可以转换为一种类似标量乘的运算。然而在R-ate对计算过程中部分参数可以省略以节省运算时间，如采用NAF算法降低参数的汉明重量以减少点加的次数。假设参数为a，则其对应的NAF形式a'如下所示。

• $a=100100000000000000000000000000000000000010000101011101100100111110$

• $a^{\prime }=1001000000000000000000000000000000000000100010\overline{1}0\overline{1}000\overline{1}0\overline{1}001010000\overline{1}0$

在计算R-ate对的过程中遍历参数，当前遍历位为零则进行步骤4，为非零位时多进行一次步骤6。在此基础之上，可以采用滑动窗口算法改进R-ate对计算方法，通过预计算存储迭代过程中所涉及到的椭圆曲线上点的坐标以及对应的Miller函数值f，进一步降低Miller循环过程中所涉及点加运算的计算次数，从而提升R-ate对的计算效率，则可以将参数a改进为适用于滑动窗口的形式a''，以窗口长度r = 3为例，参数形式可以理解为如下所示：

• $a^{″}=100100000000000000000000000000000000000010000005000700300100007030$

由此遍历参数，将遍历非零位的计算改进为与对应数值所存储参数之间的计算，进而得到一种基于滑动窗口改进的计算R-ate双线性对的快速实现方法，具体实现算法如算法4所示。其中步骤24中模幂部分采用现有方法拆分为困难部分以及简单部分，通过加法链的方法来计算困难部分。

算法4：基于滑动窗口的R-ate对的计算

输入：$P\in E\left(F_q\right)\left[r\right],Q\in E^{\prime }\left(F_{q^2}\right)\left[r\right],a=6t+2$

输出：$R_a\left(Q,P\right)$

1) 设窗口长度r > 1

2) 预计算$f_1=1,T_1=Q;f_2=g_{Q,Q}\left(P\right),T_2=\left[2\right]Q$

3) For i from 1 to $2^{r-1}-1$ do:

4) $f_{2i+1}=f_{2i-1}{f}_2{g}_{\left[2i-1\right]Q,\left[2\right]Q}\left(P\right);T_{2i+1}=T_{2i-1}+T_2$

5) end for

6) 设

7) 置f = 1, T = Q, j = L − 2

8) While j ≥ 0 do:

9) if $a_i=0$

10) 计算$f=f^2\cdot g_{T,T}\left(P\right)$ ，$T=\left[2\right]T$

11) else

12) 令t为$j-t+1\ge r$ 且$a_t=1$ 的最小整数

13) $h_j={\displaystyle {\sum }_{i=0}^{j-t}{a}_{t+i}}{2}^i$

14) For i from j – t + 1 to 0 do

15) $f=f^2\cdot g_{T,T}\left(P\right)$ ，$T=\left[2\right]T$

16) end for

• 17) $f\leftarrow f\cdot f_{T_{h_j}}{g}_{T,T_{h_j}}\left(P\right);T\leftarrow T+T_{h_j}$

18) 置j = t − 1

19) end if

20) end while

21) 计算$Q_1={\pi }_q\left(Q\right),Q_2={\pi }_{q^2}\left(Q\right)$

22) 计算$f=f\cdot g_{T,Q_1}\left(P\right),T=T+Q_1$

23) 计算$f=f\cdot g_{T,-Q_2}\left(P\right),T=T-Q_2$

24) 计算$f=f^{\left(q^{12}-1\right)/r}$

25) 输出f

采用算法4在窗口长度为3的情况下可以将R-ate对运算过程中点加运算以及直线函数的运算数量从传统Miller算法中的16次降低至7次，相对比NAF算法进一步减少4次点加运算以及稀疏乘法运算，相对地只增加了2次模乘运算。具体算法所对应运算数量如表3所示。

Table 3. The Miller loop uses the number of operations

表3. Miller循环调用模块量

方法	点加–直线函数	f∙g	f∙f
传统方法	16	16	0
NAF方法	11	11	0
滑动窗口方法	7	7	2

Figure 1. Graph of the results of the R-ate pairing simulation

图1. R-ate对仿真模拟结果图

根据表1对比计算量分析可知，滑动窗口算法相比于传统算法节省549M + 1159A，相比于NAF算法节省184M + 304A。由此说明本文提出的基于滑动窗口的R-ate对计算方法在优化效果上要优于使用NAF改进的R-ate对计算方法。为了验证提出的实现方法的正确性，本文在Vivado 2018.3平台，使用verilog语言进行编译对实现方法进行仿真验证。在底层运算均相同的情况下，采用串行运算条件下，本文算法实现后仿真模拟结果如图1所示，进行一次双线性对运算所需时间为6.865 ms，相比于优化前运算时间8.42 ms计算效率提升约为18.46%。

4.3. 基于Montgomery的点乘算法快速实现

根据第2章提到的Montgomery模乘的特性，可以采用相同方法改进SM9算法中G₁群中的倍点运算以及点加运算。将椭圆曲线上点的运算转换为蒙哥马利形式进行运算，减少采用模乘运算过程中参数形式反复转换的步骤，进而提升椭圆曲线上倍点运算以及点加运算的计算效率。具体改进后倍点运算计算过程如算法5所示，点加运算计算过程如算法6所示。

算法5：G₁群倍点运算

输入：P₁ = (X₁, Y₁, Z₁)

输出：P₃ = (X₃, Y₃, Z₃)

1) 将坐标点X₁、Y₁分别转换为Montgomery形式X₁R、Y₁R

2) T₁ = (Y₁R)²

3) T₂ = T₁·X₁R

4) T₂ = 4·T₂

5) T₃ = (T₁)²

6) T₃ = 8·T₃

7) T₁ = X₁R·X₁R

8) T₁ = 3·T₁

9) X₃ = (T₁)²

10) X₃ = X₃ − T₂

11) X₃ = X₃ − T₂

12) Y₃ = T₂ − X₃

13) Y₃ = Y₃·T₁

14) Y₃ = Y₃ − T₃，将X₃、Y₃转换为素数域形式

15) Z₃ = Y₁R·Z₁

16) Z₃ = Z₃ + Z₃

算法5中，在不改变变量之间映射关系的条件下，将原有的有限域模乘算法拆分为更小的Montgomery模乘基础运算，通过将变量在素数域以及Montgomery域之间来回转换以减少原来算法中的冗余计算，进而提高倍点运算的计算效率。对比分析可知，通过算法5可以将原有倍点运算的计算量由7M + 12A减少至5.5M + 12A。

算法6：G₁群点加运算

输入：P₁ = (X₁, Y₁, Z₁), P₂ = (X₂, Y₂)

输出：P₃ = (X₃, Y₃, Z₃)

1) 将坐标点Z₁转换为Montgomery形式Z₁R

2) T₁ = (Z₁R)²

3) T₂ = T₁·Z₁R

4) T₁ = T₁·X₂

5) T₂ = T₂·Y₂

6) T₁ = T₁ − X₁

7) T₂ = T₂ − Y₁

8) Z₃ = T₁·Z₁R

9) 将T₁转换为Montgomery形式赋值给T₅

10) T₃ = (T₅)²

11) T₄ = T₃·T₅

12) T₃ = T₃·X₁

13) T₁ = T₃ + T₃

14) 将T₂转换为Montgomery形式赋值给T₂

15) X₃ = (T₂)²

16) X₃ = X₃ − T₄

17) 将X₃转换为Montgomery形式赋值给X₃

18) X₃ = X₃ − T₁

19) T₃ = T₃ − X₃

20) T₃ = T₃·T₂

21) T₄ = Y₁·T₄

22) Y₃ = T₃ − T₄

算法6采用算法5相同原理，通过两种形式的转换，可以将点加运算的计算量由11M + 7A下降至7.5M + 7A。在此基础上，通过NAF算法对进行点乘运算的参数进行改进可以进一步减少点乘过程中调用点加运算的次数，进而提升计算效率。

5. 实验结果

本文提出的快速实现通过Vivado 2018.3在Xilinx FPGA上基于verilog语言完成综合实现。通过第3章对比可知，采用文本提出的算法在串行条件下进行R-ate双线性对的计算可以将一次R-ate对运算时间由8.42 ms下降至6.865 ms，整体效率提升约为18.46%，该方法同样可以在其他串行运算条件上使用，可以将R-ate对Miller循环部分中进行点加运算、直线函数运算以及稀疏乘法的次数由16次下降到7次，提升效果相对于原有NAF改进方法更高。

在此基础上，提出根据Montgomery模乘特性改进有限域模逆运算、G₁群中点加运算以及倍点运算，同时采用NAF方法进一步减少点乘运算参数的汉明重量，以提升计算效率。通过在Vivado平台采用相同底层算法分别实现现有文献以及本文的快速实现方法，各模块运算所需时间如表4所示。

Table 4. Operation time of each module (µs)

表4. 各模块运算所需时间(µs)

	二次扩域模逆	G1群倍点运算	G1群点加运算	G1群点乘运算
改进前	4.2145	2.6625	4.0525	1156.535
改进后	2.4775	2.0675	2.7175	731.058

最后，本文根据《SM9标识密码算法参数定义》使用verilog语言实现SM9数字签名算法。实验结果表明，采用本文优化后的椭圆曲线运算方法，可以将G₁群倍点运算和点加运算效率分别提升22.35%以及32.94%，结合NAF算法优化点乘运算后，可以将SM9数字签名过程中的点乘运算效率提升约36.79%。在此基础上结合前文提出的基于滑动窗口改进的R-ate对计算方法，可以将串行条件下SM9数字签名算法的运算效率提升约13.55%。其中签名中各运算所需时间如表5所示，数字签名实验仿真结果如图2所示。

Table 5. Each operation time in SM9 digital signature (ms)

表5. SM9数字签名中各运算所需时间(ms)

	R-ate对	模幂运算	点乘运算	SM9数字签名
改进前	8.421	5.484	1.156	15.144
改进后	6.865	5.484	0.732	13.092

Figure 2. SM9 digital signature simulation result

图2. SM9数字签名仿真模拟结果图

6. 结论

R-ate对的计算效率对SM9算法的广泛普及以及配对在密码学中的应用具有重要意义，本文通过对R-ate对的运算过程进行分析，提出一种基于滑动窗口的R-ate对快速实现方法，其中包括基于滑动窗口算法优化R-ate对计算中Miller循环部分以及使用Montgomery特性改进底层椭圆曲线上点的计算和扩域运算。通过将原有运算复杂度与改进后算法进行分析，证明本文采用实现方法对R-ate对运算效率的提升有效。实验结果证明，本文所提出的实现方法可以将R-ate对以及SM9数字签名算法的计算效率分别提升约18.46%以及13.55%。

基金项目

国家自然科学基金(62472040)；中国版权保护中心版权研究课题(BQ2024017)；北京市教委科研计划资助(No. KM202110015004)；北京市高等教育学会2022年立项面上课题(MS2022093)；北京市教育委员会科学研究计划项目资助(KM202310015002)。

NOTES

^*通讯作者。