CN112367172A - 一种基于Co-Z运算的蒙哥马利梯形算法的AVX2快速实现方法 - Google Patents

Abstract

本发明公开了一种基于Co‑Z运算的蒙哥马利梯形算法的AVX2快速实现方法：首先计算蒙哥马利梯形算法中的两个Z坐标相同的椭圆曲线点R₁＝2P和R₀＝P；然后预计算A’和T’以调整Co‑Z Ladder算法的计算步骤；然后将R₁、R₀、A’和T’的转换成字长为26位的冗余表示法并存储到AVX2的256位寄存器中；利用基于AVX2实现的双向素数域运算对Co‑Z共轭加法、Co‑Z加法运算以及Co‑Z Ladder算法中的对称运算进行双向实现，并将其应用到蒙哥马利梯形算法中；最后将随机点的标量乘法运算结果从字长为26位的冗余表示法转换成普通的多精度表示法并转化成仿射坐标。本发明利用AVX2的计算能力和并行能力达到在不增加存储空间的前提下实现轻量级、安全且高效的基于Co‑Z运算的蒙哥马利梯形算法的目的。

Description

一种基于Co-Z运算的蒙哥马利梯形算法的AVX2快速实现方法

技术领域

本发明涉及一种在椭圆曲线公钥密码算法体制中的标量乘法运算的快速、轻量级实现方 法，特别涉及一种基于Co-Z运算的蒙哥马利梯形算法的AVX2快速实现方法，属于计算机领 域中的网络空间安全领域。

背景技术

蒙哥马利梯形算法最初是为特定的蒙哥马利曲线提出的一种固定时间的标量乘法运算， 在侧信道攻击发展较快的背景下，蒙哥马利梯形算法由于具有轻量级和固定时间的特点，因 此不需要过多的预计算空间且可以抵抗侧信道攻击，因此该算法也被广泛地推广应用到普通 的Weierstrass曲线中。

Co-Z运算是椭圆曲线群运算中的一部分，主要包括Co-Z初始化倍点、Co-Z加法、Co-Z 共轭加法、Co-Z取逆运算。它是一种将椭圆曲线点的射影坐标(X,Y,Z)中的Z坐标转换为相等 值从而提高椭圆曲线群运算的运行效率的一种技术。其中Co-Z加法(R,P′)＝P+Q仅需要 5MUL+2SQR+7ADD(MUL：串行乘法；SQR：串行平方；ADD：串行加法/减法；后文中还 会出现这三种运算的并行实现符号，分别为：MUL’：双向并行乘法；SQR’：双向并行平方；ADD’：双向并行加法/减法)即可完成计算并且更新P为等价点P′使其Z坐标与R相等。而Co-Z共 轭加法(R,R′)＝(P+Q,P-Q)除了计算两点P、Q的和P+Q以外还计算了两点之差，并且维持P+Q和P-Q的Z坐标相等。Venelli和Dassance则进一步地提出在蒙哥马利梯形算法的主循环中，Co-Z共轭加法和Co-Z加法中的Z坐标由于不参与结果的X、Y坐标计算，因此提出只涉及 X、Y坐标的Co-Z运算，减少了这两个Co-Z运算中计算Z坐标的乘法运算，Co-Z加法和 Co-Z共轭加法输出的两个点的Z坐标均保持相等这个特性使得Co-Z运算可以很好地被应用 到标量乘法运算中，特别是蒙哥马利梯形算法。因此，相对于原版蒙哥马利梯形算法，利用 Co-Z运算可以降低一定数量的运算次数并提升蒙哥马利梯形算法的运行效率。

AVX2是Intel处理器从Haswell架构开始支持的一种可处理256位数据的单指令多数据型指 令集。AVX2提供了多种可处理256位整型向量数据的指令，例如加法、乘法、移位、按变量 移位、排列、组合等，这些指令可同时处理256位数据，即4个64位向量数据或者8个32位向量 数据，AVX2的并行性为椭圆曲线密码算法的快速实现提供了极大的潜能。

发明内容

发明目的：利用AVX2的并行性，设计一种基于Co-Z运算的蒙哥马利梯形算法的AVX2 快速实现，使得这种轻量级且抗侧信道攻击的标量乘法运算能够在Intel处理器上提高运行效 率。

技术方案：一种基于Co-Z运算的蒙哥马利梯形算法的AVX2快速实现方法充分利用AVX2的并行性能将Co-Z运算进行并行化，降低Co-Z运算的计算时间，进而将利用AVX2 进行加速后的Co-Z运算运用到基于Co-Z运算的蒙哥马利梯形算法中，以提升其运行效率。 该实现包含以下计算步骤：

一种基于Co-Z运算的蒙哥马利梯形算法的AVX2快速实现方法，具体步骤如下：

步骤1，由于蒙哥马利梯形算法需要标量k的最高位的一个比特为k₂₅₅＝1，因此可以利用 Co-Z初始化倍点算法XYCZ_IDBL计算蒙哥马利梯形算法中的两个Z坐标相同的椭圆曲线点 R₁＝2P＝(X₁,Y₁,Z)和R₀＝P＝(X₀,Y₀,Z)，其中P是蒙哥马利梯形算法中输入的随机椭圆曲线 点，R₀是通过XYCZ_IDBL算法更新Z坐标后与P等价的椭圆曲线点，X₀、Y₀分别表示R₀的横、纵坐标，R₁则是通过XYCZ_IDBL算法得到的与2P等价的椭圆曲线点，X₁、Y₁分别 表示R₁的横、纵坐标；

步骤2，通过椭圆曲线点R₁和R₀的横坐标X₁、X₀预计算大整数A′＝(X₀-X₁)²和 T′＝(X₀-X₁)A′；

步骤3，将椭圆曲线点R₁、R₀以及大整数A'、T'转换成字长为26位的冗余表示法，并存 储到AVX2的256位寄存器中；

步骤4，进入蒙哥马利梯形算法的主循环，主循环的索引为i∈[1,254]，i从标量k的第 254个比特开始往下遍历k的每个比特，每次循环获取k的第i个比特k_i和第(i+1)个比特k_i+1， 并令符号a＝(k_i+k_i+1)mod 2，a∈{0,1}，通过a调整R₁、R₀的次序R_a、R_1-a，与A'、T'同时作为双向并行Co-Z Ladder算法的输入，并将双向并行Co-Z Ladder算法的输出分别更新到 R_a、R_1-a、A'、T'中，循环调用双向并行Co-Z Ladder算法并不断更新R_a、R_1-a、A'、T'直到 循环结束；

步骤5，主循环结束后，获取k最低位的两个比特k₀、k₁，并令符号a＝(k₀+k₁)mod2，，通 过a调整R₁、R₀的次序R_a、R_1-a，与A'、T'同时作为双向并行Co-Z共轭加法的输入，调用 双向并行Co-Z共轭加法，并将双向并行Co-Z共轭加法的输出分别更新到R_a、R_1-a中；

步骤6，将R_1-a,R_a,P,a作为串行Co-Z取逆运算的输入进行取逆运算得到结果

