CN1607518A

CN1607518A - 在Montgomery乘法内利用SIMD指令

Info

Publication number: CN1607518A
Application number: CN200410085541.5A
Authority: CN
Inventors: P·L·蒙特格米里
Original assignee: Microsoft Corp
Current assignee: Microsoft Technology Licensing LLC
Priority date: 2003-10-15
Filing date: 2004-10-15
Publication date: 2005-04-20
Anticipated expiration: 2024-10-15
Also published as: EP1524594A3; DE602004017559D1; US20050084099A1; JP2005122141A; EP1524594B1; EP1524594A2; KR101103893B1; CN100437548C; ATE413642T1; US7532720B2; JP4662744B2; KR20050036698A

Abstract

描述了一种体系结构和方法论，用于在支持SIMD指令的计算机***上执行Montgomery乘法。

Description

在Montgomery乘法内利用SIMD指令

技术领域

本发明涉及支持SIMD(单指令，多数据)指令的计算机体系结构。

背景

一些计算机体系结构支持SIMD(单指令，多数据)指令的执行。每个这样的指令对多个数据集合执行一项运算。在执行期间，使用一项运算的多个同时发生的实例来处理不同的数据。举例来讲，一个指令可能会执行几项独立的减法或其他数学运算(但到处是相同的运算)。

来自英特尔公司的奔腾_4微处理器是能够执行SIMD指令的计算机体系结构的一个例子。该微处理器利用被称作“SSE2(流型SIMD扩展2)指令”的144个指令来作用于被存储在寄存器中的数据。该微处理器具有八个128位整数寄存器，其中的每个整数寄存器可以保存多个数据集合(例如，两个64位整数、四个32位整数、八个16位整数或十六个8位整数)。每个单独的SSE2指令对寄存器中所包含的这多个数据集合执行一项运算——可能利用来自另一个寄存器的输入。这些SSE2指令减少被用来执行特定的程序任务的指令总数，从而提高总性能。

如今，要求微处理器在几乎一瞬间执行许多乘法和除法运算。安全性是在要求计算机执行许多精确运算的情况下的一个例子。安全函数使用复杂的密码算法，这些密码算法经常使用涉及许多模数乘法和运算的求幂运算(例如，RSA算法)。寻找减少复杂算法内的计算，从而改善性能，这是不断努力的目标。

Montgomery乘法是众所周知的可避免除法的模数乘法的算法。可通过从右边(而不是从左边)减少双倍字长乘积，来实现这一点。

概述

描述了一种体系结构和方法论，用于在支持SIMD指令的计算机***上执行Montgomery乘法。在所描述的实施中，使用SSE2(流型SIMD扩展2)指令来执行Montgomery乘法所利用的机器级运算。

内容简述

图1示出支持SIMD指令的计算机体系结构。

图2示出SIMD指令在Montgomery乘法(执行模数乘法)内的运用，以及数据值在寄存器中的布置。

图3是用于在图1中的计算机体系结构上执行Montgomery乘法的过程的流程图。

本说明书中，同一标号表示相同元件和特征。

详细说明

Montgomery乘法是一种广泛使用的、免于除法的模数乘法算法。下文描述了双向SIMD(单指令，多数据)指令在Montgomery乘法内的运用。可以在支持SIMD指令的各种计算机体系结构上执行Montgomery乘法，同时，将在支持作用于两个独立的操作数的SIMD指令的计算机体系结构的上下文中描述一项特定的实施。

计算机体系结构

图1表现了计算机体系结构100，该计算机体系结构可以被配置成为Montgomery乘法而使用双向SIMD指令。计算机体系结构100包括微处理器102和存储器104。微处理器102具有一个或多个ALUs(算术逻辑部件)106(1)、…、106(N)，以管理数学运算(例如，加法和乘法)、逻辑算子(比如“或”、“与”、“异或”等)以及右移和左移。微处理器102进一步包括一个或多个寄存器108(1)、…、108(M)，以保存ALUs 106所产生的中间值和最终结果。

微处理器102的一个例子是英特尔_IA-32微处理器，它具有两个ALUs 106和八个128位整数寄存器108。寄存器108也被称作“XMM寄存器”。来自英特尔公司和微软公司的编译器支持用户对这些XMM寄存器进行访问，而不要求该用户书写汇编码。

计算机体系结构100支持SIMD指令。在我们的例子中，英特尔_IA-32微处理器支持SSE2(流型SIMD扩展2)指令。SSE2指令集主要增强音频和视频压缩过程，从而促使几项增强致力于推进MPEG 2编码和文件加密过程。SSE2指令集包括在其他事物之中支持128位SIMD整数算术运算的144个指令。这些SSE2指令也意在协助对程序进行矢量化(vectorizing)。例如，如果一个作用于32位数据，那么，每个128位XMM寄存器可以保存四个32位值。通过使用SSE2指令，可以利用单指令来处理这四个32位操作数，并获得四个32位结果。

使用一连串SIMD指令110来执行Montgomery乘法——它是关于模数乘法的算法。这些指令一般被示作驻留在微处理器102中，但它们可能在ALUs106(1)-106(N)中的任何ALUs内加以执行。也许，在密码函数内的模数求幂期间，这些指令指示该处理器执行Montgomery乘法中所涉及的计算。

模数乘法和Montgomery乘法

为了介绍模数乘法和Montgomery乘法，可考虑普通的除法问题。让n₁和n₂是整数，并且n₂＞0。通过将n₁除以n₂，欧几里得除法算法可提供唯一的整数商q和余数r，以便：

n₁＝qn₂+r 0≤r≤n₂

出于讨论的目的，数学表达式“quo”和“rem”被定义如下：

quo(n₁，n₂)＝q(商)；以及

rem(n₁，n₂)＝r(余数)。

模数乘法的数学定义以正整数模数N开始，这里假设该模数是固定的。给出两个整数a和b，则其模数乘积modmul(a，b)被定义为rem(ab，N)。模数乘法的一项可能的实施是：首先计算乘积ab，然后将这个乘积除以N，以产生商q和余数r，以便：

(1a)ab＝qN+r

这种算法输出余数r＝rem(ab，N)＝modmul(a，b)。

