模幂滑动窗口法分析及加法链在预计算中的应用

屈 晓1，孙达志1,2

(1.天津大学计算机科学与技术学院，天津 300072；2.中国科学院信息工程研究所，北京 100093)

1 概述

1978年RSA公钥密码算法被提出以来，公钥加密体制逐步成为了信息安全领域应用最广泛的技术。模幂运算作为公钥体制的核心环节，其执行速度决定着公钥体制的总体效率。然而计算机执行大数模幂运算速度慢这一问题，至今没有得到很好的解决。

大数模幂运算最基本的方法是依次计算底数的各次模幂值，然而这种方法效率低下。目前已有许多关于大数模幂运算问题的研究，提出了二进制法、滑动窗口法、重编码法等多种方法，使大数模幂运算速度慢这一问题得到了一定缓解。

本文首先对滑动窗口法的复杂度进行分析，利用马尔可夫链求出任意位任意窗口长度下的非零窗口数的一般表达式，然后提出一种应用于预计算部分的改进算法。

2 基本方法

文献[1]总结了二进制法。二进制法将指数e表示成二进制串的形式：

其中，k为指数e的位数。由左至右扫描每一位，每次先做一次模平方计算，若扫描位为1，再做一次模乘计算。二进制法的平均复杂度为3(k-1)/2。

分块法[1]作为二进制法的改进方法，采用空间换时间的思想来加速模幂运算过程。分块法预计算并存储指数长度为d的所有可能的幂值，将指数e的二进制串划分为长度为d的块，根据块值进行计算。分块法的平均复杂度为：

3 滑动窗口法及复杂度分析

3.1 滑动窗口法

文献[2]中的滑动窗口法是一种应用非常广泛的模幂方法。滑动窗口法类似分块法，将指数e划分为零窗口和非零窗口，只需一半预计算量。这里将滑动窗口法重新阐述为算法1。

算法1

输入M,e,n

输出C

(1)计算并存储Mwmodn，w=3,5,…,2d–1。

(2)根据以下策略由左至右划分指数e，得到窗口F1,F2,…,Fs：

1)初始状态S为零窗口状态，用ZW表示，非零窗口状态用NW表示。

2)当S=ZW时，扫描一位，若是0，归入ZW，S=ZW；若是1，归入新的NW，S=NW。

3)当S=NW时，扫描d–1位，将本NW结束在最后一位非0位设为i位，i+1至d–1位归入ZW，S=ZW。

(3)初始C=MF1。

(4)i从2到s做循环：

1)C=CEimod n ，这里Ei=2Li，Li为窗口Fi的长度。

2)若Fi为NW，C=CMFimod n。

(5)返回C。

3.2 复杂度分析

文献[1]提到一种带参数q,r的变长滑动窗口法，文献[3]对其进行了复杂度分析，但存在着逻辑错误。文献[4]给出了其近似分析及一个精确分析的例子，并未给出精确复杂度的一般表达式。下面利用马尔可夫状态转移矩阵对算法1进行精确的复杂度分析。

滑动窗口法复杂度可以表示为：

其中，T表示总复杂度；PRE表示预计算次数；LW表示最左窗口长度；NW表示非零窗口数量。

(1)PRE：预计算需要计算除去M1的指数不大于2d的所有指数为奇数的模幂值，加上M2，所以PRE=2d–1。

(2)LW：最左窗口最高位为1，窗口长度取决于d位中最后一个非零位，因此：

(3)NW：对于每一个扫描位定义3类状态：

1)ZO：未进行回溯即归入ZW的0。

2)Ni：归入NW 的第i位，i=1,2,…,d。

3)Zi：被NW扫描的第i位，后回溯归入ZW，i=2,3,…,d。定义状态转移矩阵的行列状态分别依次为：

状态转移矩阵P如下：

其中：

(1)ZO：已确定当前位为ZW，下一窗口状态只由下一位决定，下一位的状态只可能是ZO,N1，概率均为2–1。

(2)N1～Nd–1：当前位为i，只有当i+1～d 位全部为0，下一位转入ZW，以概率2–(d–i)转为ZO，否则以1–2–(d–i)的概率停在NW，转为Ni+1。

(3)Nd：当前位已为该窗口最后一位，下一窗口状态只由下一位决定，下一位的状态只可能是ZO,N1，概率均为2–1。

(4)Z2～Zd–1：当前位为i，由于此位经扫描后回溯归入ZW，因此i+1～d位必然全部为0，i位以1的概率从Zi转为Zi+1。

