CN115348002A - 一种基于多字长乘法指令的Montgomery模乘快速计算方法 - Google Patents

Abstract

本发明涉及公钥密码体制领域，尤其涉及一种基于多字长乘法指令的Montgomery模乘快速计算方法。所述方法包括：确定Montgomery模乘的输入参数；将Montgomery模乘的多精度乘法运算和约减运算分别拆分为三类分块；采用分块集成的方式，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，得到Montgomery模乘结果。本发明充分利用多字长大数乘法的优势以及可用寄存器的数量，可适应不同操作数长度的多字长指令，适用于不同模数长度的Montgomery模乘，具有很强的适应性；相对于其他复杂的分块方法以达到优化目的而言，本发明分块方法简洁但优化效果明显。

Description

一种基于多字长乘法指令的Montgomery模乘快速计算方法

技术领域

本发明涉及公钥密码体制领域，具体地说，尤其涉及一种基于多字长乘法指令的Montgomery模乘快速计算方法。

背景技术

公钥密码体系中目前使用最广泛的密码算法是RSA和ECC算法，其最重要的运算是有限域上的大数模乘。为了提高ECC和RSA的计算效率，有限域上大数模乘的优化成为研究的重点。1985年，P.L.Montgomery提出了Montgomery模乘算法，它是现在应用最广泛的一种模乘算法。同时蒙哥马利模乘器是也调用频率最高最耗时的模块，作为公钥加密系统和众多通信系统中的核心运算，其运算速度决定了系统的性能。其基本思想是用简单省时的加法和移位操作代替费时的求逆和除法操作。对于计算模乘P＝A*B mod N(A,B,N为m比特的二进制大整数且0<A,B<N)，该算法首先选取一个与N互素的整数R(一般取R＝2^m)，将模N的乘法运算转化为模R的乘法运算，而模R乘法可用移位代替，避免模N乘法中的除法和求逆运算。Montgomery模乘算法的计算结构见图1所示，Montgomery域变换的计算方法为A’＝A*Rmod N。定义Montgomery乘积：MonPro(A,B)＝A*B*R^-1mod N。具体分为以下三个步骤：

(1)计算A，B的N剩余类：A’＝A*R mod N＝A*R²*R^-1mod N＝MonPro(A，R²),

B’＝B*R mod N＝B*R²*R^-1mod N＝MonPro(B，R²)；

(2)计算A’和B’的Montgomery乘积P’＝A’*B’*R^-1mod N＝MonPro(A’,B’)；

(3)将P’转化为模乘乘积P，P＝A*B mod N＝A’*R^-1*B’*R^-1mod N

＝P’*R^-1mod N＝MonPro(P’，1)。

因此，Montgomery模乘算法的核心是计算Montgomery乘积MonPro(A’,B’)，其具体算法流程描述见图2所示。从图2中可以看出Montgomery乘积计算可以分为两个关键步骤：计算Montgomery多精度乘法C←A*B和计算Montgomery约减P←(C+M*N)/R。

目前，有许多研究都集中在如何实现并优化Montgomery模乘算法上。国内外也相继提出了各式的优化方法，比较集中的思路是通过合理利用流水线和硬件资源，设计Montgomery模乘操作数各字长的计算顺序。比如由胡世文等2019年提出的流水线型蒙哥马利计算方法就是将模乘中的模数R进行拆分循环以减少硬件资源的使用。还有由杨帆2019年提出的蒙哥马利不均等分块法，通过对A,B,N等操作数添加多余的位数来满足分块需求，从而进行小分块规划来减少计算步骤，但该方法有添加位的操作可能会带来新的计算量。而2020年由李杨等提出的CIPOHS方法也是利用不均等分块的思想计算Montgomery模乘，该算法将混合扫描思想应用于Montgomery模乘计算中，通过动态选取d合理利用操作数，减少嵌入式系统中内存存取数量，并与之前相关算法SOS、CIOS、FIOS、FIPS和CIHS的存储性能优劣做了对比。而本发明则是在此基础上结合多字长乘法的特点提出一种新的分块方法，进一步优化Montgomery模乘计算的存储性能及硬件资源的使用。

发明内容

本发明的目的在于克服现有技术缺陷，提出了一种基于多字长乘法指令的Montgomery模乘快速计算方法。

为了实现上述目的，本发明提出了一种基于多字长乘法指令的Montgomery模乘快速计算方法，所述方法包括：