其中 Z＝x_PY_a(X₀-X₁),λ＝y_PX_a，x_P、y_P表示P在仿射坐标下的横、纵坐标，X_a、Y_a表示R_a在 射影坐标下的横、纵坐标，X₀、X₁则分别是R₀、R₁在射影坐标下的横坐标，λ是一个大整数；

步骤7，将R₀、R₁作为双向并行Co-Z加法的输入，调用双向并行Co-Z加法，并将双向并行Co-Z加法的结果更新到R₀、R₁中，最后通过坐标转换

求得蒙哥马利梯形算法仿射坐标形式的结果；

步骤8，将蒙哥马利梯形算法的结果从字长为26位的冗余表示法转换成字长为32位的多 精度表示法。

进一步，双向并行Co-Z加法的实现步骤如下(以下符号除了椭圆曲线点P、Q、R外，均为大整数，X_P,Y_P,X_Q,Y_Q,X_R,Y_R表示椭圆曲线点P、Q、R的横纵坐标，其它符号表示的是在 计算过程中的中间值，在计算过程中可能会多次更新)：

(1)Co-Z加法运算的输入分别为椭圆曲线点P＝(X_P,Y_P)、椭圆曲线点Q＝(X_Q,Y_Q)；

(2)利用双向并行减法运算，同时计算T₁＝X_P-X_Q和T₂＝Y_Q-Y_P；

(3)对T₁、T₂利用双向并行平方运算计算其平方，分别得到A＝T₁ ²和D＝T₂ ²；

(4)利用双向并行乘法运算，同时计算B＝X_PA和C＝X_QA；

(5)利用双向并行加法运算，计算T₁＝B+C；

(6)利用双向并行减法运算，同时计算X_R＝D-T₁和T₃＝C-B；

(7)利用双向并行减法运算，计算T₁＝B-X_R；

(8)利用双向并行乘法运算，计算T₁＝T₁T₂和计算E＝Y_PT₃；

(9)利用双向并行减法运算，计算Y_R＝T₁-E；

(10)得到椭圆曲线点R＝P+Q＝(X_R,Y_R)，并更新P的坐标为(B,E)。

进一步，双向并行Co-Z共轭加法的实现步骤如下(以下符号除了椭圆曲线点P、Q、R、 R’外，均为大整数，X_P,Y_P,X_Q,Y_Q,X_R,Y_R,X_R′,Y_R′表示椭圆曲线点P、Q、R、R’的横纵坐标， 其它符号表示的是在计算过程中的中间值，在计算过程中可能会多次更新)：

(1)Co-Z共轭加法的输入分别为椭圆曲线点P＝(X_P,Y_P)、椭圆曲线点Q＝(X_Q,Y_Q)，预 计算的大整数A′＝(X_Q-X_P)²,T′＝(X_Q-X_P)A′；

(2)利用双向并行乘法运算，同时计算C＝X_QA′和E＝Y_PT′；

(3)利用双向并行减法运算，同时计算B＝C-T′和T₁＝Y_Q-Y_P；

(4)利用双向并行加法运算，同时计算T₂＝B+C和T₃＝Y_P+Y_Q；

(5)利用双向并行平方运算，同时计算

和

(6)利用双向并行减法运算，同时计算X_R＝D-T₂和X_R′＝F-T₂；

(7)利用双向并行乘法运算，同时计算T₃＝T₁T₂和T₃＝T₃T₄；

(8)利用双向并行减法运算，同时计算Y_R＝T₂-E和Y_R′＝T₃-E；

(9)，得到椭圆曲线点R＝P+Q＝(X_R,Y_R)和R′＝P-Q＝(X_R′,Y_R′)。

在执行完Co-Z共轭加法的计算流程后，不需要再更新A’和T’，因为这两项仅参与计 算Co-Z共轭加法，而在基于Co-Z运算的蒙哥马利梯形算法(如图4)中，Co-Z共轭加法后只需要再执行一次Co-Z加法，该算法不需要使用预计算项A′、T′。

进一步，双向并行Co-Z Ladder算法的实现步骤如下(以下符号除了椭圆曲线点P、Q、R、 R’外，均为大整数，X_P,Y_P,X_Q,Y_Q,X_R,Y_R,X_R′,Y_R′表示椭圆曲线点P、Q、R、R’的横纵坐标，其 它符号表示的是在计算过程中的中间值，在计算过程中可能会多次更新(包括P点的坐标值 也可能会被更新)。另外由于该算法是Co-Z共轭加法和Co-Z加法的合并算法，注意到这两 个算法在运行的过程中有多个重复的项，如B、C、D、E、F等，为了区分这些重复的项，该算法中B′，C′，D′，E′，F′表示在计算Co-Z共轭加法过程中产生的项，B，C，D，E，F 则表示在计算Co-Z加法过程中产生的项)：

(1)Co-Z Ladder算法的输入分别为椭圆曲线点P＝(X_P,Y_P)、椭圆曲线点Q＝(X_Q,Y_Q)， 预计算的大整数A′＝(X_Q-X_P)²,T′＝(X_Q-X_P)A′；

(2)利用双向并行乘法运算，同时计算C′＝X_QA′和E′＝Y_PT′；

(3)利用双向并行减法运算，同时计算B′＝C′-T′和T₁＝Y_P-Y_Q；

(4)利用双向并行加法运算，同时计算T₂＝B′+C′和T₃＝Y_P+Y_Q；

(5)利用双向并行平方运算，同时计算

和

(6)利用双向并行减法运算，同时计算X_R＝D′-T₂和X_R′＝F′-T₂；

(7)利用双向并行减法运算，同时计算T₂＝B′-X_R和T₄＝X_R′-B′；

(8)利用双向并行乘法运算，同时计算T₂＝T₁T₂和T₄＝T₃T₄；

(9)利用双向并行减法运算，同时计算Y_R＝T₂-E′和Y_R′＝T₄-E′；

(10)利用双向并行减法运算，同时计算T₁＝X_R-X′_R和T₂＝Y_R-Y_R′；

(11)利用双向并行平方运算，同时计算

和

(12)利用双向并行乘法运算，同时计算X_P＝B＝X_RA和C＝X_R′A；

(13)利用双向并行加法运算，同时计算Y₃＝T₂+B和T₄＝B+C；

(14)利用双向并行减法运算，同时计算X_Q＝D-T₄和T₁＝C-B；

(15)利用双向并行减法运算，同时计算T₄＝X_Q-X_P和T₃＝T₃-X_Q；

(16)利用双向并行平方运算，同时计算

和

(17)利用双向并行乘法运算，同时计算T′＝T₄A′和E＝Y_RT₁；

(18)利用双向并行加法运算，同时计算T₁＝D+A′和T₂＝E+E；

(19)利用双向并行减法运算，计算T₃＝T₃-T₁；

(20)利用双向并行减法运算和右移位运算，计算

(21)得到椭圆曲线点2P＝(X_P,Y_P)和P+Q＝(X_Q,Y_Q)。注意此时的(X_P,Y_P)、(X_Q,Y_Q)和输 入的点P、Q的坐标虽然符号相同，但是值已经不同，在计算过程中会更新每个坐标的值。在计算完Co-Z Ladder算法的计算流程后，A′、T′还作为输出传递到下一轮循环中参与Co-ZLadder算法或者Co-Z共轭加法的计算。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