(5)Zd：当前位已为上一非零窗口扫描过的最后一位，下一窗口状态只由下一位决定，下一位的状态只可能是ZO,N1，概率均为2–1。

由于指数e最高位为1，初始状态的概率分布为π0=[0,1,0,…,0]。

于是第i位状态概率分布为πi=π0Pi–1。

3.3 实验数据对比

为了验证上述理论的正确性，选取较短的128位、使用及分析普遍的512位、1024位和较长的4096位，窗口d选取3～6，每个单元选取10000个随机数进行实际划分，获取实际复杂度。表1中每单元第1行为实际复杂度，第2行为理论公式计算得出的复杂度。

表1 总复杂度对比

复杂度测量值说明：由于模乘所需时间远大于其他计算所需时间，因此只记模乘的次数，一次模平方记为一次模乘。总复杂度为PRE+n–LW+NW，详见3.2节。

4 预计算中的加法链应用

4.1 加法链

若序列S=(a0,a1,…,as)满足以下性质，则称S为e的一条加法链：

(1)S单调递增。

(2)a0=1,as=e。

(3)ai=aj+ak,0≤j≤k

利用最短加法链加速模幂运算是理想的方法之一，但求给定整数的最短加法链本身是一个NPC问题[5]。文献[6]介绍了几种求最短加法链的方法，但从实验可以看出，在随机选取的十进制4位整数中，求最短加法链耗时已经将近7 s。而二进制法、分块法、滑动窗口法都可以看成是可行的求指数加法链的方法。

4.2 算法分析及步骤

在滑动窗口法确定窗口大小d后，对Mwmodn,w=3,5,…,2d–1，进行预计算，加上计算M2，共2d–1次计算。

然而，对于任意给定位数的指数确定其最优窗口大小，通常需要采集大规模样本，即传统的暴力搜索法；或在一定范围内逐一测试。两者都相当耗费资源。文献[7]给出了一个对分块法窗口大小的最优化估计。

另一方面，当窗口大小d选择过大时，预计算量呈指数型增长，预计算结果很可能没有被充分利用，造成资源浪费。因此，使预计算结果得到充分的利用将是一件有意义的事，而且可以提供其他提高效率的方向，例如4.5节的部分。

利用加法链进行预计算的意思是，通过所有划分出来的需要进行预计算的非零窗口值寻找一条加法链。前文提到过寻找最短加法链是一个NPC问题，而寻找一条通过多个给定整数的最短加法链比寻找一个固定整数的加法链问题更加复杂。针对这一问题本文给出一种计算机执行简单可行的方法，来求出一条通过所有需要进行预计算值的加法链。算法2给出了算法的细节。

算法2

输入 经算法1第(2)步划分后得到的非零窗口值F1,F2,…,Fs

输出 加法链A=a0,a1,…,as

(1)对F1,F2,…,Fs按升序排列并且只保留不同值得到序列W0=w01,w02,…,w0i。

(2)保存a0=w0i，计算t0=w0i–w0i–1。

(3)若t0已在W0中出现或t0=1，则在W0中删去w0i，得到W1=w01,w02,…,w0i–1。否则在W0中用t0替换w0i，得到W1=w01,w02,…,w0i–1,t0。

(4)类似的，重复步骤(1)～步骤(3)，直至Wi中只剩下一个元素wi1。

(5)用任何一种计算可行的方法求出wi1的一条加法链，方便起见这里用二进制法，得到ai,ai+1,…,as。

(6)对A=a0,a1,…,as按升序排列。

(7)返回A。

算法2求出的A即是一条包含所有非零窗口值的加法链。之后根据A依次计算出Mai，i=0,1,…,s，即可得到所有预计算所需的模幂值。

例 选择窗口长度d=6，经过步骤(1)后，得到W0=7,19,21,35,47,55,61；

经过步骤(2)后，得到a0=61,t0=6；

经过步骤(3)后，用6替换61，并按升序排列得到W1=6,7,19,21,35,47,55；

第1次循环后，有：

第2次循环后，有：

第3次循环后，有：

…

第13次循环后，有：

最后，2的一条加法链为1,2，添加到A后按升序排列得：

根据加法链A计算所有非零窗口值共需14次计算，而传统方法则需要计算出全部值，共26–1=32次计算。

4.3 实验数据

同3.3节，此处同样选取128位、512位、1024位、4096位，窗口d选取4～7，每个单元选取10000个随机数进行实际划分，利用算法2进行预计算，结果如表2所示。表中数据为相应位数下计算划分出来所有非零窗口值所需模乘次数。传统方法的预计算量在d为4～7时，分别为固定的8、16、32、64。