确定Montgomery模乘的输入参数；

将Montgomery模乘的多精度乘法运算和约减运算分别拆分为三类分块；

采用分块集成的方式，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，得到Montgomery模乘结果。

作为上述方法的一种改进，所述确定Montgomery模乘的输入参数；具体包括：

设模N为m比特，Montgomery的系数为R＝2^m，A和B为模乘运算的乘数和被乘数且A和B均为m比特，0<A,B<N，多字长乘法指令操作数T为U比特，处理器的字长为W比特，预计算常数为N’，N’为N的模R负逆，并且满足N×N’＝-1mod R，mod表示取模运算；

计算N、A和B的字数大小均为

表示取大于等于其的最小整数运算，T的字数大小为

且U≥192比特，将数A表示为：A＝(A[n-1],…,A[2],A[1],A[0])，将数B表示为：B＝(B[n-1],…,B[2],B[1],B[0])，乘积C＝A×B表示为：C＝(C[2n-1],…,C[2],C[1],C[0])，T表示为T＝(T[t-1],…,T[2],T[1],T[0])，计算参数r＝n％t，％为取余数运算符，令d＝r+t，令s为(n-d)/t的整数部分。

作为上述方法的一种改进，所述将Montgomery模乘的多精度乘法运算和约减运算分别拆分为三类分块；具体为：

将Montgomery模乘的多精度乘法运算A×B表示成第一菱形图，所述第一菱形图的每个点表示一个乘积A[i]×B[j]，0≤i,j<n；

将约减运算M×N表示成第二菱形图，所述第二菱形图的每个点表示一个乘积M[i]×N[j]，0≤i,j<n；M[i]为约减运算的参数，根据下式计算得到：

C[x]＝A[i]×B[j],x∈[0,2n-1],i∈[0,n-1],j∈[0,n-1]

M[i]＝C[x]×N’[i],x∈[0,2n-1],i∈[0,n-1]

将第一菱形图拆分成三类分块，所述三类分块每类由若干个形状相同的小块组成，第二菱形图的拆分方式与第一菱形图相同。

作为上述方法的一种改进，所述三类分块包括第一类分块、第二类分块和第三类分块；具体为：

所述第一类分块为边长为d的菱形，其中d＝r+k×t，k为正整数，根据处理器寄存器数量V选择大小，d≤V且(n-d)能被(k×t)整除，从参与模乘运算的乘数A和被乘数B的最低位开始，用于处理n/t非整除部分，所述第一类分块包括1个小块；

所述第二类分块为边长分别为d和k×t的平行四边形，位置分布为沿着第一类分块的两边，以d为底边长，分别以k×t为步长延伸到最高位，所述第二类分块包括2×(n-d)/k×t个小块；

所述第三类分块为边长为k×t的菱形，是除了第一类分块和第二类分块后余下的部分，所述第三类分块包括[(n-d)/(k×t)]²个小块。

作为上述方法的一种改进，所述采用分块集成的方式，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，得到Montgomery模乘结果；具体为：

对于第一类分块的小块计算，非整除部分使用单字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，整除部分使用多字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，得到Montgomery模乘结果；

对于第二类分块的小块和第三类分块的小块，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，同时得到这两类的Montgomery模乘结果。

作为上述方法的一种改进，所述对于第一类分块的小块计算，非整除部分使用单字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，整除部分使用多字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，得到Montgomery模乘结果；具体为：

对于非整除部分，长度为r，使用单字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，其中，

多精度乘法运算为：

{C[0],C[1]....C[r-1]}＝{C[0],C[1]...C[r-1]}+B[i]×{A[0],A[1]....A[r-1]}；

约减运算为：

{M[0],M[1]....M[r-1]}＝{C[0],C[1]....C[r-1]}×N’mod R；

P’＝{C[0],C[1]....C[r-1]}+{M[0],M[1]....M[r-1]}×N[i]，

其中i为循环控制参数，0≤i<d；

对于整除部分，使用多字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，其中，

多精度乘法运算为：

{C[r],C[r+1]....C[r+t-1]}＝{C[r],C[r+1]...C[r+t-1]}+B[i]×{A[r],A[1+r]....A[t+r-1]}；其中，B[i]的循环遵循操作数扫描，i为循环控制参数，0≤i<d；

约减运算为：