本申请中被提名的发明者Peter Montgomery(彼得·Montgomery)观察到整数除法在许多计算机体系结构上速度缓慢，他发现了另一种关于模数乘法的算法(即Montgomery乘法)，这另一种算法避免了这些除法，并且在利用单一模数执行外延计算时速度更快。关于模数乘法的Montgomery算法方面的背景信息，可指引读者注意以下的出版物：Peter L.Montgomery的《无试除的模数乘法》(计算数学，第44卷，第170号，1985年4月，第519-521页)，Donald E.Knuth的《计算机编程技术》(第2卷，“半数字算法”，第3版本，Springer-Verlag，1998年，练习4.3.1-41)。

Montgomery乘法要求与模数N互质的固定整数R，从而意味着：固定整数R和模数N的最大公分母是一(即，gcd(N，R)＝1)。为了方便继续讨论，假设：模数N是奇数，整数R是2的幂(除了在使用R＝10的例子内以外)。提出这个假设，可存在整数N’，NN’≡1(mod R)。N’的值只取决于N(和R)，并且每一模数N只需要被计算一次。

给出整数a和b(并利用上下文已知的N、R)，Montgomery乘积montmul(a，b)可被定义为整数r’，以便0≤r’＜N以及

ab≡r′R

这个整数r’存在(并且是唯一的)，这是因为R是可逆模数N。关于Montgomery乘法的公式如下所示：

(M1)montmul(a，b)＝rem((ab-qN)/R，N)，这里q＝rem(abN’，R)

为了验证方程式(M1)，可定义r₁＝(ab-qN)/R，其中，r₁被认为是整数。那个断定从以下公式中得出：

qN≡(abN′)N＝ab(N′N)≡ab(mod R)

接下来，确认r’＝rem(r₁，N)满足方程式(1b)的条件，因为：

r′R＝rem(r₁，N)R≡r₁R＝ab-qN≡ab(mod N)

通过“rem”的定义，整数r’也满足0≤r’＜N的条件。

如果R是2的幂，并且，如果R大得足以促使|ab-qN|＜NR，则Montgomery的方程式(M1)的右侧容易在二进制计算机上评估。如果0≤a，b＜N＜R，则确保后者。

为了看到Montgomery乘积“montmul”的有效性，可为整数x定义函数“tomont”，以便tomont(x)＝rem(xR，N)。然后，

montmul(tomont(a)，tomont(b))R≡tomont(a)tomont(b)//montmul的定义

≡(aR)(bR)＝(abR)R≡tomont(ab)R(mod N)//tomont的定义

假设gcd(N，R)＝1允许用R除这个叠合，余下：

montmul(tomont(a)，tomont(b))≡tomont(ab)(mod N)

这两边必须相等，因为两者都在区间[0，N-1]内，并且它们是同余模N。已知ab≡modmul(a，b)，并且，如果x≡y(mod N)，则tomont(x)＝tomont(y)，这样可作出结论：

montmul(tomont(a)，tomont(b))＝tomont(ab)＝tomont(modmul(a，b))

相应地，如果使用modmul来定义域上的乘法并使用montmul来定义图像上的乘法，则函数“tomont”是从Z/NZ到Z/NZ的乘法同态(意味着：乘积的图像是图像的乘积)。这里的Z/NZ表示整数模数N。进一步的调查揭示了：

tomont(a)+tomont(b)≡tomont(a+b)(mod N)

tomont(a)-tomont(b)≡tomont(a-b)(mod N)

所以，函数“tomont”是环同态(具有是rem(R，N)的1的图像)。从以下公式中，

如果并且只有当a≡b(mod N)时，tomont(a)＝tomont(b)

函数“tomont”是从Z/NZ到Z/NZ的环同构，而不只是同态。

多倍精度Montgomery乘法

已介绍了Montgomery乘法，该章节展示多倍精度数字在Montgomery乘法中的运用。出于讨论的目的，考虑模数乘法在密码函数内的运用。例如，RSA密码***要求模数求幂——其中，该模数可能会是512个位、1024个位或更长。大多数计算机体系结构只支持32×32→64位或64×64→128位乘积。当乘以更长的数字时，每个操作数被分成几个较小的(通常是32位或64位)部分。

让r＝2³²或2⁶⁴，这取决于计算机文字尺寸。找到正整数k，以便r^k＞N(通常选择服从该限制的尽可能小的k)。设置R＝r^k。

本文中使用大写字母(例如，A)来表示以下的非负多倍精度数字。相当于N的值可以由k个基数-r位数来表示。利用括号中的指数所跟随的相同的字母来书写这些位数。如果0≤j≤k-1，那么，A[j]表示基数r中的A的第j个位数(第0个位数最不重要)。也就是说，

A＝∑_0≤j≤k-1A[j]r^j 0≤A[j]≤r-1

书写这个的另一种方法是A＝(A[k-1]A[k-2]…A[0])_r，其中，最后的下标表示基数r。

给出这些符号，一个关于montmul(A，B)的算法(其中，0≤A，B＜N＜R)是：

T₁：＝0；

T₂：＝0；

DISCARD：＝0；

for j from 0 to k-1 do

mul₁：＝A[j]；

T₁：＝T₁+mul₁*B； //Temporarily need k+1 digits

tlow：＝T₁[0]； //tlow：＝rem(T₁，r)

T₁：＝(T₁-tlow)/r； //T₁：＝quo(T₁，r)

DISCARD：＝DISCARD+tlow*r^j； //For presentation purposes only

mul₂：＝rem((tlow-T₂[0])*N′[0]，r)； //So new T₂ will be an integer

T₂：＝(T₂+mul₂*N-tlow)/r；

end for；

return rem(T₁-T₂，N)；

通过归纳核对，可发现在每个迭代的末尾处：

0≤T₁＜B

0≤T₂＜N

0≤DISCARD＜r^j+1

T₁r^j+1+DISCARD＝(A[j]A[j-1]...A[0])_r*B

T₂r^j+1+DISCARD≡0(mod N).

两个推论是：

-N＜T₁-T₂＜B

(T₁-T₂)r^j+1≡(A[j]A[j-1]...A[0])_r*B(mod N).

在最后的迭代(j＝k-1)完成之后，这些结果是：

-N＜T₁-T₂＜B＜N

(T₁-T₂)R≡AB(mod N).