表2 预计算量对比

4.4 上下界分析

在4.3节的实验中，可以看出，应用算法2虽然在d选择过大时，对预计算量增长影响较小，比传统方法全部计算呈指数型增长具有较大的改进，但在选择最优d的情况下，并没有减少较多的模乘次数，可以说效率基本相同。对于一个固定整数的最短加法链问题，其长度的上界为：[lbe]+H(e)–1，文献[8]给出了其下界：lbe+lbH(e)–2.13。文献[9]指出了加法序列问题，即寻找一条通过多个给定整数的最短加法链问题，与向量加法链间存在一一对应的关系。文献[10]给出了一个向量加法链的界：l(r1,r2,…,rt)=lbr+(t+o(1))lbr/lblbr，其中，r=max{r1,r2,…,rt}。

当窗口长度d取得越大，划分出来的窗口数就越少。应用加法链思想简化预计算，当d取得足够大及加法链长度足够短时，就可以在总复杂度上低于传统方法的总复杂度。文献[11]提出了一种方法，可以在一定程度上缩短加法链的长度，但其方法十分复杂，计算机执行起来存在一定困难，非常消耗资源。

4.5 其他应用

在信息发送中，会出现这样一种情况：发送方A有一条消息M，需要发送给多个接收方B,C,D等。这一问题相当于需要计算Cb=Mb,Cc=Mc,Cd=Md,…。通常情况下，会对这些模幂单独进行计算，如果应用加法链的思想，即求一条包含b,c,d,…的加法链，然后根据加法链进行模幂计算，将会大大降低需要执行的模幂运算次数。然而一个完整的指数不像窗口，对多个指数应用加法链的思想，将需要不少的空间，这种方法还有许多待改进的地方。

5 结束语

本文给出了有限位指数模幂计算应用滑动窗口法的平均复杂度精确表达式，非零窗口中的各位及零窗口中的零位都可表示成马尔可夫链中的某一状态，各状态间以固定概率相互转化，并通过实验进行误差对比，在选取的各个条件下，误差均小于0.1次模乘。提出一种计算机执行简单可行的求加法序列的算法，使之能应用于预计算部分，实验表明，随着选择窗口长度的增加，该方法比传统方法提升的效率愈加明显，使预计算的结果得到100%利用。今后将对算法2进行改进，在计算机执行简单可行的前提下，进一步缩短加法链的长度；对于同一消息多方发送问题，考虑如何只利用少量空间应用加法链提高加密效率。文献[12]指出二进制法、分块法等方法在本质上都是加法链算法，今后将研究基于幂树思想的构造加法链的算法，结合滑动窗口法减少非零窗口数，从而减少总复杂度。

[1]Knuth D E.The Art of Computer Programming:Semi-numerical Algorithms[M].New Jersey,USA:Addison-Wesley,1981.

[2]Menezes A J,Van Oorschot P C,Vanstone S A.Handbook of Applied Cryptography[M].[S.l.]:CRC Press,1996.

[3]Koc C K.Analysis of Sliding Window Techniques for Exponentiation[J].Computers&Mathematics with Applications,1995,30(10):17-24.

[4]Park H,Park K,Cho Y.Analysis of the Variable Length Nonzero Window Method for Exponentiation[J].Computers&Mathematics with Applications,1999,37(7):21-29.

[5]Downey P,Leong B,Sethi R.Computing Sequences with Addition Chains[J].SIAM Journal on Computing,1981,10(3):638-646.

[6]王晓东.最短加法链算法[J].小型微型计算机系统,2001,22(10):1250-1253.

[7]曾 皓,范明钰,王光卫,等.模幂与点乘m_ary算法中窗口大小的最优化估计[J].计算机应用研究,2007,24(10):35-36.

[8]Schönhage A.A Lower Bound for the Length of Addition Chains[J].Theoretical Computer Science,1975,1(1):1-12.

[9]OlivosJ.On Vectorial Addition Chains[J].Journal of Algorithms,1981,2(1):13-21.

[10]Yao A C C.On the Evaluation of Powers[J].SIAM Journal on Computing,1976,5(1):100-103.

[11]Bos J,Coster M.Addition Chain Heuristics[C]//Proc.of the 9th Annual International Cryptology Conference.New York,USA:Springer,1989:400-407.

[12]董付国,厉玉蓉.几种方幂模快速算法的加法链一致性分析[J].计算机工程与应用,2010,46(36):48-50.