(1)在底层素数域运算中，利用AVX2实现的双向并行素数域运算均比单向素数域运 算快，而且加速比在1.2～1.6之间，即执行用AVX2实现的双向并行素数域运算的时钟周期比 执行两次相同的单向素数域运算要快1.2～1.6倍；

(2)在Co-Z运算和基于Co-Z运算的蒙哥马利梯形算法的实现上，底层素数域运算的 效率提升也使得并行Co-Z加法的效率比普通Co-Z加法提升1.26倍，并行Co-Z共轭加法的 效率则提升了1.60倍，而这Co-Z Ladder算法的效率则提升了1.33倍；而利用并行的Co-Z 加法、共轭加法和Co-Z Ladder算法实现的基于Co-Z运算的蒙哥马利梯形算法的运行效率也 提升了1.31倍；

(3)基于Co-Z运算的蒙哥马利梯形算法在提升效率的同时并没有增加存储空间的需 求，而且其执行时间也是固定时间，可以在保证轻量级和安全的前提下很好地提高算法的运 行效率，这都归功于AVX2指令集的合理运用。

附图说明

图1是并行Co-Z加法的AVX2实现算法；

图2是并行Co-Z共轭加法的AVX2实现算法；

图3是并行Co-Z Ladder的AVX2实现算法；

图4是基于Co-Z运算的蒙哥马利梯形算法的AVX2快速实现算法；

图5是双向并行乘法的AVX2实现算法；

图6是SM2素数域的快速约减算法的AVX2实现；

图7是快速约减算法中的SimpleAdd算法的AVX2实现；

图8是快速约减算法中根据SM2素数p建立起来的等价关系；

图9是SM2素数域的系数约减算法的AVX2实现；

图10是本发明的方法流程图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

为了利用AVX2实现轻量级、安全且高效的基于Co-Z运算的蒙哥马利梯形算法，本发明 提出的方案包含三个部分：1.基于SM2素数域使用AVX2实现的双向并行素数域运算；2.重 新调度Co-Z加法、Co-Z共轭加法以及Co-Z Ladder算法，使用双向并行的素数域运算同时 处理两组对称的运算；3.利用并行的Co-Z加法、Co-Z共轭加法以及Co-Z Ladder算法实现 基于Co-Z运算的蒙哥马利梯形算法。

下面先介绍第2、3部分，再介绍底层素数域运算的实现。

一、只涉及(X,Y)坐标的双向并行Co-Z运算

根据Meloni提出的Co-Z运算，令P＝(X_P,Y_P,Z),Q＝(X_Q,Y_Q,Z),其中X_P,Y_P,Z分别是椭 圆曲线点P的横、纵、Z坐标，X_Q,Y_Q,Z分别是椭圆曲线点Q的横、纵、Z坐标，则Co-Z加 法R＝P+Q可由以下公式进行计算：

X_R＝D-(B+C),Y_R＝(Y_Q-Y_P)(B-X_R)-E,Z_R＝Z(X_Q-X_P)

A＝(X_Q-X_P)²,B＝X_PA,C＝X_QA,D＝(Y_Q-Y_P)²,E＝Y_P(C-B)

公式中的每一项字母A、B、C、D、E、F都是一个256位的大整数，X_R,Y_R,Z_R分别是Co-Z加法得到的椭圆曲线点R的横、纵、Z坐标。Co-Z加法不仅计算了R＝P+Q＝(X_R,Y_R,Z_R)，还 更新P的坐标为(X_P′,Y_P′,Z_R)＝(B,E,Z_R)使得P的Z坐标与R相等。

而共轭加法(R,R′)＝((X_R,Y_R,Z_R),(X_R′,Y_R′,Z_R))＝(P+Q,P-Q)，除了需要进行以上运 算之外，还需要额外进行以下运算以计算R’＝P-Q，并且还需将两点的Z坐标转化为相等的。

X_R′＝F-(B+C),Y_R′＝(Y_P+Y_Q)(X_R′-B)-E

F＝(Y_P+Y_Q)²,A,B,C,D,E与以上公式相等。

观察到以上计算公式中存在多组数据(如B和C、A和D等)的计算具有对称的特征，本发 明在仔细地研究计算公式的数据相关性之后，将以上两个公式的计算顺序进行重新调度，并 得出图1～图3的Co-Z加法、Co-Z共轭加法、Co-Z Ladder算法的并行实现，注意Z坐标的 计算由于可以在蒙哥马利梯形算法最后恢复，因此可以减少计算Z坐标的步骤，图1～图3中 每行的运算都是利用AVX2指令集同时计算的，具体实现方法将在素数域运算部分详细讲解。

Co-Z加法并行实现(SIMD_XYCZ_ADD)步骤如图1所示：

以下符号除了椭圆曲线点P、Q、R外，均为大整数，其它符号表示的是在计算过程中的 中间值，在计算过程中可能会多次更新。

(1)Co-Z加法运算的输入分别为椭圆曲线点P＝(X_P,Y_P)，椭圆曲线点Q＝(X_Q,Y_Q)；

(2)利用双向并行减法运算同时计算T₁＝X_P-X_Q和T₂＝Y_Q-Y_P；

(3)再对T₁、T₂利用双向并行平方运算计算其平方，分别得到A＝T₁ ²和D＝T₂ ²；

(4)利用双向并行乘法运算同时计算B＝X_PA和C＝X_QA；

(5)利用双向并行加法运算计算T₁＝B+C；

(6)利用双向并行减法运算同时计算X_R＝D-T₁和T₃＝C-B；

(7)利用双向并行减法运算计算T₁＝B-X_R；

(8)利用双向并行乘法运算同时计算T₁＝T₁T₂和E＝Y_PT₃；

(9)最后，利用双向并行减法运算计算Y_R＝T₁-E；

至此，Co-Z加法中所有项均已求得，得到椭圆曲线点R＝P+Q＝(X_R,Y_R)并更新P的坐标为 (B,E)。

考虑到蒙哥马利梯形算法的主循环中需要先执行一次Co-Z共轭加法和一次Co-Z加法运 算，因此本发明还将这两种加法的运算进行合并实现(即Co-Z Ladder算法)，以减少在两种 算法中利用AVX2进行优化时出现冗余运算的现象(即部分运算只能对一组数据进行计算， 这样会浪费部分AVX2的计算能力。如Co-Z加法中(图1)第5、7步只能对一组数据进行 计算)，提高运算效率。Co-Z Ladder算法的并行实现如图3所示，其中大整数B’～F’分别 为Co-Z Ladder算法中计算Co-Z共轭加法的中间结果，大整数A～E则为Co-Z Ladder算法中 计算Co-Z加法的中间结果。由于Co-Z共轭加法在蒙哥马利梯形算法中是在Co-ZLadder算 法后进行的，而且Co-Z共轭加法需要使用到Co-Z Ladder算法的中间结果A’和T’，因此 首先介绍Co-Z Ladder算法的过程。