所以，函数montmul(A，B)将会是T₁-T₂(当T₁≥T₂时)或T₁-T₂+N(当T₁＜T₂时)。

利用较小的修改，可以使对变量T₁和T₂的这些更新显得十分类似。利用tlow和mul₂的初始化中的小变化，主循环成为：

for j from 0 to k-1 do

mul₁：＝A[j]；

tlow：＝rem(T₁[0]+mul₁*B[0]，r)；

mul₂：＝rem((tlow-T₂[0])*N′[0]，r)；

T₁：＝(T₁+mul₁*B-tlow)/r；

T₂：＝(T₂+mul₂*N-tlow)/r；

end for；

在此重写之后，T₁更新增加B的倍数，而T₂更新增加N的倍数。T₁、T₂、B和N数组表示多倍精度数字≤N＜R＝r^k，其中，当在基数r中书写这些数字时，全部都具有长度k。乘数mul₁和mul₂是单精度(即，0≤mul₁，mul₂＜r)。更新T₁所需要的这些指令十分类似于更新T₂所需要的指令。这样，用于更新T₁和T₂的该指令序列(考虑B和N是多倍精度)如下所述：

[更新代码1]

sum₁：＝T₁[0]+mul₁*B[0]； //0 ≤sum₁≤r²-r

//Use rem(sum₁，r)＝tlow to compute

mul₂(formula suppressed here)

sum₂：＝T₂[0]+mul₂*N[0]；

sum₁：＝quo(sum₁，r)； //0≤sum₁≤r-1

sum₂：＝quo(sum₂，r)；

for i from 1 tok-1 do

sum₁：＝sum₁+Tm[i]+mul₁*B[i]； //0≤sum₁≤r²-1

sum₂：＝sum₂+T₂[i]+mul₂*N[i]；

T₁[i-1]：＝rem(sum₁，r)；

T₂[i-1]：＝rem(sum₂，r)；

sum₁：＝quo(sum₁，r)； //0≤sum₁≤r-1

sum₂：＝quo(sum₂，r)；

end for；

T₁[k-1]：＝sum₁；

T₂[k-1]：＝sum₂；

当代的计算机体系结构支持流水线技术和指令级平行，其中，与更新数组T₁有关的运算后面是与更新数组T₂有关的运算。为了进一步改善性能，需要通过每一对更新使用一个指令，来在i上的循环内一起更新数组T₁和T₂。

关于Montgomery乘法的SIMD指令

图1中所展示的计算机体系结构100使用SIMD指令110(以及在我们特定的例子中是SSE2指令)来同时更新Montgomery乘法montmul(A，B)中所使用的数组T₁和T₂。这个过程部分涉及：按从右到左的方式来处理A的位数；并且，交错涉及B(即，mul₁*B[i])和模数N(即，mul₂*N[i])的乘法的结果，B和模数N被用来更新数组T₁和T₂。以下将更加详细地描述这一点。

如上所述，SSE2指令有助于程序的矢量化，其中，可以将多个数据集合保存在单一寄存器中。为了展示这个要点，关于涉及32位数据的操作，每个128位寄存器108可以保存四个32位值。然后，循环可以被书写如下：

int32x[100]，y[100]，z[100]；//假设int32数据类型占据32个位

for i from 0 to 99 do

z[i]：＝x[i]+y[i]；

end for；

这个循环可以被(程序设计员或编译器)局部展开，所以，生成的代码类似：

int32 x[100]，y[100]，z[100]；

for i4 from 0 to 96 by 4 do

z[i4]：＝x[i4]+y[i4]；

z[i4+1]：＝x[i4+1]+y[i4+1]；

z[i4+2]：＝x[i4+2]+y[i4+2]；

z[i4+3]：＝x[i4+3]+y[i4+3]；

end for；

一个SSE2指令将来自邻接的存储位置的四个32位数据值(x[i4+3]，x[i4+2]，x[i4+1]，x[i4])载入XMM寄存器。另一个SSE2指令将(y[i4+3]，y[i4+2]，y[i4+1]，y[i4])载入不同的XMM寄存器。单一的PADDD指令形成所有四个32位总和，如下所示：

(x[i4+3]+y[i4+3]，x[i4+2]+y[i4+2]，x[i4+1]+y[i4+1]，x[i4]+y[i4])。

然后，将这四个总和置入z数组的邻接的位置。通过每次处理四个数组元素，来按四的因数减少循环迭代的数目(从100到25)。另一个效果是：编译器可以保留四个32位寄存器，以保存循环指数i4和数组地址x、y、z，而不是弄乱具有数据值(例如，x[i4])的那些寄存器。

在Montgomery乘法的情况中，使用SSE2指令来一起更新数组T₁和T₂。而在最后的例子中，同时处理四个32位整数值；这里，使用SSE2指令来处理两个64位值。在继续的例子中，假设：基数r＝2³²。将参照图2来描述该实施，它展示了范例寄存器和SSE2命令。

观察以上所述的伪代码[更新代码1]，sum₁和sum₂的值始终在区间[0，r²-1]＝[0，2⁶⁴-1]内。使用一个XMM寄存器108(1)来保存这两个64位值(sum₁，sum₂)。其他的XMM寄存器108(2)、108(4)和108(5)(或存储位置)将保存以下各对64位值：

(mul₁，mul₂)

(B[i]，N[i]) (0≤i≤k-1)

(T₁[i]，T₂[i]) (0≤i≤k-1)

虽然mul₁、mul₂、B[i]、N[i]、T₁[i]、T₂[i]的数值都适合32个位，但是，为每个数值分配64个位。也使用常数(r-1，r-1)，它在十六进制中是00000000FFFFFFFF00000000FFFFFFFF。这个常数被保存在寄存器108(3)中。

一旦选择这些数据格式，伪代码[更新代码1]的内循环就容易转变。如所述，八个XMM寄存器中的三个XMM寄存器108(1)、108(2)和108(3)保存各个对(sum₁，sum₂)、(mul₁，mul₂)和(r-1，r-1)。在开始内循环之前，对所有三个寄存器进行初始化。如图2中的两个装入指令INST1：LOAD和INST2：LOAD所表现的，内循环主体将各个对(B[i]，N[i])和(T₁[i]，T₂[i])载入XMM寄存器108(4)和108(5)。

