高基Montgomery模乘阵列结构设计与实现*

邬贵明，谢向辉，吴 东，郑 方，严忻恺

(数学工程与先进计算国家重点实验室， 江苏 无锡 214125)

高基Montgomery模乘阵列结构设计与实现*

邬贵明，谢向辉，吴 东，郑 方，严忻恺

(数学工程与先进计算国家重点实验室， 江苏 无锡 214125)

提出了两种高基Montgomery模乘线性阵列结构。两种线性阵列结构分别利用两种不同的并行化开发方法，沿不同的循环维度进行任务分配和调度，都能够充分开发算法的流水线并行。在Xilinx XC5VLX330 FPGA上实现了两种256位宽、基为216的模乘阵列结构。实验结果表明，两种结构具有84个时钟周期的延迟，吞吐率分别为1/17和1/21，与相关结构相比吞吐率更高。两种结构在性能和实现代价间能够达到合理平衡。

模乘；线性阵列结构；现场可编程阵列；流水化

1 引言

最常用的模乘算法是Montgomery模乘算法[1]，在软件实现和硬件实现两方面都非常高效。它使用简单的移位操作来代替更费时的除法操作，使用部分结果加模数倍数的方法来代替部分结果减模数倍数的方法，其好处是加法操作在被乘数的最低位处理后就可以立即开始，不用等到所有位处理完再开始，这一特性有利于开发并行性。

目前，学术界提出了大量Montgomery模乘硬件加速结构，这些结构可分为以下几类：基2模乘结构、多字基2模乘结构、高基模乘结构和采用大整数乘法器的模乘结构。1999年Tenca A F和 Koç Ç K[2]首次提出一种多字基2的Montgomery乘法算法和硬件实现，当操作数为n位时，一次Montgomery模乘需要2n个时钟周期。为降低计算延迟，大量研究人员基于Tenca A F和Koç Ç K的工作进行改进，其中Huang M[3]的工作最为突出，利用最高位的两种可能性预先对部分结果进行计算，将一次Montgomery模乘运行时间降为n个时钟周期，性能提升了1倍。基于经典的多字基2硬件结构，还有研究人员提出了多字基4的结构，但由于这些算法和结构的基比较低，计算延迟比较长，无法满足高性能结构的需求。在高基Montgomery模乘结构方面，McIvor C等[4]提出了高基Montgomery模乘算法和Systolic阵列结构，但该结构不能流水化连续处理多个模乘运算，仍然存在资源无法充分利用的问题。Zhou L[5]在Huang M[3]所提结构基础上扩展实现了高基模乘结构，由于关键路径上多了一个加法操作和一个乘法操作，导致运行性能较低。

本文提出了两种高基Montgomery模乘流水化阵列结构，两种结构分别在吞吐率和实现代价方面各有优势。

2 Montgomery算法

2.1 算法概要