{M[r],M[r+1]...M[r+t-1]}＝{C[r],C[r+1]....C[r+t-1]}×N’mod R，

P’＝{C[r],C[r+1]....C[r+t-1]}+N[i]×{M[r],M[r+1]...M[r+t-1]}。

作为上述方法的一种改进，所述对于第二类分块的小块和第三类分块的小块，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，同时得到这两类的Montgomery模乘结果；具体为：

对于第二类分块的小块和第三类分块的小块，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，分为外循环和内循环；其中，

第二类分块的外循环控制参数为j，0≤j<s，用于计算第二类分块的相邻下一块，在j的一个外循环体内，依次进行d次内循环的多精度乘法运算，再依次进行d次内循环的约减运算；i为内循环的循环控制参数，0≤i<d，用于控制每个小块；直至外循环重复s次，得到第二类分块的Montgomery模乘结果；

第三类分块的外循环控制参数为j，0≤j<s，用于计算第三类分块的相邻下一块，在j的一个外循环体内，包括两层循环嵌套，中循环控制参数为u,0≤u<j，在u的一个循环体内，依次进行v次内循环的多精度乘法运算，再依次进行v次内循环的约减运算；v为内循环的循环控制参数，0≤v<t，用于控制每个小块；中循环重复j次后，外循环控制参数j加1，再重复中循环和内循环；直至外循环重复s次，得到第三类分块的Montgomery模乘结果。

与现有技术相比，本发明的优势在于：

1、本发明充分利用多字长大数乘法的优势以及可用寄存器的数量，可适应不同操作数长度的多字长指令，适用于不同模数长度的Montgomery模乘，具有很强的适应性；

2、本发明以减少中间结果存储、减少计算步骤为目的设计分块计算方式，最大限度地利用多字长乘法指令单次运算长度优势，将Montgomery模乘的多精度乘法运算和约减运算分别拆分为三类分块进行计算，相对于其他复杂的分块方法以达到优化目的而言，本发明分块方法简洁但优化效果明显；

3、本发明根据算法的需求重新设计扫描方式，结合算法分块类型设计适合算法的扫描方式，尽可能共享操作数以及节省存取次数，整体上减少了内存存取次数，同时减少乘法指令执行次数，提高了Montgomery模乘算法实现效率。

附图说明

图1是现有的Montgomery模乘计算结构示意图；

图2是现有的Montgomery模乘计算的流程图；

图3是本发明基于多字长乘法指令的Montgomery模乘快速计算方法的流程图；

图4是本发明基于多字长乘法指令的Montgomery模乘快速计算方法的结构示意图。

具体实施方式

Montgomery模乘算法的实现都可以由一系列操作：乘法mul，加法add，加载load，存储store等实现。高性能的实现算法主要集中在优化这些操作上面。在系统中，由于可用寄存器的数量有限，load/store等存储操作显得尤为重要。本发明的目的在于，为了充分利用部分CPU(如Cavium Octeon系列MIPS64 CPU)大数乘法指令提供的优良性能，提高Montgomery模乘算法的计算效率，提出了一种利用多字长大数乘法指令设计分块计算Montgomery模乘算法中的多精度乘法C←A*B和Montgomery约减P←(C+M*N)/R的优化算法MWBCHS(Multi-Words Block Coarsely Integrated Product and Operand HybridScanning)，充分利用多字长大数乘法指令的优势，合理安排分块运算顺序，用混合扫描计算的方式计算Montgomery模乘，减少指令执行次数，同时充分利用处理器的可用寄存器数量，通过在操作数扫描时共享操作数减少操作数存取数目，整体上减少了内存存取次数，提高了算法实现效率。

下面结合附图对本发明的技术方案进行详细的说明。

对于多字长算法的设计，首先要考虑的是操作数字数长是否为多字长乘法被乘数长度t的整数倍。若为非整数时，怎样处理才能更好的利用多字长乘法的优势以减少存取次数。实际设计分块的时候必须依照t的大小来确定分块的最小边长，其至少有一边不小于t的大小。且为了充分发挥该指令优势，减少存取次数和指令调用次数，分块时应尽可能多的使分块边长有一边为多字长乘法被乘数字数长t的倍数。另外该算法分块长度与操作数字数长n是否是多字长乘法操作数字数长t的整数倍有关，但区别仅在第一类分块的非整除部分，计算时若没有非整除部分，不计算该部分即可。故以下描述时不单独描述整除情况，以(n％t)作为描述参量，整除时取0即可。下面描述具体的操作步骤(以下步骤中的分块参数t在处理器寄存器数量充足时均可替换为k*t，k为正整数)：