一个SSE2指令“PMULUDQ”(它乘以两个32位值，以产生64位无符号的乘积(即，32×32→64位无符号的乘积))执行两个运算mul₁*B[i]和mul₂*N[i](即图2中的INST3：PMULUDQ)。然后，使用两个SSE2指令“PADDQ”(它增加两个64位值)来产生：

sum₁：＝sum₁+T₁[i]+mul₁*B[i]；

sum₂：＝sum₂+T₂[i]+mul₂*N[i]；

将被更新的(sum₁，sum₂)放置在另一个XMM寄存器108(6)中(即图2中的INST4：PADDQ和INST5：PADDQ)。复制操作(即INST6：MOVDQA)制作(sum₁，sum₂)的复制品。然后，一个SSE2指令“PAND”(即128位“与”，它等同于具有(r-1，r-1)的两个64位“与”)(即INST7：PAND)和存储操作(即INST8：STORE)执行以下函数：

T₁[i-1]：＝rem(sum₁，r)

T₂[i-1]：＝rem(sum₂，r)

PAND毁坏(sum₁，sum₂)的一个副本，但保留另一个副本。然后，按照32个位的SSE2运算“PSRLQ”(即，右移具有填零的64位数据)提供以下结果(即INST9：PSRLQ)：

sum₁：＝quo(sum₁，r)

sum₂：＝quo(sum₂，r)

排斥关于地址计算和循环控制的指令(使用32位寄存器，而不是XMM寄存器)，内循环主体只需要两次装载、一次乘法、两次加法、一次拷贝、一次按位“与”、一次存储、一次移位——九个SSE2指令，来更新数组T₁的一个元素和数组T₂的对应的元素。微软Visual Studio.NET和英特尔^_IA-32编译器具有允许用户在不使用汇编语言的条件下生成这种类型的代码的固有本质。

注意，按从i＝0到i＝k-1的顺序来处理B[i]，从而虑及A中的各个位数的从右到左的处理，A是数组B的乘数。A上的循环的每个迭代再处理全部的B。此外，注意，当i增进时，乘法mul₁*B[i]可以与乘法mul₂*N[i]交错，从而导致这些计算的重排。然而，该过程有效率地产生所需的最终结果，以下将更加详细地展示这一点。

存储布局

以上所示的伪代码[更新代码1]假设B[i]和N[i]在邻近的64位文字中，因为SSE2取指令(MOVDQA)要求正在被装载的数据在存储器中是邻接的。这个存储器也应该在16字节边界上对准。但当调用montmul(A，B)时，通过B[k-1]的输入B[0]通常将会在邻接的32位文字的一个数组中(在4字节边界上对准)，通过N[k-1]的N[0]将在另一个这样的数组中。那不是所需的存储模式。相应地，可以实施计算机来执行额外的循环初始化，以创建所需的数组布局。

通过将B数组和N数组交替***具有128位元素的临时数组(并对其进行适当的对准)，可避免该问题。可以早期在montmul中制作这些副本。当k很大时，这些副本的成本是可以忽略的，这是因为以后将访问被拷贝的每一对元素达k次。

同样，数组T₁和T₂的这些元素可以在具有128位元素的数组内交错。最后的rem(T₁-T₂，N)计算将按标准形式(即邻近的32位文字)来存储其输出。

运算

图3表现了用于在Montgomery乘法内使用SIMD指令的过程300。例如，可以由计算机体系结构100来执行过程300，并且，将参照图1和图2来描述过程300。过程300被展示为一系列方框，其中的每个方框表示由计算机体系结构执行的运算或泛函性。这些运算或函数可以在软件、固件和/或硬件中加以执行。照此，可以使用这些方框来表示被存储在介质上的指令，这些指令在被执行时可实行所述泛函性。

在方框302处，对在计算函数montmul(A，B)的过程中所使用的值进行初始化。这些值包括：整数输入A、B；临时数组T₁、T₂；以及模数N。在方框304处，临时数组T₁、T₂被设置为零。将关于输入B和模数N的值交替***矢量，并且，来自这些数组的一对元素将适合这些寄存器之一(例如，图2中的寄存器108(4))。

在方框306处，迭代循环开始从右到左地处理输入A的位数。在方框308处，假设将利用输入B的倍数来更新数组T₁(乘数mul₁是来自数组A的位数)，则处理器确定：模数N的什么倍数mul₂允许更新T₁、T₂以相同的位数结束。在方框310处，执行(A的位数)乘以B的乘法(即，伪代码“更新代码1”中的mul₁*B[i])，以及乘以模数N的确定倍数的乘法(即，伪代码“更新代码1”中的mul₂*N[i])。在方框312处，更新数组T₁、T₂。通过使用两个数组T₁和T₂，乘法AB与Montgomery乘法中的乘法qN交错。如决定方框314所表现的，循环继续进行，直到输入A中的所有位数都已被处理为止。在方框316处，返回结果T₁-T₂(mod N)。

Montmul计算内的循环重排的展示

现在将利用现实、但虚构的数据来描述Montgomery函数montmul的例子。首先，导出montmul计算的所需答案，随后讨论如何使用以上所描述的技术来获得这个结果。

设A＝123，B＝456，该问题是计算montmul(123，456)——R＝1000(十进制)，N＝789。通过定义，输出将会是整数r’，以便0≤r’＜789，并且

123*456≡1000r′(mod 789)

解NN’≡1(mod R)，发现整数N’是109(即，789 N’≡1(mod 1000)，得出N’＝109)。如果以前使用过模数N＝789，则可以一次解答关系789 N’≡1(mod 1000)，以获得N’＝109。为了检验该结果，将N’＝109***关系789*109＝86001，其mod R＝1000结束于001中。注意，如果这是模数789的第一次使用，则可以使用被扩展的欧几里得算法。

接下来，通过使用以上导出的Montgomery乘法函数[M1](即，montmul(A，B)＝rem(AB-qN)/R，N)，其中的q＝rem(AB N’，R))，来计算乘积AB，以产生56,088(即，123*456＝56,088)。***前导零(即056088)，以形成六位数，以便保存两个三位数数字的乘积。在模数R＝1000处，结果中的最后三个位数是088。为了发现商q，将结果088乘以整数N’＝109，以产生9592(即，088*109＝9592)。再次，在模数R＝1000处，这个结果的末端的三位数是592(即，q＝592)。将商q＝592乘以模数N＝789，并从123*456的乘积中减去(即，AB-qN)，可提供：

