大数模幂运算与椭圆曲线标量乘运算的实现及数据仿真*

刘芳芳，刘增芳，芦殿军

（青海师范大学 数学与统计学院，青海 西宁 810008）

0 引 言

随着公钥密码体制的不断完善和发展，RSA密码体制和椭圆曲线密码体制（Elliptic Curve Cryptosystem，ECC）被越来越多的应用于密码学领域[1-3]。RSA 密码体制由 Riverst，Shamir和Adleman[4]于1978年提出，它的安全性建立在大数因子分解的困难性上，其核心运算为大数模幂运算，但大数模幂运算会耗费大量时间。2008年，陈勇涛等人[5]提出一种流水线结构的大数模运算单元，他们基于Montgomery模乘算法，得到了一种可同时进行模乘和模加减运算的大数模运算单元。2010年，李云飞等人[6]提出了一种RSA算法的改进方案，将RSA解密时的部分运算转移至加密方，使得大数模幂运算的模位数与指位数减小，进一步提升了RSA算法的性能。陈财森等人[7]于2014年提出一种针对从右到左平方乘算法实现的RSA故障分析算法，建立了故障模型，利用密钥空间搜索方法恢复密钥片断，最终恢复完整密钥。椭圆曲线密码体制分别由 Kibitz[8]和 Miller[9]独立提出，其安全性建立在椭圆曲线离散对数困难问题上，与RSA密码体制相比，ECC具有密钥长度短、计算量小、计算速度快、安全性高、存储空间小等特点，在信息安全领域具有广泛应用[10-12]。

在模素数的椭圆曲线密码体制中，标量乘运算是其核心运算。2009年，田敏[13]提出了一种计算多点标量乘的新算法，对标量从左到右进行编码，合并了编码阶段与多点标量乘的主计算阶段，有效地提高了Shamir多点乘算法的效率。徐明等人[14]于2018年优化了椭圆曲线密码系统的群运算和标量乘运算，使椭圆曲线密码系统的整体性能得到提升。同年，邬贵明等人[15]采用修正的Jacobian坐标的点加和倍点算法以及Montgomery模逆的算法，构造了一种实现素域上椭圆曲线标量乘算法的硬件结构。

本文首先基于RSA密码体制，对能快速实现大数模幂运算的平方乘算法进行分析；其次讨论平方乘算法在椭圆曲线密码体制中的对偶算法，即倍数和差算法，该算法可快速进行模素数椭圆曲线上的标量乘运算，通过利用NAF（Non-adjacent Form）算法改进标量的二进制表示方法，以提高倍数和差算法的运算效率；最后我们将对这两种算法分别进行数据仿真。

1 平方乘算法

1.1 平方乘算法的引入

在RSA密码体制中，加密和解密都存在模指数运算xcmodn。如果用一般方法直接计算xcmodn的值，即 x·x··x（modn）需要 c-1 次模乘来实现。通常情况下，指数c都非常大，这会使其运算效率低下。

考虑一个例子：计算x8，最简单的方法就是进行7次乘法运算：

而更快捷的方法只需要3次平方运算：

再看一个更一般的例子，计算x24。最简单的方法就是计算23次乘法，而更有效的方法是进行两次乘法操作（乘以x），之后再进行三次平方操作，即5次操作即可得到结果：

1.2 平方乘算法的定义

平方乘算法可以把计算xcmodn所需模乘的次数降低为最多2ℓ次，其中ℓ是c的 二进制表示的比特 数。

要计算z=xcmodn，首先要将指数c用二进制表示，且处理二进制的顺序是从高位到低位。其算法如下：

1.3 平方乘算法的实现

根据上述算法，相应的C语言程序如下所示：

#include ＜stdio.h＞

#define L 999

longSqure_and_Multiply（intx,intc,int n）｛ long z; int i,len,c_2[L]; len=0; z=1;

while（c＞0）｛ c_2[len]=c%2; len=len+1; c=c/2; ｝for（i=len-1;i＞=0;i--）

｛printf（“%d %d %d*%d”,i,c_2[i],z,z）;z=（z*z）%n; if（c_2[i]==1）｛ z=（z*x）%n;

printf（“*%d=%d ”,x,z）; ｝else printf（“=%d ”,z）; ｝return z;｝

int main（）｛intx,c,n; long data=0; printf（“input x,c,n: ”）;scanf（“%d%d%d”,&x,&c,&n）;

printf（“i ci z ”）;Squre_and_Multiply（x,c,n）;return 0; ｝

1.4 平方乘算法的数据仿真

若令n=20 001，公开加密指数c=2 929，则可通过运行以上程序，计算5 6732929mod 20 001来加密明文5 763，过程如表1所示。