1)由需快速计算的算法(如椭圆曲线密码算法ECC、RSA、国密SM2等)确定Montgomery大数模乘的输入参数，包括：m比特模N，Montgomery系数R＝2^m，2个m比特操作数A、B(A,B分别代表模乘运算的乘数和被乘数且0<A,B<N)，预计算常数N’(满足NN’＝-1(modR)，N’为N的模R负逆)。还有两个参数：多字长乘法指令操作数(被乘数T)比特长度U，处理器的字长W比特。将上述参数转换为以处理器字长整数倍表示的数。N、A、B的字数大小是

T的字数大小是

数A、B可以被表示为：A＝(A[n-1],…,A[2],A[1],A[0])，B＝(B[n-1],…,B[2],B[1],B[0])，乘积C＝A*B可以表示为：C＝(C[2n-1],…,C[2],C[1],C[0])，被乘数T可以表示为T＝(T[t-1],…,T[2],T[1],T[0])，计算参数r＝n％t，d＝r+t，s＝[(n-d)/t]，其中r表示A、B字数n除以多字长乘法指令操作数被乘数字数t的余数部分；

2)Montgomery模乘的多精度乘法运算A*B和约减运算M*N表示成菱形图，菱形结构的每个点表示一个部分乘积A[i]·B[j]，最右边的点表示运算字数的最低位的部分乘积A[0]·B[0]，最左边表示最高位部分乘积A[n-1]·B[n-1]，M[i]·N[j]同理。根据1)计算的参数r、d和s，将菱形拆分成三类分块。三类分块分别为：第一类分块是边长为d的菱形块，从参与模乘运算的乘数和被乘数的最低位开始，用于处理n/t非整除部分(若n/t可整除，则d＝t)，小块个数为1；第二类分块是边长为d和t的平行四边形块，位置分布为沿着第一类分块的两边，以d为底边长，分别以t为步长延伸到最高位，小块个数为2*(n-d)/t；第三类分块是边长为t的菱形块，是划分了除第一类和第二类分块后剩下未划分的边长为(n-d)的大菱形块，小分块个数为[(n-d)/t]²；

3)第一类小块计算：处理长度为r的部分，使用常规单字长乘法指令进行运算，计算多精度乘法部分：

{C[0],C[1]....C[r-1]}＝{C[0],C[1]...C[r-1]}+B[i]*{A[0],A[1]....A[r-1]}和约减部分{M[0],M[1]....M[r-1]}＝{C[0],C[1]....C[r-1]}*N’mod R与

P’＝{C[0],C[1]....C[r-1]}+{M[0],M[1]....M[r-1]}*N[i]；

使用多字长乘法指令，计算长度为t的部分：多精度乘法部分

{C[r],C[r+1]....C[r+t-1]}＝{C[r],C[r+1]...C[r+t-1]}+B[i]*{A[r],A[1+r]....A[t+r-1]}，约减部分{M[r],M[r+1]...M[r+t-1]}＝{C[r],C[r+1]....C[r+t-1]}*N’mod R，

P’＝{C[r],C[r+1]....C[r+t-1]}+N[i]*{M[r],M[r+1]...M[r+t-1]}，其中B[i]的循环遵循操作数扫描OS(Operands scan)。其中i为循环控制参数，0≤i<d，i++；

4)第二类小块与第三类小块的计算：采用多字长大数乘法指令，分为内循环和外循环，内循环为以下a,b,c,d,e步骤顺次计算，其中j为外循环参数，0≤j<s，在a,b,c,d,e步骤中j保持不变。

(a)计算第二类小块的多精度乘法部分(对应图4中标号7，8)：

{C[z],C[z+1]....C[z+t-1]}＝{C[z],C[z+1]...C[z+t-1]}+B[i]*{A[t*j+d],A[t*j+1+d]....A[t*j+t+d-1]}和

{C[z],C[z+1]....C[z+t-1]}＝{C[z],C[z+1]...C[z+t-1]}+A[i]*{B[t*j+d],B[t*j+1+d]....B[t*j+t+d-1]}，其中B[i]的循环遵循操作数扫描，其中：0≤i<d，循环次数为d次，i为循环变量，i++，z＝i+t*j+d；