056088-592*789＝056088-467088＝-411000

通过构建，这个结果结束于000中。将该值-411000除以R(即，形成“rem”函数(AB-qN)/R中的第一个元素)，最终结果是-411。将三个结果插回到Montgomery函数montmul中可得出最后的所需答案：

montmul(123，456)rem(-411，789)＝378

参考所需答案，本章节的剩余部分通过使用SIMD指令，在以上所描述的多倍精度实施中运用这些范例值。本说明中的多倍精度版本(r＝10，k＝3)的目的是：在余数的底部三次***一个零位数，而不是同时***所有三个零。它部分通过从右到左地处理输入A中的位数(始于个位数，终于百位数)，来完成这一点。

由于N’＝109，其低位位数是N’[0]＝9。乘积AB或056088的低位位数是8。所以，乘积056088*N’结束于与8*9＝72相同的位数中。最后的位数是2，所以，从乘积AB(即056088)中两次减去模数N(即789)：

056088-2*789＝056088-1578＝054510

下一个步骤考虑所产生的值054510中的十位数1。乘积05451*N’(在抑制尾随零之后)结束于与1*9＝9相同的位数中。这样，决定从结果054510中九次减去模数N(即789)——具有一个尾随零占位符，以确保十位数(即7890)的处理；从而可得出：

054510-9*7890＝054510-71010＝-016500

为该结果中的百位数5再一次重复这个步骤，乘积-0165*N’(抑制这两个尾随零)结束于与5*9＝45相同的位数中，它具有最后的位数5。从结果-016500中五次减去模数N(即789)——具有两个尾随零占位符(即78900)，这样可得出：

-016500-5*78900＝-016500-394500＝-411000

早先，输出是montmul(123，456)＝rem(-411，789)＝378。这里，我们从056088中减去2*789、9*7890和5*78900，从而给出与我们早先从056088中减去592*789时相同的结果。

通过使用两个临时数组T₁和T₂，乘法AB(即，123*456＝056088)可以与缩减AB-qN(即，056088-592*789)交错，以获得结果-411000。将为数组T₁加入456的倍数，用于计算乘法AB；同时，将为数组T₂加入789的倍数，用于计算qN的缩减。T₁和T₂将仍然是非负的。对T₁＝T₂＝0进行初始化，乘法AB或123*456在123的低位数或3处开始：

T₁＝T₁+3*456＝0+1368＝1368

当前的T₂＝0结束于零中。为了使数组T₂象数组T₁＝1368的当前最后位数一样结束于8中，可确定将模数N(即789)乘以2的因数(即，2*9＝18，它结束于8中，它在被加入0时等于8)：

T₂＝T₂+2*789＝0+1578＝1578

放弃两个数组T₁和T₂中的尾随8的位数(即，关于循环不变量，DISCARD＝8)，从而余下：

T₁＝136 T₂＝157

A＝123中的下一个位数(从右边运作)是2。

T₁＝T₁+2*456＝136+912＝1048

当前的T₂＝157结束于7中。为了使数组T₂象数组T₁＝1048的当前最后位数一样结束于8中，可确定将模数N(即789)乘以9的因数(即，9*9＝81，它结束于1中，它在被加入7时等于8)：

T₂＝T₂₁+9*789＝157+7101＝7258

放弃两个数组T₁和T₂中的尾随8的位数(即，关于循环不变量，DISCARD＝88)，从而余下：

T₁＝104 T₂＝725

A＝123中的最后一个位数(从右边运作)是1：

T₁＝T₁+1*456＝104+456＝560

当前的数组T₂＝725中的最后的位数是5。为了使数组T₂象数组T₁＝560的当前最后位数一样结束于0中，可确定将模数N(即789)乘以5的因数(即，5*9＝45，它结束于5中，它在被加入5时等于10，它结束于零中)：

T₂＝T₂+5*789＝725+3945＝4670

放弃两个数组T₁和T₂中的尾随0的位数(即，关于循环不变量，DISCARD＝088)，从而余下：

T₁＝56 T₂＝467

观察到：T₁＝56具有123*456＝56088的前导位数，而T₂＝467具有592*789＝467088的前导位数。共享的末端位数在DISCARD＝088中。输出是：

montmul(123，456)＝rem(56-467，789)＝rem(-411，789)＝378

注意，虽然这些计算已被重新安排，但是，该过程产生相同的最终结果。按不同的顺序来加入六个被加数：

003*456＝001368

020*456＝009120

100*456＝045600

-002*789＝-001578

-090*789＝-071010

-500*789＝-394500

-411000

利用以上所描述的计算机体系结构100，并通过使用SIMD指令，可一起更新数组T₁和T₂。例如，可以同时执行以下两个公式：

T₁＝136+2*456＝1048

T₂＝157+9*789＝7258

首先观察T₁更新的末端(即单元)数位：6+2*6＝18。这结束于8中，而157结束于7中，所以，T₂更新中的789的乘数将会是(8-7)N’[0]≡9(mod 10)。一旦确定这个乘数，就发现T₂更新的末端位数是7+9*9＝88。放弃18和88中的(共享的)尾随8，从而将进位矢量初始化成高位数(1，8)。

进行T₁和T₂更新，但使用SIMD指令来一起实行它们，下一个步骤作用于这些输入的十位数(下标1)：

临时变量＝进位+(T₁[1]，T₂[1])+(2，9)*(B[1]，N[1])

＝(1，8)+(3，5)+(2，9)*(5，8)

＝(14，85)

临时变量的这两个元素是从0到99的整数。末端位数(4，5)被保存为(T₁[0]，T₂[0])的新的值，上位数被放置在进位矢量中(即，进位＝(1，8))，这偶然与进位矢量的老的值相同。

对百位数同样这样实行。

临时变量＝进位+(T₁[2]，T₂[2])+(2，9)*(B[2]，N[2])

＝(1，8)+(1，1)+(2，9)*(4，7)

＝(10，72)

末端位数(0，2)被保存为(T₁[1]，T₂[1])的新的值，上位数被放置在进位矢量中(即，进位＝(1，7))。