Co-Z Ladder算法并行实现(SIMD_XYCZ_ADDC_ADD)步骤，如图3所示：

以下符号除了椭圆曲线点P、Q、R外，均为大整数，其它符号表示的是在计算过程中的 中间值，在计算过程中可能会多次更新。另外由于该算法是Co-Z共轭加法和Co-Z加法的合 并算法，注意到这两个算法在运行的过程中有多个重复的项，如B、C、D、E、F等，为了 区分这些重复的项，该算法中B′、C′、D′、E′、F′表示在计算Co-Z共轭加法过程中产生的项， B，C，D，E，F则表示在计算Co-Z加法过程中产生的项。

(1)首先观察到等式Y_Q＝(Y_R-Y_R′)(B-X_Q)-E中可以通过2Y_Q＝(Y_R′-Y_R+B- X_Q)²-(Y_R′-Y_R)²-A′-2E替换，其中(Y_R′-Y_R)²、A’都是在计算Y_Q前必须计 算的项，通过这种方法可以重复利用这两项来减少一次乘法，并将该乘法替换成1 个平方和5个加法，此时Co-Z Ladder算法的计算量降低到8MUL+6SQR+23ADD， 由于所需乘法次数和平方次数均为偶数，因此利用双向素数域运算最好可以将计 算复杂度降低一半；

(2)该算法输入分别为椭圆曲线点P＝(X_P,Y_P)，椭圆曲线点Q＝(X_Q,Y_Q)，预计算的大整 数A′＝(X_Q-X_P)²,T′＝(X_Q-X_P)A′,大整数A’和T’的预计算可以使得调整运 算次序后的Co-Z Ladder算法无过多的冗余计算，计算量可以降低一半到 4MUL’+3SQR’。而直接地并行实现Co-Z共轭加法和Co-Z加法并不能将算法 计算量降低到最低，最低只能降低到5MUL’+2SQR’。由于Co-Z Ladder算法在 主循环中需要循环254次，因此通过预计算两项的优化方案来并行化Co-Z Ladder 算法能带来极大的效率提升；

(3)利用双向并行乘法运算同时计算C′＝X_QA′和E′＝Y_PT′；

(4)利用双向并行减法运算同时计算B′＝C′-T′和T₁＝Y_P-Y_Q；

(5)利用双向并行加法运算同时计算T₂＝B′+C′和T₃＝Y_P+Y_Q；

(6)利用双向并行平方运算同时计算

和

(7)利用双向并行减法运算同时计算X_R＝D′-T₂和X_R′＝F′-T₂；

(8)利用双向并行减法运算同时计算T₂＝B′-X_R和T₄＝X_R′-B′；

(9)利用双向并行乘法运算同时计算T₂＝T₁T₂和T₄＝T₃T₄；

(10)利用双向并行减法运算同时计算Y_R＝T₂-E′和Y_R′＝T₄-E′；

(11)利用双向并行减法运算同时计算T₁＝X_R-X_R′和T₂＝Y_R-Y_R′；

(12)利用双向并行平方运算同时计算

和

(13)利用双向并行乘法运算同时计算X_P＝B＝X_RA和C＝X_R′A；

(14)利用双向并行加法运算同时计算T₃＝T₂+B和T₄＝B+C；

(15)利用双向并行减法运算同时计算X_Q＝D-T₄和T₁＝C-B；

(16)利用双向并行减法运算同时计算T₄＝X_Q-X_P和T₃＝T₃-X_Q；

(17)利用双向并行平方运算同时计算

和

(18)利用双向并行乘法运算同时计算T′＝T₄A′和E＝Y_RT₁；

(19)利用双向并行加法运算同时计算T₁＝D+A′和T₂＝E+E；

(20)利用双向并行减法运算计算T₃＝T₃-T₁；

(21)利用双向并行减法运算计算

至此，Co-Z Ladder算法中的所有项均已求得，得到椭圆曲线点2P＝(X_P,Y_P)和 P+Q＝(X_Q,Y_Q)，注意此时的(X_P,Y_P)、(X_Q,Y_Q)和输入的点P、Q的坐标虽然符号相同，但是值 已经不同，在计算过程中会更新每个坐标的值。在计算完Co-Z Ladder算法的计算流程后， A′、T′还作为输出传递到下一轮循环中参与Co-Z Ladder算法或者Co-Z共轭加法的计算。

Co-Z共轭加法并行实现(SIMD_XYCZ_ADDC)步骤，如图2所示：

以下符号除了椭圆曲线点P、Q、R、R′，均为大整数，其它符号表示的是在计算过程中 的中间值，在计算过程中可能会多次更新。

(1)如图4所示，Co-Z共轭加法在主循环结束后才需要执行，而Co-Z Ladder算法在主循环的最后一轮循环还更新了A′和T′，这两项也可以参与到Co-Z共轭加法中以 调整计算次序，由于A′和T′的参与，Co-Z共轭加法中所有的计算均可以充分利用 AVX2的并行性两两进行同时运算；

(2)该算法输入分别为椭圆曲线点P＝(X_P,Y_P)，椭圆曲线点Q＝(X_Q,Y_Q)，预计算的大整 数A′＝(X_Q-X_P)²,T′＝(X_Q-X_P)A′；

(3)利用双向并行乘法运算同时计算C＝X_QA′和E＝Y_PT′；

(4)利用双向并行减法运算同时计算B＝C-T′和T₁＝Y_Q-Y_P；

(5)利用双向并行加法运算同时计算T₂＝B+C和T₃＝Y_P+Y_Q；

(6)利用双向并行平方运算同时计算

和

(7)利用双向并行减法运算同时计算X_R＝D-T₂和X_R′＝F-T₂；

(8)利用双向并行乘法运算同时计算T₃＝T₁T₂和T₃＝T₃T₄；

(9)利用双向并行减法运算同时计算Y_R＝T₂-E和Y_R′＝T₃-E；

至此，Co-Z加法中的所有项均已求得，得到椭圆曲线点R＝P+Q＝(X_R,Y_R)和R′＝P-Q＝ (Y_R′,Y_R′)。在执行完Co-Z共轭加法的计算流程后，不需要再更新A′和T′，因为这两项仅参与 计算Co-Z共轭加法，而在基于Co-Z运算的蒙哥马利梯形算法(图4)中，Co-Z共轭加法后 只需要再执行一次Co-Z加法，该算法不需要使用预计算项A′和T′。

实现效率：

关于单向和双向并行的素数域运算的时钟周期比较如表1所示，双向并行的素数域运算 的效率比单向素数域运算要提高1.2～1.6倍之间，这也为Co-Z加法、Co-Z Ladder算法的效率 提升做出了贡献，如表1所示，并行Co-Z加法(Co-Z ADD)、共轭加法(Co-ZADDC)以 及Co-Z Ladder算法(Co-Z ADDC_ADD)比顺序版本实现分别快了1.26倍、1.60和1.33倍。

表1 并行和顺序的Co-Z运算和基于Co-Z运算的蒙哥马利梯形算法的时钟周期(单位:clock cycles)比较