(b)计算第二类小块的约减部分(对应图4中标号9.10)：

{M[d+t*j],M[d+t*j+1]...M[d+t*j+t-1]}＝{C[d+t*j],C[d+t*j+1]...C[d+t*j+t-1]}*N’mod R，P’＝N[i]*{M[d+t*j],M[d+t*j+1]...M[d+t*j+t-1]}及

P’＝{C[d+t*j],C[d+t*j+1]...C[d+t*j+t-1]}+M[i]*{N[d+t*j],N[d+t*j+1]...N[d+t*j+t-1]}，其中：0≤i<d，循环次数为d次，i为循环变量，i++，z＝i+t*j+d；

(c)计算第三类分块，分内外两层循环，本步骤和后续d步骤为内层循环：计算单个小块的多精度乘法部分(对应图4中标号11)：

{C[z],C[z+1]...C[z+t-1]}＝

{C[z],C[z+1]...C[z+t-1]}+B[d+v+t*i]*{A[t*i+d],A[t*i+1+d]...A[t*i+t+d-1]}，

{C[z],C[z+1]...C[z+t-1]}＝

{C[z],C[z+1]...C[z+t-1]}+A[d+v+t*i]*{B[t*i+d],B[t*i+1+d]...B[t*i+t+d-1]}，其中，0≤v<t，v为循环控制变量，循环次数为t次，v++，0≤u<j，u为c步骤和后续d步骤第三类分块外层循环共同控制变量，c)，d)步骤中i保持不变，z＝v+2t*u+2d；

(d)计算第三类单个小块的约减部分(对应图4中标号12)：

P’＝N[k+t*u+d]*{M[t*u+d],M[t*u+1+d]...M[t*u+t+d-1]}，

P’＝{C[d+t*u],C[d+t*u+1]...C[d+t*u+t-1]}+

M[k+t*u+d]*{N[t*u+d],N[t*u+1+d]...N[d+t+t*u-1]}，其中，0≤v<t，0≤u<j，循环次数为t次，v++，i同步骤c，z＝v+2t*u+2d；

(e)计算第三类分块的外层循环：对步骤(c),(d)顺次循环执行u次。

(a),(b),(c),(d)，(e)顺次执行s次循环(外循环)。之后，P＝P’/R，如果P≥N则P←P-N；返回结果P。算法整体流程图参照说明书附图3。

为了清晰的表现出算法的结构，同时给出伪代码的方式对算法进行描述：

首先算法1给出Montgomery算法流程及参数，算法2给出细化伪代码。

算法1.计算Montgomery乘积算法流程输入：m比特模N；Montgomery系数R＝2^m；2个操作数0<A,B<N；

预计算常数N’，以满足NN’＝-1(modR)，N’为N的模R负逆。

输出：Montgomery乘积P＝MonPro(A,B)＝A*B*R^-1mod N.

1)C←A*B；//采用普通乘法指令进行多精度计算

2)M←C*N’mod R；

3)P←(C+M*N)/R；

4)If P≥N then P←P-N；

5)Return P

算法2.多字长乘法计算Montgomery乘积伪代码：

对于第一类分块的小块计算，非整除部分使用单字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，整除部分使用多字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，得到Montgomery模乘结果；

对于第二类分块的小块和第三类分块的小块，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，同时得到这两类的Montgomery模乘结果。

输入：m比特模N；Montgomery系数R＝2^m；2个操作数0<A,B<N；

预计算常数N’，以满足NN’＝-1(modR)；

输出：Montgomery乘积P＝MonPro(A,B)＝(A*B*R)^-1mod N。

1)

//(m为大数的比特位长度，W是处理器的字的比特位长)。

2)C←0；//C为A*B的值，C[z]为A[i]*B[j]计算C[i+j]的中间值,(例C[2]＝A[0]*B[2]+A[1]*B[1]+A[2]*B[0],C[z]＝A[0]*B[2]或A[0]*B[2]+A[1]*B[1]等中间运算结果),C[]为计算完成C值([]中为具体参数值)。

3)r＝n％t；d＝n％t+t，s＝[(n-d)/t]；//r表示n除以t的余数部分

4)操作数中用{}表示该步骤使用多字长乘法，一次运算多个字长。以下为三类分类计算过程的伪代码：