我们已用完输入A＝123中的所有位数。进位矢量(即，进位＝(1，7))的内容被存储为(T₁[2]，T₂[2])的新的值。这些新的值(其最后的8已经被放弃)是T₁＝104和T₂＝725。

到临时变量计算的输入是从0到9的整数，输出的范围从0到9。关于硬件寄存器108，这对应于0至2³²-1的输入以及0至2⁶⁴-1的输出。当从存储器中参考(T₁[2]，T₂[2])时，这些操作数位于邻近的位置，因为这是它们如何被存储在先前的迭代上。更明确地说，即使它们的值适合32个位，它们也被存储为邻近的64位位置。象(B[2]，N[2])之类的各个对早先被拷贝到邻近的位置，并且，当处理A＝123的每个位数时，再使用这些对。

范例源代码

以下提供了用于执行伪代码[更新代码1]的范例源代码。

一些宏在这个代码内被加以参考，但没有在这里加以定义。ptr_roundup_m128i(临时变量)采用在4字节边界上对准的临时数组地址，并且返回在16字节边界上对准的临时变量、临时变量+4、临时变量+8、临时变量+12中的任何一个。给出128位值x，SSE2_digit2(x)宏返回x的位63-32中的该值，并且，SSE2_digit0(x)宏返回x的位31-0中的该值。

digit_t、digit_tc、DWORDREG、DWORDREGC都是32位类型。尾随“c”或“C”识别恒定的操作数。

一些变量对应是：

文本 modmul_from_right_SSE2

A，B b，a

N 模数

N’[0] minv

k lng

(B[i]，N[i]) a_modulus[i]

(T₁[i]，T₂[i]) temp12[i]

(2³²-1，2³²-1) mask02

(mul₁，mul₂) mul12

(sum₁，sum₂) prod12或carry12

mul₂——将要被加入T₂的N的倍数的计算将通常类似：

mul₂＝(T₁[0]+mul₁*B[0]-T₂[0])*N′(mod r)

这个公式化要求：在可开始乘以N’的乘法之前，结束乘法mul₁*B[0]。而该源代码是使用：

mul₂＝(T₁[0]-T₂[0])*N′+mul₁*(B[0]*N′)(mod r)

其B[0]*N’(mod r)被预先计算。如果硬件安排这类计算，则修订的mul₂公式允许这两个整数乘法重叠。

This is the first of several modmulchx files with algorithms for

    modular multiplication，both FROM_LEFT and FROM_RIGHT.

    All procedures have an argument list

        (digit_tc *a，*b，               //Two numbers to multiply

                                         //0＜＝a，b＜modulus

        digit_t  *c                      //Product.May overlap a or b.

        mp_modulus_tc*pmodulo，          //struct with information about modulus

        digit_t  *temps)                 //Temporaries(length

                                         //pmodulo-＞modmul_algorithm_temps)

        FROM_LEFT codes return

            c＝＝(a*b)mod modulus

        FROM_RIGHT codes return the Montgomery product

            c＝＝a*b/RADIX^lng(mod pmodulo-＞modulus).

    where lng＝pmodulo-＞length.

        This file starts with

            modmul_from_left_default

            modmul_from_right_default

    which work on all architectures and lengths.It also has

            modmul_from_right_SSE2

    which is specific to X86 hosts with SSE2 instructions(e.g.，Pentium 4).
*/
/******************************************************************************/
static BOOL WINAPI modmul_from_left_default

            (digit_tc  *a，               //IN

            digit_tc  *b，                //IN

            digit_t   *c，                //OUT

            mp_modulus_tc *pmodulo，      //IN

            digit_t  *temps)              //TEMPORARIES，at least 2*lng
/*

        This implements ordinary modular multiplication.

    Given inputs a，b with 0＜＝a，b＜modulus＜RADIX^lng，

    and where lng＞0，we form

        c＝＝(a*b)mod modulus.
*/
{

    BOOL OK＝TRUE；

    DWORDREGC lng＝pmodulo-＞length；

    digit_t*temp1＝temps；      //Length 2*lng

    assert(pmodulo-＞modmul_algorithm_temps＝＝2*lng)；

    OK＝OK &amp;&amp; multiply(a，lng，b，lng，temp1)；    //Double-length product

    OK＝0K &amp;&amp; divide(temp1，2*lng，pmodulo-＞modulus，lng，

                   &amp;pmodulo-＞left_reciprocal_1，digit_NULL，c)；

    return OK；
}   //modmul_from_left_default
				
				<dp n="d16"/>
/******************************************************************************/
static BOOL WINAPI modmul_from_right_default

       (digit_tc  *a，             //IN

        digit_tc  *b，             //IN

        digit_t  *c，              //OUT

        mp_modulus_tc*pmodulo，    //IN

        digit_t  *temps)           //TEMPORARIES，at least 2*lng