二、基于Co-Z运算的蒙哥马利梯形算法的AVX2实现

蒙哥马利梯形算法的功能是计算标量乘法kP，令椭圆曲线点P的坐标为P＝(X_P,Y_P,Z_P),椭 圆曲线点Q的坐标为Q＝(X_Q,Y_Q,Z_P)，为了方便后续公式的演示，我们将P、Q分别赋值给 R_a＝(X_a,Y_a,Z_a)和R_1-a＝(X_1-a,Y_1-a,Z_1-a)，其中a∈{0,1}，即R_a＝P，R_1-a＝Q，定义 R_a,R_1-a,a是为了在P和Q之间进行选择时方便表达。则蒙哥马利梯形算法的算法原理是在 主循环中维持以下等式：

(X_1-a,Y_1-a,Z_1-a)＝(X_1-a,Y_1-a,Z_1-a)+(X_a,Y_a,Z_a),(X_a,Y_a,Z_a)＝2(X_a,Y_a,Z_a) 在每个循环之后，两个点R_a,R_1-a的差值的绝对值均为P，而基于Co-Z运算的蒙哥马利梯形 算法则是用以下等式替换了原始等式：

((X_1-a,Y_1-a,Z),(X_a,Y_a,Z))＝((X_a,Y_a,Z)+(X_1-a,Y_1-a,Z),(X_a,Y_a,Z)-(X_1-a,Y_1-a,Z)),

((X_a,Y_a,Z),(X_1-a,Y_1-a,Z))＝(X_1-a,Y_1-a,Z)+(X_a,Y_a,Z)

基于Co-Z运算的蒙哥马利梯形算法在主循环中先做一个Co-Z共轭加法，然后执行一个Co-Z 加法。基于Co-Z运算的蒙哥马利梯形算法与原始蒙哥马利梯形算法的主要区别在于其椭圆曲 线层运算不同，前者基于Co-Z运算(如Co-Z加法、Co-Z共轭加法以及Co-ZLadder算法)， 而后者基于射影坐标上的椭圆曲线层运算(射影坐标上的椭圆曲线点加、倍点)。由于Co-Z 运算建立在射影坐标的Z坐标相同的情况下，部分涉及Z坐标的椭圆曲线运算因此可以减少 计算量，因此基于Co-Z运算的运行效率高，相应地基于Co-Z运算的蒙哥马利梯形算法的运 行效率也比原始的蒙哥马利梯形算法的好。因此本发明针对基于Co-Z运算的蒙哥马利梯形算 法进行优化。

如图4和10所示，基于Co-Z运算的蒙哥马利梯形算法的算法步骤为：

(1)该算法输入椭圆曲线点P和大整数标量k，且标量k的最高位的一个比特为k₂₅₅＝1；

(2)令(R₁,R₀)＝XYCZ_IDBL，其中XYCZ_IDBL为初始化倍点算法。为了使用Co-Z 运算，首先需要使用初始化倍点算法XYCZ_IDBL初始化椭圆曲线点R₁,R₀为2P,P， 并且两点的Z坐标相等，由于XYCZ_IDBL的运算对称性不好，因此该算法并未 使用AVX2进行优化，而只是实现串行版本的算法，该算法计算流程如下；

C＝A²,D＝3(B-1),E＝2((X_P+A)²-B-C),

X₁＝D²-2E,Y₁＝D(E-X_Q)-8C,Z₁＝2Y_P

此时R₁＝2P＝(X₁,Y₁,Z₁)，并更新P的坐标为并存储到R₀＝P＝(X₀,Y₀,Z₀)且

X₀＝E,Y₀＝8C,Z₀＝Z₁；

(3)然后，通过椭圆曲线点R₁和R₀的横坐标X₁,X₀预计算大整数A′＝(X₀-X₁)²和 T′＝(X₀-X₁)A′；

(4)再后，将椭圆曲线点R₁、R₀以及大整数A'、T'转换成字长为26位的冗余表示法，并存储到AVX2的256位寄存器中；

(5)再后，进入蒙哥马利梯形算法的主循环，主循环的索引为i∈[1,254]，i从k的第 254个比特开始往下遍历k的每个比特。每次循环获取大整数k的第i个比特k_i和 第(i+1)个比特k_i+1并令符号a＝(k_i+k_i+1)mod 2，通过a调整R₁,R₀的次序 R_a、R_1-a与A'、T'同时作为双向并行Co-Z Ladder算法的输入，并将双向并行Co-Z Ladder算法的输出分别更新到R_a、R_1-a、A'、T'中，反复调用双向并行Co-Z Ladder 算法不断更新R_a、R_1-a、A'、T'直到循环结束；

(6)再后，主循环结束后，获取大整数k的最低位的两个比特k₀、k₁并令符号a＝(k₀+k₁)mod2，通过a调整R₁,R₀的次序R_a、R_1-a与A'、T'同时作为双向并 行Co-Z共轭加法的输入调用双向并行Co-Z共轭加法，并将双向并行Co-Z共轭 加法的输出分别更新到R_a、R_1-a中；

(7)然后将R_1-a,R_a,P,a作为串行Co-Z取逆运算的输入进行取逆运算得到结果

其中 取逆运算的计算公式为Z＝x_PY_a(X₀-X₁),λ＝y_PX_a，x_P、y_P表示椭圆曲线点P仿 射坐标下的横纵坐标，X_a、Y_a表示椭圆曲线点R_a射影坐标下的横纵坐标，X₀、X₁则 分别是R₀,R₁射影坐标下的横坐标，λ是一个大整数；

(8)再将R₀,R₁作为双向并行Co-Z加法的输入调用双向并行Co-Z加法，并将双向并行 Co-Z加法的结果更新到R₀,R₁中，最后通过坐标转换

求 得蒙哥马利梯形算法仿射坐标形式的结果；

(9)最后，将蒙哥马利梯形算法的结果从字长为26位的冗余表示法转换成字长为32 位的多精度表示法。

实现效率：本发明中实现的基于Co-Z运算的蒙哥马利梯形算法(Co-ZMont.Ladder)， 其时钟周期对比如表1所示，其运行效率比原始版本快1.31倍。

三、底层双向并行素数域运算

该部分将介绍针对SM2素数域实现的底层双向并行素数域运算，素数域实现部分作为本 发明的基础运算部分，支撑着上层的Co-Z运算和蒙哥马利梯形算法的实现，是发明的重要组 成部分之一。SM2作为我国椭圆曲线公钥密码算法商业标准，未来几年在国内具有广泛的应 用前景，因此本发明根据SM2的椭圆曲线参数对底层素数域进行实现。需要注意的是，如果 能够实现相对应的底层双向并行素数域运算，本发明中高效安全的基于Co-Z运算的蒙哥马利 梯形算法也适用于其它Weierstrass椭圆曲线公钥密码算法标准，如NISTP-256。

1.AVX2中的大整数表示

SM2素数域运算主要包括加法、减法、乘法、平方、快速约减、取逆以及系数约减运算。 为了充分利用AVX2的并行性，我们选择使用位长为26位的冗余表示法来表示素数域中的元 素，这样可以保证每个系数的运算都具有独立、无相关性的特点，可以在AVX2指令集中并 行执行就算任务，使用26位字长保证了在每个系数上进行运算的结果不会超过64位。