为了进一步说明本算法优势，以图4为例我们分析本专利优化算法内存存取的数目的优势。图4上面的菱形表示计算多精度乘法C←A*B，下面的表示计算Montgomery约减P←(C+M*N)/R。对于上面的多精度乘法，菱形结构使用灰色表示，对于下面的Montgomery约减，菱形结构使用深灰色表示。以下算法设计多精度大数乘法指令以Marvell/Cavium公司基于MIPS64扩展的v3mulu指令为例，该指令每次可完成192*64bits及以上乘法运算。图4算法操作数为n＝8即512bit*512bit(W＝64bit)乘法。由图示菱形图中每一个箭头经过的点代表一次运算的中间运算结果，箭头的方向表示操作数扫描顺序(单字长和多字长顺序相同，每次运算区域均为从右至左运算)。每一个小模块编号代表计算的顺序，每个色块所包涵的点代表一次连续运算所计算的点。

本算法分块分三大类，如图4(图4示例中以t＝3为例)，右下角深灰色和浅灰色的小菱形框住部分为第一类非整除分块，分块大小为(n％t+t)*(n％t+t)，采用了少量的单字长乘法以配合计算(n/t)非整除部分，整除部分用多字长乘法计算如图4的5，6分块。第二类为右边和下边沿边长为(n％t+t)*t的边长为多字长乘法指令运算长度及第一类分块长度的平行四边形，如图4中的7，8，9，10分块。最后一部分是左上角(n-n％t-t)*(n-n％t-t)的规则菱形且边长为多字长乘法指令运算长度分块，为第三类分块，该分块是边长为t的整数倍的菱形块，计算时分成t*t的小分块进行计算，如图4中的11，12分块。

具体分析算法的存取数量，按上述描述分三类分别阐述，为使描述简洁其中(n％t+t)参量用参数d代替。同时为具体对比存储数量以下v3mulu指令的操作数的字数长t＝3，算法MWBCHS分块和存储数量见表1。

首先分析第一类右下角分块，右下角分块单指令部分运算所需要的LD次数为2*d+2，即d个寄存器存储A，一个寄存器依次存储B，加载两次B。上面一小块多字长运算所需要的LD次数为3，因为连续计算A不用重新加载，只加载B即可。所以第一类分块的总LD数量为2d+5。

第二类为右边和下边d*t分块，t为3，每个小分块所需LD为(3+d)。共有2*((n-d)/3)个第二类小块，所以第二类小分块总的LD次数为2*((n-d)/3)*(3+d)。

第三类分块是t*t分块，共有((n-d)/3)²个，每个小块需要6次LD。所以第三类分块总的LD数量为6*((n-d)/3)²次。

所以本算法总LD的数量为[6*((n-d)/3)²+2*((n-d)/3)*(3+d)+2d+5]，即[(2n²-2nd+6n+15)/3]，而在整个算法中需要存储的是n个字数长的M(M←C*n₀’mod r)和n+1个字数长的最终结果C，故store的总量是n+n+1＝2n+1；所以内存存取(load和store)的总量为[(2n²-2nd+6n+15)/3]+2n+1。

目前业内已有的几种主要的改进Montgomery算法：SOS、CIOS、FIOS、FIPS、CIHS和CIPOHS的存储对比如下表1：

表1本方法和几种常用Montgomery模乘方法内存存取复杂度对比分析

几种算法中CIPOHS是目前存储性能最优的算法。下面将本算法的存储性能和CIPOHS算法进行对比，本算法借鉴了李扬等在一种适用于嵌入式系统的Montgomery模乘计算方法CIPOHS中的粗粒度扫描，但分块方式则结合多字长乘法进行了优化，存储性能和运算指令执行次数均有明显提升，见表2：

表2本方法和CIPOHS算法的性能对比

由表2可知，本算法存储性能明显优于CIPOHS算法，且当操作数n值越大，本算法相对CIPOHS算法运算次数明显减少，节省存储空间，存储性能更好。同时本算法所用寄存器第一类分块最多为d+1+2d+2＝3d+2个，其他小块计算点数均小于第一类分块，所以每个小块所需寄存器不超过3d+2个，也优于CIPOHS算法。