/*

        This implements Montgomery(FROM_RIGHT)multiplication.

    Let lng＝pmodulo-＞length＞0.

    Given inputs a，b with 0＜＝a，b＜pmodulo-＞modulus＜RADIX^lng，we form

           c＝＝a*b/RADIX^lng(mod pmodulo-＞modulus)。

        At the start of the loop on i，there exists templow

    (formed by the discarded values of

    LOW_DIGIT(prod1)＝LOW_DIGIT(prod2))such that

           0＜＝temp1，temp2＜modulus

           temp1*RADIX^j+templow＝b[0：j-1]  *a

           temp2*RADIX^j+templow＝＝0(mod modulus)

    when j＝lng，we exit with c＝＝temp1-temp2(mod modulus)
*/
(

    BOOL OK＝TRUE；

    DWORDREGC lng＝pmodulo-＞length；

    digit_t *temp1＝temps，*temp2＝temps+lng；  //Both length lng

    digit_tc *modulus＝pmodulo-＞modulus；

    digit_tc minv＝pmodulo-＞right_reciprocal_1；

    digit_tc minva0＝minv*a[0]； //mod RADIX

    DWORDREG i，j；

    digit_t carry1，carry2，mul1，mul2；

    dblint_t prod1，prod2；

    assert(pmodulo-＞modmul_algorithm_temps＝＝2*lng)；
//  Case j＝0 of main loop，with temp1＝temp2＝0 beforehand。

    mul1＝b[0]；

    mul2＝minva0*mul1；    //Modulo RADIX

    carry1＝HPRODUU(mul1，a[0])；          //Top of product

    carry2＝HPRODUU (mul2，modulus[0])；

    asert(mul1*a[0]＝＝mul2*modulus[0])；  //mod RADIX

    for(i＝1；i！＝lng；i++) {

    prod1＝MULTIPLY_ADD1(mul1，      a[i]，carry1)；

    prod2＝MULTIPLY_ADD1(mul2，modulus[i]，carry2)；

    temp1[i-1]＝LOW_DIGIT(prod1)；

    temp2[i-1]＝LOW_DIGIT(prod2)；

    carry1＝HIGH_DIGIT(prod1)；

    carry2＝HIGH_DIGIT(prod2)；
}
temp1[lng-1]＝carry1；
temp2[lng-1]＝carry2；
for(j＝1；j ！＝lng；j++) (

    mul1＝b[j]；

    mul2＝minva0*mul1+minv*(temp1[0]-temp2[0])；  //Modulo RADIX

    prod1＝MULTIPLY_ADD1(mul1，a[0]，temp1[0])；

    prod2＝MULTIPLY_ADD1(mul2，modulus[0]，temp2[0])；

    //Replace temp1 by(temp1+b[j]*a-LOW_DIGIT(prod1))/RADIX

    //Replace temp2 by(temp2+mul2*modulus-LOW_DIGIT(prod2))/RADIX

    assert(LOW_DIGIT(prod1)＝＝LOW_DIGIT(prod2))；

    carry1＝HIGH_DIGIT(prod1)；

    carry2＝HIGH_DIGIT(prod2)；
				
				<dp n="d17"/>
　　   for(i＝1；i!＝lng；i++){

　　      prod1＝MULTIPLY_ADD2(mul1，      a[i]，temp1[i]，carry1)；

　　      prod2＝MULTIPLY_ADD2(mul2，modulus[i]，temp2[i]，carry2)；

　　      temp1[i-1]＝LOW_DIGIT(prod1)；

　     　 temp2[i-1]＝LOW_DIGIT(prod2)；

　　      carry1＝HIGH_DIGIT(prod1)；

　　      carry2＝HIGH_DIGIT(prod2)；

　　   }     

　　   tamp1[lng-1]＝carry1；

　　   temp2[lng-1]＝carry2；

　　}

　　OK＝OK&amp;&amp;sub_mod(temp1，temp2，c，modulus，lng)；

　　return OK；
}//modmul_from_right_default；
/*********************************************************/
#if TARGET＝＝TARGET_IX86&amp;&amp;USESSE2
static BOOL WINAPI modmul_from_right_SSE2