SM2采用的素数域是256位的素数p＝2²⁵⁶-2²²⁴-2⁹⁶+2⁶⁴-1，由于采用字长为26 位的大整数表示方法，因此该素数域中的元素A由10个26位的系数组成，可以表示为

i为遍历素数域元素A每个系数的索引，a_i表示组成A的第i个系数。为了在AVX2中表示SM2素数域中的元素并且实现双向并行的素数域运算，我们使用交错元组<A,B> 来表示两个元素A和B，其中

a_i、b_i分别表示元素A、B的第i个 系数。交错元组共需要5个256位的AVX2寄存器存储，每个AVX2寄存器可表示为 <A,B>_i＝[a_2i+1,a_2i,b_2i+1,b_2i]，i为5个寄存器组的索引值，a_2i+1,a_2i,b_2i+1,b_2i分别表示元 素A、B的对应索引的系数。每个AVX2寄存器分别存储元素A、B的两个系数，每个系数 占用256位AVX2寄存器的64位，AVX2寄存器中4个系数的字长可以超过26位但不能超 过64位，否则单个256位的寄存器无法正确存储元素的4个系数，我们称交错元组<A,B>中 某个系数的字长超过26位的元组为未约减元组。

2.多精度乘法

素数域的运算中加法和减法比较简单，只要对两个元组中的5对寄存器分别进行加法或 减法运算即可。乘法算法和平方算法本发明采用的是两种传统的操作数扫描法和乘积扫描法， 其中平方算法相对于乘法用移位指令来替代了重复的乘法指令，因此平方算法相对较快一点。 多精度乘法运算如图5所示，具体实现步骤如下：

(1)双向并行的多精度乘法运算的输入为两组元组<A,B>、<C,D>，其中A、B、C、D 是SM2素数域中的256位的元素，双向并行的多精度乘法运算同时计算两组乘法 并输出<E,F>＝<A×C mod p,B×D mod p>；

(2)首先，利用AVX2的ALIGN指令对元组<C,D>进行重组，图5中第2-4行将元组 <C,D>的每个寄存器<C,D>_i＝[c_2i+1,c_2i,d_2i+1,d_2i]重组成[c_2i+2,c_2i+1,d_2i+2,d_2i+1] 并存储到寄存器<C′,D′>_i，5个寄存器<C′,D′>_i构成一个新元组<C’,D’>；

(3)然后，利用AVX2的SHUF指令对元组<A,B>进行重组，图5中第6行将元组<A,B>的每个寄存器<A,B>_i＝[a_2i+1,a_2i,b_2i+1,b_2i]重组成[a_2i,a_2i,b_2i,b_2i]并存储到寄存 器U中；图5中的第7-9行，利用AVX2的MUL和ADD指令将寄存器U和元组 <C,D>中的每个寄存器<C,D>_j进行相乘、累加运算，将结果存储到寄存器Z_i+j中；

(4)再后，利用AVX2的SHUF指令对元组<A,B>进行重组，图5中第10行将元组<A,B>的每个寄存器<A,B>_i＝[a_2i+1,a_2i,b_2i+1,b_2i]重组成[a_2i+1,a_2i+1,b_2i+1,b_2i+1]并存储 到寄存器V中；图5中的第11-13行，利用AVX2的MUL和ADD指令将寄存器 V和元组<C’,D’>中的每个寄存器<C′,D′>_j进行相乘、累加运算，将结果存储 到寄存器Z_i+j+1中；

(5)再后，图5中的第14-16行，利用AVX2的MUL指令将寄存器V和元组<C’,D’>的 第5个寄存器<C′,D′>₄进行乘法运算，并将结果通过BLEND和ADD指令调整 后分别累加到寄存器Z_i,Z_i+5中；遍历完元组<A,B>的5个寄存器<A,B>_i后，得到 需要10个AVX2寄存器存储的中间乘积(用Z表示)，为了将该中间乘积Z转换 成SM2素数域中的元素，需要对Z进行快速约减算法(fast_reduction)获得最乘 积结果<E,F>＝<A×C mod p,B×D mod p>。

3.快速约减算法

多精度乘法计算结束后，得到一个需要10个256位寄存器存储的中间乘积Z，为了将该 中间乘积Z转换成SM2素数域中的元素，本发明还针对SM2的素数p设计了在AVX2上快速实现的快速约减算法来将乘法和平方的中间乘积约减成SM2素数域中的元素。中间乘积Z的每个系数由于是由两个26位系数相乘得到的乘积累加的结果，因此每个系数的字长都大于 26位，属于未约减元组，未约减元组包含两个大整数，我们称这两个大整数为未约减元素。 对于未约减元素A中的第i个系数a_i，我们可以将a_i分成三个部分a_i＝h_i||m_i||l_i(由于a_i占据 AVX2寄存器的64位，因此我们可以将该64位的系数划分成3个部分，其中l_i表示a_i的低位 部分，m_i表示a_i的中间部分，h_i表示a_i的高位部分，||表示连接，即a_i由这三个部分连接而成)， 其中每个部分的比特长度为：|m_i|＝|l_i|＝26,|h_i|≤12。

快速约减算法如图6、7所示，具体实现步骤如下：

(1)快速约减算法的输入为一个需要10个AVX2寄存器存储的中间乘积Z，且其中的每个寄存器存储的系数的字长都超过52位(两个26位的系数的乘积字长为52位)， 因为在快速约减的过程中，需要将中间乘积Z的高5个寄存器进行一系列的移位 和累加到低5个寄存器中，因此需要先将中间乘积Z的高5个寄存器进行系数约 减，将高5个寄存器的系数字长约减到28位，这样可以保证在快速约减算法的后 续操作中不会出现溢出(即某个系数的字长超过64位)；

(2)首先，图6的第1-8行，利用AVX2的AND、SHR指令将中间乘积Z的高5个寄 存器Z₄-Z₉划分成三个部分L_i，M_i，H_i，其中L_i表示AVX2寄存器Z_i中的4个系数 的低位部分，M_i表示AVX2寄存器Z_i的4个系数的中间部分，H_i表示AVX2寄存 器Z_i的4个系数的高位部分。这三个符号L_i，M_i，H_i与上述的l_i，m_i，h_i略有不同， 后者表示的是单个系数的划分，是一个26位的系数；而前者则表示的是一个AVX2 寄存器中4个系数的划分，前者是一个AVX2寄存器；

(3)然后，图6的第10行，利用AVX2的ALIGN函数，将寄存器M_i和M_i-1进行重组 成

记作M′_i，其中

分别表示系数a_2i,a_2i-1的中间 部分

分别表示系数b_2i,b_2i-1的中间部分，并将M′_i、L_i和H_i-1进行累加得 到高5个新Z_i，经过(2)、(3)步骤，新Z_i里的未约减系数被大大约减，这可以使 得后续约减算法(SimpleADD)可以安全执行，而不会导致溢出(即某个系数的 字长超过64位)；