综上所述，本发明为了充分利用部分CPU大数乘法指令提供的优良性能，提高Montgomery模乘算法的计算效率，提出了一种利用多字长大数乘法指令设计分块计算Montgomery模乘中的多精度乘法C←A*B和约减运算P←(C+M*N)/R的优化计算方法MWBCHS(Multi-Words Block Coarsely Integrated Product and Operand Hybrid Scanning)，充分利用多字长大数乘法的优势以及可用寄存器的数量，可适应不同操作数长度的多字长指令，适用于不同模数长度的Montgomery模乘，具有很强的适应性。方法以减少中间结果存储、减少计算步骤为目的设计分块计算方式，最大限度地利用多字长乘法指令单次运算长度优势，将Montgomery模乘的多精度运算和约减运算分别拆分为三类分块进行计算，相对于其他复杂的分块方法以达到优化目的而言，本专利分块方法简洁但优化效果明显。且本方法根据算法的需求重新设计扫描方式，结合算法分块类型设计适合算法的扫描方式，尽可能共享操作数以及节省存取次数，整体上减少了内存存取次数，同时减少乘法指令执行次数，提高了Montgomery模乘算法实现效率。

本发明可适应不同操作数长度的多字长指令和适用于不同模数长度的Montgomery模乘。采用分块集成的方式交替计算Montgomery多精度乘法和Montgomery约减两部分，根据所需计算的操作数实际大小和可用寄存器数量设计分块方式，充分利用多字长大数乘法指令优势减少Montgomery多精度乘法和约减的计算步骤，优化存储性能。相对于其他专利复杂的分块方法以达到优化目的而言，本方法分块方法简洁但优化效果明显。且本算法根据算法的需求重新设计扫描方式，结合算法分块类型设计适合算法的扫描方式，尽可能的共用操作数，减少内存和寄存器访问次数，减少运算指令的调用次数，提高了Montgomery模乘算法实现效率。

本发明将多字长乘法指令操作数长度和Montgomery模乘的A,B乘数作为参数，可采用不同操作数长度的多字长指令，适应不同模数长度的Montgomery模乘，可应用不同算法如ECC、RSA、国密SM2算法等计算中用到的大数模乘部分，适用性强。且重新设计Montgomery模乘算法分块方式及计算扫描顺序，优化存储性能，使用多字长大数乘法减少Montgomery多精度乘法和约减计算所需的指令实现步骤，减少系统中内存存取数据数量及存取次数(其中内存存取操作次数是指load和store操作次数的总和)，减少乘法运算指令调用次数，提高Montgomery模乘算法实现效率。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (7)

1.一种基于多字长乘法指令的Montgomery模乘快速计算方法，所述方法包括：

确定Montgomery模乘的输入参数；

将Montgomery模乘的多精度乘法运算和约减运算分别拆分为三类分块；

采用分块集成的方式，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，得到Montgomery模乘结果。

2.根据权利要求1所述的基于多字长乘法指令的Montgomery模乘快速计算方法，其特征在于，所述确定Montgomery模乘的输入参数；具体包括：

设模N为m比特，Montgomery的系数为R＝2^m，A和B为模乘运算的乘数和被乘数且A和B均为m比特，0<A,B<N，多字长乘法指令操作数T为U比特，处理器的字长为W比特，预计算常数为N’，N’为N的模R负逆，并且满足N×N’＝-1mod R，mod表示取模运算；

计算N、A和B的字数大小均为

表示取大于等于其的最小整数运算，T的字数大小为

且U≥192比特，将数A表示为：A＝(A[n-1],…,A[2],A[1],A[0])，将数B表示为：B＝(B[n-1],…,B[2],B[1],B[0])，乘积C＝A×B表示为：C＝(C[2n-1],…,C[2],C[1],C[0])，T表示为T＝(T[t-1],…,T[2],T[1],T[0])，计算参数r＝n％t，％为取余数运算符，令d＝r+t，令s为(n-d)/t的整数部分。

3.根据权利要求1所述的基于多字长乘法指令的Montgomery模乘快速计算方法，其特征在于，所述将Montgomery模乘的多精度乘法运算和约减运算分别拆分为三类分块；具体为：

将Montgomery模乘的多精度乘法运算A×B表示成第一菱形图，所述第一菱形图的每个点表示一个乘积A[i]×B[j]，0≤i,j<n；

将约减运算M×N表示成第二菱形图，所述第二菱形图的每个点表示一个乘积M[i]×N[j]，0≤i,j<n；M[i]为约减运算的参数，根据下式计算得到：

C[x]＝A[i]×B[j],x∈[0,2n-1],i∈[0,n-1],j∈[0,n-1]