　　  (digit_tc*a，    //IN

　　  digit_tc*b，     //IN

　　  digit_t*c，      //OUT

　　  mp_modulus_tc*pmodulo，    //IN

　　  digit_t   *temps)    //TEMPORARIES，at least 8*lng+3
/*

　　  With modmul_from_right_default(above)，the computations on
   (a，carry1，mul1，prod1，temp1)are very similar to those on
   (modulus，carry2，mul2，prod2，temp2).
   This module takes advantage thereof，using X86 SSE2 instructions.

　　  The last code often fetches a[i] and modulus[i] together.
   We construct one temporary array with lng 128-bit words (hence 4*lng digit_t 
   elements)，whose i-th element contains (0，a[i]，0，modulus[i]).
   Another array of this length has (0，tamp1[i]，0，tamp2[i]).
*/
{
   BOOL OK＝TRUE；
   DWORDREGC lng＝pmodulo-＞length；
   digit_tc*modulus＝pmodulo-＞modulus；
   digit_tc minv＝pmodulo-＞right_reciprocal_1；
   digit_tc minva0＝minv*a[0]；//mod RADIX
   const_m128i mask02＝_mm_set_epi32(0，RADIXM1，0，RADIXM1)；
   _m128i*temp12＝ptr_roundup_m128i(temps)；//Length lng
   _m128i*a_modulus＝temp12+lng；           //Length lng
   _m128i carry12，mul12，prod12；
   DWORDREG i，j；
   assert(pmodulo-＞modmul_algorithm_temps＝＝8*lng+3)；
   for (i＝0；i!＝lng；i++)(

　　  a_modulus[i]＝_mm_get_epi32(0，a[i]，0，modulus[i])；

　　  tamp12[i]＝_mm_setzero_si128()；
   )
   for(j＝0；j!＝lng；j++)(

　　  digit_tc mul1＝b[j]；

　　  digit_tc mul2＝minva0*mul1

　　      +minv*(SSE2_digit2(temp12[0]-SSE2_digit0(temp12[0]))；

　　                                            //Modulo RADIX

　　  mul12＝_mm_set_epi32(0，mul1，0，mul2)；

　　  prod12＝_mm_add_epi64(temp12[0]，_mm_mul_epu32(mul12，a_modulus[0]))；

　　  //Replace temp1 by  (temp1+b[j]*a-LOW_DIGIT(prod1))/RADIX

　　  //Replace temp2 by  (temp2+mul2*modulus-LOW_DIGIT(prod2))/RADIX

　　  assert (SSE2_digit0(prod12)＝＝SSE2_digit2(prod12))；

　　  carry12＝_mm_srli_epi64(prod12，32)；

　　         //Upper halves of products：(0，carry1，0，carry2)
				
				<dp n="d18"/>
    for(i＝1；i！＝lng；i++)  {

        prod12＝_mm_add_epi64(_mm_add_epi64(carry12，temp12[i])，

                               _mm_mul_epu32(mul12，a_modulus[i]))；

        temp12[i-1]＝_mm_and_sil28(mask02，prod12)；  //Lower halves

        carry12＝_mm_srli_epi64(prod12，32)；         //Upper halves

    }

    temp12[lng-1]＝carry12；
   }//for j
//   Do the equivalent of sub_same，but with non-contiguous input.
   {

     digit_t borrow＝0；

     for (i＝0；i ！＝lng；i++)  (

        digit_tc ai＝SSE2_digit2 (temp12[i])；

        digit_tc bi＝SSE2_digit0 (temp12[i])；

        digit_tc ci＝ai-bi-borrow；

        c[i]＝ci；

        borrow＝ai^((ai^bi)|(ai^ci))；  //MAJORITY(～ai，bi，ci)

        borrow＞＞＝(RADIX_BITS-1)；

    }

    if(borrow！＝0)borrow-＝add_same(c，modulus，c，lng)；

    assert(borrow＝＝0)；
  }
  return OK；
}  //modmul_from_right_SSE2
#endif//TARGET_IX86

结论

虽然已用针对结构特点和/或方法论动作的语言描述了本发明，但是，将会理解：所附权利要求书中所定义的本发明不一定局限于所描述的特殊的特点或动作。相反，这些特殊的特点和动作被揭示为执行所声明的发明的示范形式。

Claims

1.一种计算机***，其特征在于，它包括：

存储器；以及

支持SIMD指令的处理器，该处理器被配置成：使用SIMD指令来执行Montgomery乘法。

2.如权利要求1所述的计算机***，其特征在于，该处理器正在执行密码函数，并且使用Montgomery乘法来计算该密码函数中的幂。

3.如权利要求1所述的计算机***，其特征在于，该处理器保持两个数组，以保存来自Montgomery乘法的中间计算；并且使用这些SIMD指令来同时更新这两个数组。

4.如权利要求1所述的计算机***，其特征在于，所述Montgomery乘法涉及输入数组的第一个乘法和模数数组的第二个乘法，并且使用这些SIMD指令来同时执行这第一个和第二个乘法。

5.如权利要求1所述的计算机***，其特征在于，所述Montgomery乘法具有指令循环，并且该循环的每个迭代涉及：除复制操作外，使用不多于八条SIMD指令。

6.如权利要求5所述的计算机***，其特征在于，这些SIMD指令包括两条加载指令、一条乘法指令、两条加法指令、一条复制指令、一条按位“与”指令、一条存储指令和一条移位指令。

7.一种处理***，其特征在于，包括：

具有一组寄存器的处理器，该处理器被配置成支持SIMD指令；以及

一组SIMD指令，它们可由该处理器执行，以执行所述Montgomery乘法：

montmul(A，B)＝rem((AB-qN)/R，N)，其中，q＝rem(AB N’，R)

式中，A和B是整数，q是商，N是模数，R是与模数N互质的整数，N”是整数，使得NN’≡1(mod R)。

8.如权利要求7所述的处理***，其特征在于，这些SIMD指令包括同时执行乘法AB和qN的各个部分的单SIMD指令。

9.如权利要求7所述的处理***，其特征在于，整数B和模数N作为数组来加以执行，并且，使用至少一条SIMD指令，以便利用用于计算AB的B的倍数来更新第一数组T₁，并利用用于计算qN的N的倍数来更新第二数组T₂。

10.如权利要求9所述的处理***，其特征在于，使用单SIMD指令来同时更新第一数组T₁和第二数组T₂。

11.如权利要求9所述的***，其特征在于：

第一寄存器保存B和N数组的元素；

第二寄存器保存第一数组T₁的元素和第二数组T₂的元素；以及

第三寄存器用来保存第一数组T₁和第二数组T₂的结果，第一数组T₁利用B的倍数来加以更新，第二数组T₂利用N的倍数来加以更新。

12.一种计算机可读介质，其特征在于，包括计算机可执行SIMD指令，当被执行时，这些计算机可执行SIMD指令指示处理器执行Montgomery乘法。

13.如权利要求12所述的计算机可读介质，其特征在于，包括：

第一SIMD指令，用于将数组B和N的元素载入第一寄存器；

第二SIMD指令，用于将数组T₁和T₂的元素载入第二寄存器；

第三SIMD指令，用于将数组B中的元素乘以第一个倍数，并将数组N中的元素乘以第二个倍数；

第四和第五SIMD指令，用于将第三SIMD指令的结果加入由第二SIMD指令装载的这些数组元素上，并加入从先前迭代中被加以保存的任何进位；

第六和第七SIMD指令，用于将第五SIMD指令的每个输出分成两个大小缩减的结果，一个适合数组T₁和T₂，另一个表示关于下一个迭代的进位；以及

第八SIMD指令，用于在存储器中更新数组T₁的元素和数组T₂的元素。

14.如权利要求12所述的计算机可读介质，其特征在于，这些SIMD指令包括SSE2指令。

15.一种用于计算Montgomery乘法

montmul(A，B)＝rem((AB-qN)/R，N)，其中，q＝rem(AB N’，R)

的方法，其中，A和B是整数，q是商，N是模数，R是与模数N互质的整数，N”是整数，以便NN’≡1(mod R)，其特征在于，该方法包括：

从右到左地为整数A的每个位数迭代地执行：

利用由输入B乘以整数A的该位数的乘积来更新的数组T₁，确定模数N的什么倍数允许被更新的数组T₁、T₂以相同的位数结束；

将输入B乘以整数A的该位数，并将模数N乘以这个被确定的倍数；以及

更新数组T₁、T₂。

16.如权利要求15所述的方法，其特征在于，该执行包括：使用SIMD指令。

17.如权利要求15所述的方法，其特征在于，由单SIMD指令来执行该乘法。

18.如权利要求15所述的方法，其特征在于，进一步包括：在所述执行之前，对数组T₁、T₂和模数N进行初始化。

19.一个或多个计算机可读介质，其特征在于，存储计算机可执行指令，当被执行时，这些计算机可执行指令执行如权利要求15所述的方法。

20.一种方法，其特征在于，包括：

利用在执行Montgomery乘法的过程中所使用的值，来对一组寄存器进行初始化；以及

利用关于存储在这些寄存器中的这些值的SIMD指令，来计算Montgomery乘法。

21.如权利要求20所述的方法，其特征在于，该计算包括：使用Montgomery乘法来计算密码函数中的幂。

22.如权利要求20所述的方法，其特征在于，该计算包括：使用SSE2指令。

23.如权利要求20所述的方法，其特征在于，Montgomery乘法具有指令循环，并且使用不多于九条SIMD指令来执行该循环的每个迭代。

24.如权利要求23所述的方法，其特征在于，这九条SIMD指令包括两条加载指令、一条乘法指令、两条加法指令、一条复制指令、一条按位“与”指令、一条存储指令和一条移位指令。