给定两个整数M和R，R>M且R和M互素，R为2n。两个整数X和Y(X

步骤1t=XY；

步骤2m=t(-M-1modR) modR；

步骤3u=(t+mM)/R；

步骤4若u>M，则输出u-M，否则输出u。

模R和除R运算只需通过简单的截取和移位就能实现，避免了耗时的除法运算。

2.2 基2 Montgomery模乘

上述步骤的可实现性较差，算法1为实现性更强的基2 Montgomery模乘，算法中Y为被乘数，X为乘数，n为M的位数，y0为Y的最后一位，z0为Z的最后一位，⊕为逻辑异或操作。文献[3]给出了该算法的推导过程。

算法1基2 Montgomery模乘

Output:Z=XYR-1modM。

1:Z=0

2:fori=0ton-1do{

3: qi=(xiy0)⊕z0

4: Z=(Z+xiY+qiM)/2}

5:IfZ≥Mthen

6: Z=Z-M

BatinaL和WalterCD等人[6,7]证明，当X, Y<2M且R>4M时，XY2-n(modM) < 2M。算法中最后一步的减法可以避免，乘法结果可以直接作为下一个乘法的输入，注意这时的R已不再等于2n。这里假设R=2l，则l ≥ n+2，一般会选择l=n+2。算法2[8]给出了R =2n+2时不需要最后减法的Montgomery模乘算法。算法2中xn+1= 0，其输入X和Y、输出Z都在[0, 2M-1]范围内。

算法2不需要减法的基2Montgomery模乘

Output:Z=XYR-1mod2M。

1:Z=0

2:fori=0ton+1do{

3: qi=(xiy0)⊕z0

4: Z=(Z+xiY+qiM)/2}

2.3 多字基2Montgomery模乘

算法2中语句4的大整数Z、Y和M可以进一步切分成多个字进行表示，进一步减少实现代价。TencaAF和Koç ÇK[2]首次提出了多字基2的Montgomery模乘算法。

算法3中语句8将Z(e)赋为0，是考虑到Z(e)的最低位要赋给Z(e-1)的最高位，不能将i上次迭代的结果误传给这次迭代。可以进一步证明进位C只需两位就可以满足计算要求。该算法R =2n时，并不能处理输入X和Y都在[0, 2M-1]范围内的情况。

算法3多字基2Montgomery模乘

1:Z=0

2:fori=0ton-1do{

5:forj=1toedo{

8: Z(e)=0}

2.4 高基Montgomery模乘

高基Montgomery模乘可以进一步降低多字基2的Montgomery模乘的算法复杂度。算法4给出了高基Montgomery模乘算法，算法中字长为w，基为2w。由于R=2we> 4M，该算法能够处理X, Y<2M的输入。

算法4中的M′需要预先进行计算。我们将语句3和4定义成任务A(i)，语句6和7定义成任务B(i,j)，则该算法的数据相关图如图1所示。同一列上的两任务数据相关性由C引起，两列之间的两任务数据相关性由Z(j)引起。

算法4高基Montgomery模乘

1:Z=0

2:fori=0toe-1do{

3: qi=(Z(0)+X(i)Y(0))M′mod2w

4: (C,Z(0))=x(i)Y(0)+qiM(0)+Z(0)

5:forj=1toe-1do{

6: (C,Z(j))=C+X(i)Y(j)+qiM(j)+Z(j)

7: Z(j-1)=Z(j)}

8: Z(e-1)=C}

Figure 1 Data dependency graph of high radix Montgomery modular multiplication图1 高基Montgomery模乘数据相关图

3 素域高基Montgomery模乘阵列结构设计

3.1 并行性开发

高基Montgomery模乘算法的并行性可通过两种方法进行开发。如果任务按照i循环进行划分，i循环的迭代将由处理单元PE完成，即图1中每一列的计算将由一个PE完成，并行开发方法如图2a所示。在该方法中，每个PE中用同一个X(i)参与计算，不同时间完成对Z不同段的计算。

Figure 2 Two approaches of exploiting parallelism图2 两种开发并行的方法

如果任务按照j循环进行划分，j循环的迭代(即j相同时i不同的所有迭代)将由PE完成，并行开发方法如图2b所示。在这个方法中，每个PE中使用同一个Y(j)、M(j)和Z(j)参与计算，不同时间使用不同段的X(i)来参与对Z的同一段的更新。

3.2 素域高基Montgomery模乘阵列结构

本文提出两种线性阵列结构分别实现图2中两种并行开发方法的调度过程，实现第一种调度的结构称为结构1，实现第二种调度的结构称为结构2，如图3所示。

结构1中的每个处理单元结构相同，处理单元结构见图4。多个PE组织成流水线结构形式，每个PE接收前一个PE传递来的数据和结果，并把数据和自身的计算结果传递给下一个PE。PE间传递的是Y(j)、M(j)和Z(j)，其中Y(j)、M(j)是传递来的输入数据，而Z(j)是前一个PE的计算结果。

Figure 3 Two linear arrays图3 两种线性阵列结构

Figure 4 Processing element of architecture 1图4 结构1的处理单元

该PE结构不同时刻形成不同的数据通路来完成不同任务。在该结构中，任务A(i)需要三个节拍完成，第一拍完成P=(Z(0)+X(i)Y(0))运算，第二拍完成qi= PM′，第三拍才开始执行(C, Z(0))=X(i)Y(0)+ qiM(0)+ Z(0)，之后执行任务B(i,j)，需要一拍完成，这样PE之间的延迟为四拍，这是高基Montgomery模乘的缺点。多字基2Montgomery模乘算法的语句3是简单的逻辑运算，而高基Montgomery模乘算法的语句3变成了两个w字节的乘法和一个w字节的加法，计算延迟由忽略不计变成了两个时钟节拍。在TencaAF和Koç ÇK的结构中,PE间的延迟为2拍，由于使用了预先计算的方法，HuangM的结构变为1拍。而本文的结构不可能使用预先计算的方法，在HuangM的结构中，只有1位需要假设，而高基算法需要有w位进行假设，代价很大不可实现。

结构2中除第一个PE外，其它每个处理单元结构相同。PE0完成任务A(i)，而其它PE完成任务B(i,j)。PE间组织成流水线结构，PE间由前往后传递的是输入数据X(i)和计算得到的进位C，由后往前传递的是Z(j-1)。第一个PE用于计算qi和Z(0)，可以采用结构1的PE结构，而为了流水化连续处理多个模乘，在原来的结构上增加一个乘法器并修正原来的计算通路；其它PE仅完成任务B(i,j)，原来的PE结构用于完成任务A(i)和B(i,j)，因此结构在原来的基础上可以得到进一步精简。

PE0仅完成任务A(i)，需要三个节拍完成。前一个PE任务完成，后一个PE才能开始计算。第二个PE仅完成任务B(i,j)，计算延迟为一拍，因此第三个PE在第五个节拍才能开始计算。除第一个PE外，PEj在第3+j个节拍开始计算。若第一个PE不能流水化处理多个模乘，则第一个PE在第五拍时才能重新开始计算，与第一次计算相隔四拍。这样会导致阵列中所有PE每四拍才运算一次，计算效率非常低。实际上，可充分利用所有PE空闲的三拍进行其它三个模乘的运算。另外，对于最后一个PE，输出C直接连接到输入Z(j-1)上。

下面讨论如何在结构2上使其能够流水化处理四个模乘。每个处理单元的结构几乎不用改动，主要需要考虑如何对来自四个不同输入的Y(j)、M(j)和Z(j)以及输出Z(j-1)进行有效的调度。这里采用循环移位器来实现对四个不同输入及中间结果的调度。图5以Y为例说明了四组输入(A、B、C和D)的调度情况，Y的每一段在不同时刻进入循环移位器，四个移位器中第一个便连接在各自的PE输入上。输出Z(j-1)的调度更简单点，第二个PE的输出要送给第一个PE进行计算，而该计算在下一拍便开始，因此输出直接连接在第一个PE的输入Z(j-1)上。第三个及以后的PE的输出要在三拍之后才能用上，因此只需在输出和前一个PE的输入Z(j-1)间采用两个移位寄存器即可。

Figure 5 Scheduling of four Y inputs图5 四个Y输入的调度

4 实验结果

两种线性阵列结构都采用Verilog HDL进行描述，以Xilinx XC5VLX330为目标器件进行综合，基于ModelSim 6.5d工具进行实验评测，综合结果由Xilinx ISE 13.3给出。

对于结构1，PE的计算核心是(C,Z(j))=C+X(i)Y(j)+qiM(j)+Z(j)，该计算的关键路径在于C的反馈，为降低关键路径的计算延迟，可以使X(i)Y(j)+qiM(j)+Z(j)在C产生之前提前计算，得到A=X(i)Y(j)+qiM(j)+Z(j)，这样关键路径上只有一个加法运算(C,Z(j))=C+A。另外，中间结果Z还可以使用进位保留的形式，这样可以进一步加速加法运算并节省资源。这里为了简单化实现，我们没有采用提前计算和进位保留的方法。对于结构2，PE(除第一个PE)的计算核心仍然是(C,Z(j))=C+X(i)Y(j)+qiM(j)+Z(j)，但关键路径已不在于C，C不参与反馈，直接传递给下一个PE，关键路径转移到了X(i)到C和qi到C。

另外，两种结构中的乘法去处均采用了FPGA中的DSP模块实现，若采用进位保留的乘法器，将有利于性能的提升。在时序约束为200 MHz时，针对模数位宽n=256、字长w=16、字数e=17的情况，所设计的两种结构的综合结果如表1所示。

Table 1 Synthesis result when targeting Xilinx XC5VLX330表1 目标器件为Xilinx XC5VLX330时的综合结果

表1括号中的数据表示所占FPGA资源的百分比。从表1可以看出,结构1的运行频率低于结构2，实现代价高于结构2。两种结构的延迟虽然一致，结构2的实现代价却低于结构1的1/2，吞吐率仅降低了20%。

表2比较了本文提出的结构与其它两种典型结构的性能，设计参数基本相同。文献[4]的结构为高基Montgomery模乘结构，延迟为38拍，优于本文的结构，然而吞吐率却不高；文献[8]的结构为位级Systolic结构，需要更多的节拍数。这两种典型结构无法处理输入X和Y在[0, 2M-1]的情况，而本文的结构可以处理。

Table 2 Performance comparison表2 性能比较

5 结束语

本文提出了两种素域高基Montgomery模乘流水化阵列结构，两种结构分别在吞吐率和硬件实现代价方面各有优势。相对于长整数乘法，所提出的结构延迟太长，下一步工作将研究采用大整数乘法部件实现Montgomery模乘结构。

[1] Montgomery P L. Modular multiplication without trial division[J]. Mathematics of Computation, 1985, 44(170):519-521.

[2] Tenca A F, Koç Ç K. A scalable architecture for Montgomery multiplication[C]∥Proc of Cryptographic Hardware and Embedded Systems, 1999:94-108.

[3] Huang M, Gaj K, El-Ghazawi T. New hardware architectures for Montgomery modular multiplication algorithm[J]. IEEE Transactions on Computers, 2011, 60(7):923-936.

[4] McIvor C, McLoone M, McCanny J V. High-radix systolic modular multiplication on reconfigurable hardware[C]∥Proc of IEEE International Conference on Field-Programmable Technology, 2005:13-18.

[5] Zhou L, Huang M, Smith S C. High-performance and area-efficient hardware design for radix 2k Montgomery multipliers[C]∥Proc of International Conference on Computer Design, 2011:1.

[6] Batina L, Muurling G. Montgomery in practice:How to do it more efficiently in hardware[C]∥Proc of Topics in Cryptology, the Cryptographer’s Track at the RSA Conference, 2002:40-52.

[7] Walter C D. Montgomery’s multiplication technique:How to make it smaller and faster[C]∥Proc of Cryptographic Hardware and Embedded Systems, 1999:80-93.

[8] Örs S B, Batina L, Preneel B, et al. Hardware implementation of elliptic curve processor over GF(p)[C]∥Proc of IEEE International Conference on Application-Specific Systems, Architectures and Processors, 2003:433-443.

WUGui-ming,born in 1981,post doctor，engineer,his research interests include high-performance computer architecture, and reconfigurable computing.

DesignandimplementationofhighradixMontgomerymodularmultiplicationarraystructures

WU Gui-ming，XIE Xiang-hui，WU Dong，ZHENG Fang，YAN Xin-kai

(State Key Laboratory of Mathematical Engineering and Advanced Computing,Wuxi 214125,China)

Two linear arrays for high radix Montgomery modular multiplication are proposed. They use two different parallelization methods, both of which can exploit pipelined parallelism through task assignment and task scheduling along different loop dimensions. The two linear arrays for 256-bit modular multiplication using the radix of 216, are implemented on Xilinx XC5VLX330 FPGA. The experimental results show that both linear arrays have the latencies of 84 cycles, and the throughput of 1/17 and 1/21, respectively. Compared with the related work, our designs have higher throughput. Moreover, the balance between performance and hardware overhead can be achieved.

modular multiplication;linear array;FPGA;pipeline

2013-08-05;

：2013-10-26

中国博士后科学基金资助项目；国家863计划资助项目(2013AA010105)

1007-130X(2014)02-0201-05

TP302

：A

10.3969/j.issn.1007-130X.2014.02.002

邬贵明(1981-),男,四川隆昌人，博士后，工程师，研究方向为高性能计算机体系结构和可重构计算。E-mail:wu.guiming@meac-skl.cn

通信地址：214125 江苏省无锡市滨湖区楝泽路30号山水城软件园数学工程与先进计算国家重点实验室Address:State Key Laboratory of Mathematical Engineering and Advanced Computing,Shanshuicheng Software Park,30 Lianze Rd,Binhu District,Wuxi 214125,Jiangsu,P.R.China