(4)再后，由于每个256位的大整数在本文中是由10个系数构成的，因此两组大整数 A、C和B、D的中间乘积都是由20个系数构成的，在多精度乘法运算中，我们 将A×C、B×D的中间乘积存储在10组寄存器Z中，令寄存器组中的每个寄存器

其中

表示A×C中间乘积的第2i+1、2i个系数，

表示B×D中间乘积的第2i+1、2i个系数。针对中间乘积Z的高5个寄 存器，里面包含两组系数

分别为A×C和B×D中间乘积 的高10个系数。为方便表示，将这两组系数的上标ac和bd省略，用z₁₀-z₁₉表 示其中一组大整数乘积的高10个系数。通过分析z₁₀-z₁₉的权重，利用SM2的模 数p建立起图8所示的等价关系，如z₁₀的权重是2²⁶⁰，由于p＝2²⁵⁶-2²²⁴-2⁹⁶+ 2⁶⁴-1，则2²⁵⁶≡2²²⁴+2⁹⁶-2⁶⁴+1modp，那么2²⁶⁰≡2²²⁸+2¹⁰⁰-2⁶⁸+ 2⁴modp，所以z₁₀2²⁶⁰≡z₁₀(2²²⁸+2¹⁰⁰-2⁶⁸+2⁴)modp，所以只需要在对应 权重(2⁴、2⁶⁸、2¹⁰⁰、2²²⁸)的系数中加/减系数z₁₀即可，这部分的算法如图7 (SimpleADD算法)所示。图7中的第1-7行，将图8中所有权重在2⁰～2⁵²的系 数经过移位或者重组累加到<E,F>₀中；图7中的第8-15行，将图8中所有权重 在2⁵²～2¹⁰⁴的系数经过移位或者重组累加到<E,F>₁中；图7中的第16-21行， 将图8中所有权重在2¹⁰⁴～2¹⁵⁶的系数经过移位或者重组累加到<E,F>₂中；图7 中的第22-25行，将图8中所有权重在2¹⁵⁶～2²⁰⁸的系数经过移位或者重组累加到 <E,F>₃中；图7中的第26-31行，将图8中所有权重在2²⁰⁸～2²⁶⁰的系数经过 移位或者重组累加到<E,F>₄中；

(5)经过SimpleADD算法得出的结果将10个AVX2寄存器约减到只需要5个寄存器存储的乘积<E,F>＝<A×C、B×D>，其中每个寄存器存储的系数在累加运算中字 长都超过了26位，因此<E,F>依然是未约减元组，因此如果后续还要进行乘法、 平方运算时，还需要对<E,F>进行进一步的约减——系数约减算法。

4.未约减元组的系数约减算法

在元组上进行的每个运算都需要保证元组中在每个系数上进行的运算结果不会超过64位， 否则单个256位的寄存器无法正确存储元素的4个系数。对于未约减元组，如果AVX2寄存 器中各个64位部分的系数字长超过在各个运算中定义的阈值，则需要进行系数约减。系数约 减算法如图9所示，该算法将未约减元组的各个系数的字长约减为28位，这个字长保证了该 元组在本发明中的所有SM2素数域中的运算都不会产生溢出。该算法如图9所示，其算法步 骤如下：

(1)系数约减算法的输入为一组未约减元组<A,B>，该元组的系数字长过大导致其无法 用于乘法、平方运算；

(2)首先，图9中的第1-4行，我们首先把其中权重超过260的系数a₈,a₉或者b₈,b₉进行约减，其中a₈超过52位的部分和a₉超过26位的部分的权重都是260，获取权重 超过260的部分存储到Q后，根据等价条件Q2²⁶⁰≡Q(2²²⁸+2¹⁰⁰-2⁶⁸+ 2⁴)modp，将Q加/减到具有相应权重(2²²⁸、2¹⁰⁰、2⁶⁸、2⁴)的系数上。我们先对 其进行约减的原因是可以使得元组<A,B>的低4个寄存器系数在进行约减时，低位 传到a₈,a₉的进位不会使得其字长太长(小于等于28位)，这样约减算法只需要进 行一轮约减，否则则需要进行两轮约减，这样可以提升算法的效率；

(3)与快速约减算法类似，图9中的第8-13行，利用AVX2的AND、SHR指令将元 组<A,B>的每个寄存器<A,B>_i,i∈[0,4]划分成三个部分L_i，M_i，H_i，其中L_i表示 AVX2寄存器<A,B>_i里的4个系数的低位部分，M_i表示AVX2寄存器<A,B>_i的 4个系数的中间部分，H_i表示AVX2寄存器<A,B>_i的4个系数的高位部分；

(4)然后，图9的第15行，利用AVX2的ALIGN函数，将寄存器M_i和M_i-1进行重组 成

记作M′_i，其中

分别表示系数a_2i,a_2i-1的中间 部分

分别表示系数b_2i,b_2i-1的中间部分，并将M′_i、L_i和H_i-1进行累加得 到新<A,B>_i,i∈[0,4]，至此<A,B>_i充分地被约减，且最大字长小于等于28位， 这个字长保证了在本发明的所有素数域运算中是不会溢出的；

实现效率：利用冗余表示法、交错元组和AVX2实现的固定时间的SM2素数域运算的时 钟周期对比表2所示，本发明实现的双向并行素数域运算的速度均比单向素数域运算快，具 体地，加法和减法运算(add/sub)快1.2倍，乘法运算(mul)快1.52倍，平方运算(sqr)快1.60倍， 系数约减运算(coeff.red.)则快1.53倍。由此可见AVX2的并行性可切实地为SM2素数域 运算提供效率上的提升，而且也为Co-Z运算和基于Co-Z运算的蒙哥马利梯形算法的效率提 升作出重要贡献。在不增加存储空间的前提下，本发明利用AVX2的高效和并行性，为基于 Co-Z运算的蒙哥马利梯形算法的安全高效实现提供了有效可行的方案。本发明除了可以被 SM2应用之外，在对素数域运算进行一定的修改后，本发明还可以有效地应用到其它 Weierstrass椭圆曲线公钥密码算法标准中。

表2 双向并行和单向素数域运算的时钟周期(单位:clock cycles)比较

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟 悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明 的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (4)

1.一种基于Co-Z运算的蒙哥马利梯形算法的AVX2快速实现方法，其特征在于，具体步骤如下：

步骤1，利用Co-Z初始化倍点算法XYCZ_IDBL计算蒙哥马利梯形算法中的两个Z坐标相同的椭圆曲线点R₁＝2P＝(X₁，Y₁，Z)和R₀＝P＝(X₀，Y₀，Z)，其中P是蒙哥马利梯形算法中输入的随机椭圆曲线点，R₀是通过XYCZ_IDBL算法更新Z坐标后与P等价的椭圆曲线点，X₀、Y₀分别表示R₀的横、纵坐标，R₁则是通过XYCZ_IDBL算法得到的与2P等价的椭圆曲线点，X₁、Y₁分别表示R₁的横、纵坐标；