M[i]＝C[x]×N’[i],x∈[0,2n-1],i∈[0,n-1]

将第一菱形图拆分成三类分块，所述三类分块每类由若干个形状相同的小块组成，第二菱形图的拆分方式与第一菱形图相同。

4.根据权利要求3所述的基于多字长乘法指令的Montgomery模乘快速计算方法，其特征在于，所述三类分块包括第一类分块、第二类分块和第三类分块；具体为：

所述第一类分块为边长为d的菱形，其中d＝r+k×t，k为正整数，根据处理器寄存器数量V选择大小，d≤V且(n-d)能被(k×t)整除，从参与模乘运算的乘数A和被乘数B的最低位开始，用于处理n/t非整除部分，所述第一类分块包括1个小块；

所述第二类分块为边长分别为d和k×t的平行四边形，位置分布为沿着第一类分块的两边，以d为底边长，分别以k×t为步长延伸到最高位，所述第二类分块包括2×(n-d)/k×t个小块；

所述第三类分块为边长为k×t的菱形，是除了第一类分块和第二类分块后余下的部分，所述第三类分块包括[(n-d)/(k×t)]²个小块。

5.根据权利要求4所述的基于多字长乘法指令的Montgomery模乘快速计算方法，其特征在于，所述采用分块集成的方式，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，得到Montgomery模乘结果；具体为：

对于第一类分块的小块计算，非整除部分使用单字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，整除部分使用多字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，得到Montgomery模乘结果；

对于第二类分块的小块和第三类分块的小块，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，同时得到这两类的Montgomery模乘结果。

6.根据权利要求5所述的基于多字长乘法指令的Montgomery模乘快速计算方法，其特征在于，所述对于第一类分块的小块计算，非整除部分使用单字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，整除部分使用多字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，得到Montgomery模乘结果；具体为：

对于非整除部分，长度为r，使用单字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，其中，

多精度乘法运算为：

{C[0],C[1]....C[r-1]}＝{C[0],C[1]...C[r-1]}+B[i]×{A[0],A[1]....A[r-1]}；

约减运算为：

{M[0],M[1]....M[r-1]}＝{C[0],C[1]....C[r-1]}×N’mod R；

P’＝{C[0],C[1]....C[r-1]}+{M[0],M[1]....M[r-1]}×N[i]，

其中i为循环控制参数，0≤i<d；

对于整除部分，使用多字长乘法指令进行Montgomery模乘的多精度乘法运算和约减运算，其中，

多精度乘法运算为：

{C[r],C[r+1]....C[r+t-1]}＝{C[r],C[r+1]...C[r+t-1]}+B[i]×{A[r],A[1+r]....A[t+r-1]}；其中，B[i]的循环遵循操作数扫描，i为循环控制参数，0≤i<d；

约减运算为：

{M[r],M[r+1]...M[r+t-1]}＝{C[r],C[r+1]....C[r+t-1]}×N’mod R，

P’＝{C[r],C[r+1]....C[r+t-1]}+N[i]×{M[r],M[r+1]...M[r+t-1]}。

7.根据权利要求5所述的基于多字长乘法指令的Montgomery模乘快速计算方法，其特征在于，所述对于第二类分块的小块和第三类分块的小块，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，同时得到这两类的Montgomery模乘结果；具体为：

对于第二类分块的小块和第三类分块的小块，使用多字长乘法指令交替计算Montgomery模乘的多精度乘法运算和约减运算，分为外循环和内循环；其中，

第二类分块的外循环控制参数为j，0≤j<s，用于计算第二类分块的相邻下一块，在j的一个外循环体内，依次进行d次内循环的多精度乘法运算，再依次进行d次内循环的约减运算；i为内循环的循环控制参数，0≤i<d，用于控制每个小块；直至外循环重复s次，得到第二类分块的Montgomery模乘结果；

第三类分块的外循环控制参数为j，0≤j<s，用于计算第三类分块的相邻下一块，在j的一个外循环体内，包括两层循环嵌套，中循环控制参数为u,0≤u<j，在u的一个循环体内，依次进行v次内循环的多精度乘法运算，再依次进行v次内循环的约减运算；v为内循环的循环控制参数，0≤v<t，用于控制每个小块；中循环重复j次后，外循环控制参数j加1，再重复中循环和内循环；直至外循环重复s次，得到第三类分块的Montgomery模乘结果。

Images