表1 5 6732929mod 20 001的计算结果

所以，由此程序可以得到明文5 763所对应的密文为1 356。

2 倍数和差算法

2.1 倍数和差算法的引入

2.1.1 模素数的椭圆曲线定义

定义1[16]p＞3是素数， ℤp上的同余方程y2≡x3+ax+b（mod p）的所有解 (x,y)∈ ℤp×ℤp，连同一个特殊的点O，即无穷远点，共同构成ℤp上的椭 圆 曲 线y2=x3+ax+b。 其 中 a,b ∈ℤp是 满 足4a3+27b2≠0的常量。

2.1.2 椭圆曲线上的加法运算

定义2[16]ℤp上的椭圆曲线E的加法运算定义如下：假设 P=（x1,y1）以及 Q=（x2,y2）是 E 上的点。则可以分为以下三种情况：

情形1 若x1=x2且y2=-y1，则P+Q=O（无穷远点），即 （x,y）+（x,-y）=O，（x,y）∈ E。

情 形2 若x1=x2且y1=y2， 即P=Q， 则P+Q=（x3,y3）， 其 中 λ=（3x21+a）（2y1）-1，x3=λ2-x1-x2，y3=λ（x1-x3）-y1。 情 形 3 若 x1≠ x2， 即 P ≠ Q， 则P+Q=（x3,y3），其中 λ=（y2-y1）（x2-x1）-1，而 x3,y3的求法与情形2相同。

此外，对于所有的P∈E，定义P+O=O+P=P。在ℤp上的椭圆曲线加法运算中，（E,+）是一个阿贝尔群。

2.2 倍数和差算法的定义

设P是n阶椭圆曲线上的点。若给定c的任一带 符 号 的 二 进 制 表 示 （cℓ-1, … ,c0）， 其 中

0≤c≤n，则可以通过一系列的倍运算，以及两点之间的和差运算计算出椭圆曲线上的点P的倍数cP。其算法如下：

算 法 2 Double-And-（Add O r Subtract）（P, （ cℓ-1,… ,c0））

Q←O

算法2中的减法Q-P,应该先算出P的加法逆-P，再将所得结果与Q相加。

2.3 倍数和差算法的实现

2.3.1 模逆运算

在模素数的椭圆曲线加法运算中，求斜率λ时涉及到模逆运算m=a-1modb， a ,b ∈ℤn，因此我们首先需要利用Euclidean算法求出a,b的最大公因子，如果gcd（a,b）=1，则说明a,b互素。由于模逆运算所得值不唯一，我们取满足条件的最小值。

2.3.2 求点P的倍数

设E:y2=x3+x+6是 ℤ11上的椭 圆曲线，若令点P=（5,2），则可以通过前面所描述的倍数和差算法求得点P的倍数cP。相应的C语言如下：

#include＜stdio.h＞

#include＜math.h＞

int invert（int a1,int b1）

｛intgcd（inta,int b）;inti,s;

if（gcd（a1,b1）==1）

｛for（i=1;i＜=100000;i++）

｛if（a1＞0） ｛if（（i*b1+1）%a1==0） s=（i*b1+1）/a1;s=s%b1; ｝

else ｛if（（-1）*i*b1+1）%a1==0） s=（-1）*i*b1+1）/a1;s=s+b1; s=（b1+s%b1）%b1; ｝｝return s;｝

elseprintf（“waring ! waring ! wrong input: ”）; ｝intgcd（inta,int b）

｛intm;if（a＞0）｛m=a%b;while（m!=0） ｛a=b; b=m;m=a%b; ｝return b;｝

else｛a=a+b; m=a%b; while（m!=0） ｛ a=b; b=m;m=a%b; ｝return b; ｝｝

int temp（int x1,int y1,int x2,int y2,int *x3,int *y3）｛int invert（ int a1, int b1）; int p=11,r,m;

if（x1==0 && y1==0） ｛ *x3=x2; *y3=y2;printf（“（x3,y3）=（%d,%d） ”,*x3,*y3）; return 0;｝

if（x2==0 && y2==0） ｛ *x3=x1; *y3=y1;printf（“（x3,y3）=（%d,%d） ”,*x3,*y3）; return 0;｝

if（（x1==x2）&&（y1==-y2）） ｛*x3=0; *y3=0;printf（“（x3,y3,1）=（%d,%d） ”,*x3,*y3）;

return 0;｝

else if（（x1==x2 ）&& （y1==y2 ））｛m=invert（2*y1）%p,p）; r=（3*x1*x1+1）*m）%p;