步骤2，通过椭圆曲线点R₁和R₀的横坐标X₁、X₀预计算大整数A′＝(X₀-X₁)²和T′＝(X₀-X₁)A′；

步骤3，将椭圆曲线点R₁、R₀以及大整数A′、T′转换成字长为26位的冗余表示法，并存储到AVX2的256位寄存器中；

步骤4，进入蒙哥马利梯形算法的主循环，主循环的索引为i∈[1，254]，i从标量k的第254个比特开始往下遍历k的每个比特，每次循环获取k的第i个比特k_i和第(i+1)个比特k_i+1，并令符号a＝(k_i+k_i+1)mod 2，a∈{0，1}，通过a调整R₁、R₀的次序R_a、R_1-a，与A′、T′同时作为双向并行Co-Z Ladder算法的输入，并将双向并行Co-Z Ladder算法的输出分别更新到R_a、R_1-a、A′、T′中，循环调用双向并行Co-Z Ladder算法并不断更新R_a、R_1-a、A′、T′直到循环结束；

步骤5，主循环结束后，获取k最低位的两个比特k₀、k₁，并令符号a＝(k₀+k₁)mod 2，，通过a调整R₁、R₀的次序R_a、R_1-a，与A′、T′同时作为双向并行Co-Z共轭加法的输入，调用双向并行Co-Z共轭加法，并将双向并行Co-Z共轭加法的输出分别更新到R_a、R_1-a中；

步骤6，将R_1-a，R_a，P，a作为串行Co-Z取逆运算的输入进行取逆运算得到结果

其中Z＝x_PY_a(X₀-X₁)，λ＝y_PX_a，x_P、y_P表示P在仿射坐标下的横、纵坐标，X_a、Y_a表示R_a在射影坐标下的横、纵坐标，X₀、X₁则分别是R₀、R₁在射影坐标下的横坐标，λ是一个大整数；

步骤7，将R₀、R₁作为双向并行Co-Z加法的输入，调用双向并行Co-Z加法，并将双向并行Co-Z加法的结果更新到R₀、R₁中，最后通过坐标转换

求得蒙哥马利梯形算法仿射坐标形式的结果；

步骤8，将蒙哥马利梯形算法的结果从字长为26位的冗余表示法转换成字长为32位的多精度表示法。

2.如权利要求1所述的一种基于Co-Z运算的蒙哥马利梯形算法的AVX2快速实现方法，其特征在于，双向并行Co-Z加法的实现步骤如下：

(1)Co-Z加法运算的输入分别为椭圆曲线点P＝(X_P，Y_P)、椭圆曲线点Q＝(X_Q，Y_Q)；

(2)利用双向并行减法运算，同时计算T₁＝X_P-X_Q和T₂＝Y_Q-Y_P；

(3)对T₁、T₂利用双向并行平方运算计算其平方，分别得到A＝T₁ ²和D＝T₂ ²；

(4)利用双向并行乘法运算，同时计算B＝X_PA和C＝X_QA；

(5)利用双向并行加法运算，计算T₁＝B+C；

(6)利用双向并行减法运算，同时计算X_R＝D-T₁和T₃＝C-B；

(7)利用双向并行减法运算，计算T₁＝B-X_R；

(8)利用双向并行乘法运算，计算T₁＝T₁T₂和计算E＝Y_PT₃；

(9)利用双向并行减法运算，计算Y_R＝T₁-E；

(10)得到椭圆曲线点R＝P+Q＝(X_R，Y_R)，并更新P的坐标为(B，E)。

3.如权利要求1所述的一种基于C_o-Z运算的蒙哥马利梯形算法的AVX2快速实现方法，其特征在于，双向并行Co-Z共轭加法的实现步骤如下：

(1)Co-Z共轭加法的输入分别为椭圆曲线点P＝(X_P，Y_P)、椭圆曲线点Q＝(X_Q，Y_Q)，预计算的大整数A′＝(X_Q-X_P)²，T′＝(X_Q-X_P)A′；

(2)利用双向并行乘法运算，同时计算C＝X_QA′和E＝Y_PT′；

(3)利用双向并行减法运算，同时计算B＝C-T′和T₁＝Y_Q-Y_P；

(4)利用双向并行加法运算，同时计算T₂＝B+C和T₃＝Y_P+Y_Q；

(5)利用双向并行平方运算，同时计算

和

(6)利用双向并行减法运算，同时计算X_R＝D-T₂和X_R′＝F-T₂；

(7)利用双向并行乘法运算，同时计算T₃＝T₁T₂和T₃＝T₃T₄；

(8)利用双向并行减法运算，同时计算Y_R＝T₂-E和Y_R′＝T₃-E；

(9)椭圆曲线点R＝P+Q＝(X_R，Y_R)和R′＝P-Q＝(X_R′，Y_R′)。

4.如权利要求1所述的一种基于Co-Z运算的蒙哥马利梯形算法的AVX2快速实现方法，其特征在于，双向并行Co-Z Ladder算法的实现步骤如下：

(1)Co-Z Ladder算法的输入分别为椭圆曲线点P＝(X_P，Y_P)、椭圆曲线点Q＝(X_Q，Y_Q)，预计算的大整数A′＝(X_Q-X_P)²，T′＝(X_Q-X_P)A′；

(2)利用双向并行乘法运算，同时计算C′＝X_QA′和E′＝Y_PT′；

(3)利用双向并行减法运算，同时计算B′＝C′-T′和T₁＝Y_P-Y_Q；

(4)利用双向并行加法运算，同时计算T₂＝B′+C′和T₃＝Y_P+Y_Q；

(5)利用双向并行平方运算，同时计算

和

(6)利用双向并行减法运算，同时计算X_R＝D′-T₂和X_R′＝F′-T₂；

(7)利用双向并行减法运算，同时计算T₂＝B′-X_R和T₄＝X_R′-B′；

(8)利用双向并行乘法运算，同时计算T₂＝T₁T₂和T₄＝T₃T₄；

(9)利用双向并行减法运算，同时计算Y_R＝T₂-E′和Y_R′＝T₄-E′；

(10)利用双向并行减法运算，同时计算T₁＝X_R-X_R′和T₂＝Y_R-Y_R′；

(11)利用双向并行平方运算，同时计算

和

(12)利用双向并行乘法运算，同时计算X_P＝B＝X_RA和C＝X_R′A；

(13)利用双向并行加法运算，同时计算T₃＝T₂+B和T₄＝B+C；

(14)利用双向并行减法运算，同时计算X_Q＝D-T₄和T₁＝C-B；

(15)利用双向并行减法运算，同时计算T₄＝X_Q-X_P和T₃＝T₃-X_Q；

(16)利用双向并行平方运算，同时计算

和

(17)利用双向并行乘法运算，同时计算T′＝T₄A′和E＝Y_RT₁；

(18)利用双向并行加法运算，同时计算T₁＝D+A′和T₂＝E+E；

(19)利用双向并行减法运算，计算T₃＝T₃-T₁；

(20)利用双向并行减法运算和右移位运算，计算

(21)得到椭圆曲线点2P＝(X_P，Y_P)和P+Q＝(X_Q，Y_Q)。

Images