*x3=（（r*r-x1-x2）%p）%p; *y3=（（r*（x1-*x3）-y1）%p）%p;

if（*x3＜0） *x3=*x3+p; if（*y3＜0） *y3=*y3+p;

printf（“（x3,y3,m,r,2）=（%d,%d,%d,%d） ”,*x3,*y3,m,r）;return 0; ｝

else｛m=invert（x2-x1,p）; r=（p+（y2-y1）*m）%p）%p;

*x3=（（r*r-x1-x2）%p）%p; *y3=（（r*（x1-*x3）-y1）%p）%p;

if（*x3＜0） *x3=*x3+p; if（*y3＜0） *y3=*y3+p;

printf（“（x3,y3,m,r,3）=（%d,%d,%d,%d） ”,*x3,*y3,m,r）;return 0; ｝｝

main（）｛intp=11,x1,y1,x3,y3,i,j,c[99]; printf（“please input point P:”）; scanf（“%d%d”,&x1,&y1）;

printf（“P=（%d,%d） ”,x1,y1）; printf（“input c ”）; for（j=0;j＜90;j++）｛scanf（“%d”,&c[j]）;if（c[j]==9） break; ｝x3=0; y3=0; for（i=0;i＜=j-1;i++）｛printf（“（c[%d]）=（%d） ”,i,c[i]）;

temp（x3,y3,x3,y3,&x3,&y3）; if（c[i]==1）｛temp（x3,y3,x1,y1,&x3,&y3）; ｝

else if（c[i]==-1）｛ temp（x3,y3,x1,（-y1+p）%p,&x3,&y3）; ｝else continue; ｝

printf（“（x3,y3）=（%d,%d） ”,x3,y3）;｝

2.3.3 倍数和差算法的数据仿真

令c=19，此时c的带符号的二进制表示为（1,01,0,-1），当 P=（5,2）时，通过运行以上程序，可得如表2所示。

表2 19P的计算结果

所 以， 由 此 程 序 可 得 点 P=（5,2）的 倍 数19P=（2,4）。因此完成了利用倍数和差算法求椭圆曲线E:y2=x3+x+6上的点P倍数的数据仿真。

2.4 倍数和差算法的改进

一般来讲，一个整数有多个二进制表示，例如：13=8+4+1=16-2-1。所以 （c4,c3,c2,c1,c0）=（0,1,1,0,1）或者（1,0,0,-1,-1），两者都是13的带符号的二进制表示。

在算法2中，整数c的带符号的二进制表示为（cℓ-1,… ,c0）。由于一个整数有多个二进制表示，由算法2可知，无论ci为何值时，都需要进行倍运算Q←2Q。而当ci为非零数时，则需要在倍运算的基础上再进行加法运算Q←Q+P或减法运算Q←Q-P。为了减少算法花费的时间，可以找到整数c的一种二进制表示 （cℓ-1,… ,c0），使其中的非零系数 尽可能地少。

定义3 整数c的一个带符号的二进制表示为（cℓ-1,… ,c0），如果其中的非零系数ci均是非相邻的（Non-a djacent Form，NAF），则将这样的二进制表示称为NAF表示。

下面将描述将整数c的普通的二进制表示转换成NAF表示的过程。首先这个转换的基础是将（0,1,…,1,1）转换为 （1,0,…,0,-1），由于等式 2i+2i-1+…+2j=2i+1-2j,i＞j成立，显然这种替换不会改变c的值。替换将从最右侧（即低位）的比特开始向左，按照需要重复进行。下面给出实例说明此转换过程：

显然，任何非负整数均有且只有一个NAF表示。由于在一个非负整数的NAF表示中不会出现两个相邻的非零系数，因此平均来讲，一个ℓ比特整数的二进制表示中含有个零，它的NAF表示中含有个零。由此可知NAF（c）在算法2中只需要计算次加法（或减法）运算，而普通的整数二进制表示在算法2中则需要计算次加法（或减法）运算，所以两种不同的表示在算法2中花 费时间的比率近似为

利用此技巧，可以平均提高约11%的速度。

3 结 语

针对RSA密码体制中的大数模幂运算及椭圆曲线密码体制中的标量乘运算，将大数模幂运算xcmodn对应到平方乘算法中，指数用二进制表示运算简化为计算z←z2modn及z←（z×x）modn，提高了大数模幂运算的速度；将椭圆曲线密码体制中标量乘运算cP对应到倍数和差算法中，整数c用带符号的二进制表示为运算化简为计算Q←2Q以及Q←Q±P，提高了ECC密码体制中标量乘运算的速度，且通过利用NAF表示这一技巧使倍数和差算法的速度又提高近11%。