基于250位模乘平台的Tate对最终模幂算法的改进

王晓静

(同济大学电子与信息工程学院，上海 201804)

0 引言

近年来，双线性对越来越广泛地被运用在密码学和各种加密方案中。由于它的独特性质，可以用来构造许多的加密方案和协议。例如群签名、短签名［1］、三方密钥协商［2］、可搜索加密等都要用到双线性对。

双线性对最早被应用在攻击椭圆曲线有理点群上的离散对数问题。攻击方法就是将椭圆曲线中的离散对数问题归约到有限域乘法群中的离散对数问题［3］。使用这种方法可以使得某些椭圆曲线有理点群中的离散对数问题变得容易。

双线性对的输入变量的系数可以灵活变化，它的双线性这一特点可以用来构造许多应用和协议方案。最常见的双线性对包括Weil［4］对和Tate对，以及Tate对的变种 Eta［5］对、Ate［6］对等。本文以研究 Tate 对为主，一般而言，计算Tate对的主要算法是Miller算法［7］。

在二元扩域GF(2m)中，双线性对选用超奇异曲线［7］Eb:y2+y=x3+x+b，b∈F2。对于曲线 y2+y=x3+x+1而言，目前m的推荐参数是283、353、367、457。本文Tate对计算都是在二元扩域GF(2457)中实现的。

模乘运算在硬件中实现很占资源，很多平台(如智能卡)无法实现高位的模乘运算。Jean-Claude Bajard等人在文献［8］中，介绍了一种小特征值的Montgomery模乘计算方法。本文结合文献［8］中介绍的改进的Montgomery算法和中国剩余定理，使得在只支持250 bits以内模乘运算的硬件平台上实现457 bits的二元扩域中Tate对最终模幂运算。

1 双线性对和Miller算法介绍

1.1 Tate对的基础知识

双线性对是一种如下形式的映射:

е:G1× G2→GT

其中G1和G2是加法群，GT是乘法群。

双线性对的基本性质:

(1)双线性:对任意的 P∈G1，Q∈G2，n∈Z，有 e( nP，Q)=e( P，nQ)=e(P，Q)n。

(2)非退化性:一定存在某个P∈G1和Q∈G2，使得e( P，Q)≠1，其中1代表乘法群中的单位元。

(3)可计算性:存在有效的多项式时间算法计算出双线性对的值。

设E是在域K上的椭圆曲线E/K，n∈Z，定义集合:

即是有理点集的n绕点。假设q=pm，m为整数，Fq是含有q个元素的有限域，则p是Fq的特征，且m是Fq/Fp的扩张次数，#E(Fq)表示椭圆曲线上点的个数。n是满足n|#E(Fq)的和Fq特征互素的整数，关于n的嵌入次数k满足:

传统的双线性Tate对定义［9］如下:

1.2 二元扩域Tate对的Miller算法及最终模幂运算

定义在Fq上的椭圆曲线E，满足:#E( Fq)=q+1-t(其中t是E( Fq)的迹 )，且p是有限域Fq的特征。如果t≡0 mod p，那么就认为E(Fq)是超奇异的椭圆曲线。一般而言，当Fq的特征是2时，嵌入次数最大为4。

在二元扩域GF(2m)中，双线性对常使用如下形式的超奇异曲线:

Eb:y2+y=x3+x+b，b∈F2

为了有效地计算Eb曲线上的Tate对，定义扭映射 ψ(Distortion Map)［10-11］如下:

ψ:Eb→Eb，ψ( x，y)=(x+s2，y+sx+t)

其中s2+s+1=0和t2+t+s=0。

设P、Q是椭圆曲线Eb(F2m)上的点，Tate双线性对的公式如下:

el=fl，P(ψ(Q))(qk-1)/l

其中l为椭圆曲线的阶，且ψ(Q)∈F24m是扭映射后的点，有限域F24m是F2m的扩域。

定义 fP=。Tate 对的计算公式如下［11］:

τl(P，Q)=fP(ψ(Q))22m-1

经过对文献［11］中算法优化和推导可获得如下Tate对的算法1，其推导过程略。

算法1 Tate双线性对Miller算法。

其中 P，Q∈Eb(F2m)，且 xP、yP、xQ、yQ都是基域 F2m的元素。步骤3是扩域中的运算，是最终模幂算法，由一次扩域模幂C22m和一次扩域模逆C-1完成。

扩域模幂C22m可以通过两次Frobenius模幂变换［12］算法 2 获得。

算法2 Frobenius模幂变换算法。

定义扩域元素:

C=C0+C1t+C2t2+C3t3

当m%4=3，Frobenius映射可以转化为以下等式:

当m%4=1，Frobenius映射可以转化为以下等式:

扩域模逆C-1可以通过扩域模逆算法3获得。

算法3 扩域模逆算法。

(1)利用Frobenius变换公式计算:

(2)计算扩域F24m模乘:

Ar=A23m·A22m·A2m·A

(3)计算基域模逆:

Inv=(Ar)-1

(4)计算扩域F24m模乘:

A-1=Inv·(Ar-1)

步骤1通过重复调用Frobenius模幂变换算法完成。步骤2通过3次扩域全模乘(算法4)完成。步骤3是基域F2m的模逆运算，可以使用二进制扩展Euclidean算法［13］来完成。步骤4是基域元素和扩域元素相乘，可调用算法4。其中Ar-1=A23m·A22m·A2m，步骤2已经计算，可以避免重复运算。

算法4 扩域全模乘算法［11］。

输入参数:

不可约多项式:t4+t+1

输出参数:

算法3步骤2和步骤4调用算法4，且在计算457 bits的 c0、c1、c2、c3时，使用改进的 Montgomery 算法5和CRT算法6。其中的模乘运算都转换为250 bits以内的模乘运算，调用硬件完成运算。

2 改进的Montgomery模乘算法和中国剩余定理

算法5 改进的Montgomery模乘［8］。

输入参数:

两个多项式A，B ∈Fp［X］，且 deg (A)，d eg (B)≤k-1;多项式 N ∈Fp［X］，且 deg(N )=k是首一的不可约多项式(阶k的系数是1);多项式阶deg(Ψ)，deg(Ψ')≥k，且满足 gcd(Ψ，Ψ')=gcd(Ψ，N)=1，即最大公约数为1。deg运算值是多项式的阶。

输出参数:

MM(A，B)=ABΨ-1mod N

运算步骤:

和一般的Montgomery模乘比较，引入了Ψ'。改进算法的思想是将约减多项式Ψ分解为两个不可约多项式，多项式A模Ψ1和Ψ2。这种方法将多项式A通过模运算分解为两个多项式表达式。其表达式如下:

假设计算A×B mod N，其中A、B、N都是二元扩域GF(2457)中的元素，也可以认为是457项的多项式。

算法6 中国剩余定理(CRT算法)［14］。CRT方程:

其中ni间互素，记N=n1n2…nk;Ni=N/ni;ai≡N-1imod ni;使用二进制扩展 Euclidean算法［13］实现模逆。可得:

将457 bits模乘运算转化为250 bits以内的模乘。所以只需要两项CRT方程，即利用250 bits的多项式运算转换为457 bits模乘运算，如下所示:

x≡x1(mod n1)

x≡x2(mod n2)

n1和n2是250 bits以内的不可约多项式。

3 Miller算法中最终模幂的改进分析

在智能卡等嵌入式平台上，由于资源的限制及模乘运算的复杂性，硬件可能无法支持高位的模乘运算。在支持250 bits的模乘运算平台上，采用改进方案来进行二进制域GF(2m)中457 bits的Tate对最终模幂运算。本文中的模乘运算都调用250 bits以内硬件模乘实现。

Tate对最终模幂算法流程如图1所示。

图1 最终模幂的计算流程图

如图1所示，457 bits的Tate对最终模幂C(22m-1)由一次扩域模幂部分和一次扩域模逆部分组成。通过Frobenius变换算法5进行简化，扩域模幂部分直接调用Frobenius变换得到。457 bits扩域模逆部分由457 bits的全模乘算法计算得到(457 bits全模乘由457 bits的模乘组成)。通过改进的Montgomery算法和CRT算法可以将457 bits的模乘转化为250 bits内的模乘运算，调用硬件来完成250 bits内的模乘运算。计算参数和步骤如下(本文中后续计算以以下的计算参数为准，也可以采用其他满足下列条件的参数)。

计算参数:

在 GF(2457)中，C=c3t3+c2t2+c1t+c0;c3，c2，c1，c0是基域 F2457中的元素。

定义扩域中元素:C=(x451)t3+x。

算法5中的Ψ分解为不可约多项式Ψ1和Ψ2，即Ψ=Ψ1×Ψ2。

不可约三项式:

Ψ'分解为不可约多项式Ψ'1和Ψ'2，即Ψ'=Ψ'1×Ψ'2。

不可约三项式:

GF(2457)的不可约多项式:

N=x457+x61+1

通过上述定义，预计算出 M=Ψ2mod N、N-1mod Ψ1等包含 N，Ψ，Ψ'，Ψ1，Ψ2，Ψ'1，Ψ'2的常量。

3.1 计算最终模幂扩域模幂部分

最终模幂的(C22m)的运算(m=457)使用两次Frobenius变换公式计算可得:

3.2 计算最终模幂扩域模逆部分

最终模幂的扩域模逆部分C-1的运算使用扩域模逆算法3计算，如下所示:

调用全模乘算法4计算Cr:

同理计算出调用全模乘算法Cr-1=C23m·C22m·C2m。

算法3的步骤3可以使用二进制扩展Euclidean算法［13］来完成基域 F2m模逆 Inv=(Cr)-1。

再次调用全模乘算法 4，完成 C-1=Inv·( Cr-1)。

3.3 计算最终模幂

最后再次调用全模乘算法4，获得C(22m-1)=C22m·C-1=c'0+c'1t+c'2t2+c'3t3，c'0，c'1，c'2，c'3的结果过长，省略。

3.4 改进的Montgomery模乘方案Sage仿真验证

Sage仿真中使用上文中的计算参数进行验证，将457 bits的模乘计算使用250 bits的模乘完成。同时，定义HardMul为250 bits内的硬件模乘，使用程序对改进算法进行测试。Sage自带CRT的计算功能。

如上述程序所示，经过计算后得到:

R1←(A×Β×phi-1)mod N

再次调用上述程序，预先计算好457 bits的M(phi2mod N)作为新的输入参数A'，前面计算的R1作为新的输入参数B'，计算得R2:

通过Sage验证，可以证明通过上述的方法能够使用250 bits的模乘运算实现457 bits的模乘，且上述结果和实际结果一致。通过Sage的仿真，亦可证明改进算法的正确性。

3.5 改进的Montgomery模乘运算计算步骤

全模乘算法4中的模乘运算都可以采用如下步骤完成，将457 bits的模乘转化为250 bits以内的硬件模乘，所有参数都参照上述的计算参数。

目标:A×B mod N

输入参数:A，B

固定参数:N，Ψ，Ψ'，Ψ1，Ψ2，Ψ'1，Ψ'2

例子:计算C23m·C22m调用全模乘算法4。其中:

改进方案计算步骤如下:

按照上述步骤，可以计算出R= (A×B×Ψ-1)mod N。

对应用而言，参数 N，Ψ，Ψ'，Ψ1，Ψ2，Ψ'1，Ψ'2都是固定的。所以可以预先计算好N-1mod Ψ1，Ψ2mod N等不变的参数减少开销，等需要时再载入相关参数。

调入预先计算好457 bits的M(Ψ2mod N)作为新的输入参数A'，前面计算的R作为新的输入参数B'。按照上述改进方案计算步骤继续运算得R':

按照上述步骤计算得:

如上所述，例子中计算扩域模乘C23m·C22m的系数C0和实际结果一致。重复调用上述计算步骤计算出其他的系数(C1，C2，C3)，和预期值一致。计算C0中的R3和前面的Sage仿真结果一致。由于是在二元扩域中的运算，所以在上述计算步骤3中计算Q的负号可以省略。

最终模幂算法在扩域模逆部分需要使用改进的Montgomery模乘方案。在扩域模逆算法3的步骤2和步骤4需要使用到全模乘算法。全模乘算法中的457 bits模乘通过改进的算法替换为250 bits的硬件模乘运算。

4 结束语

本文的改进方法可以使得只支持250 bits硬件模乘的平台完成457 bits的模乘，最终支持457 bits二进制域双线性对Tate对最终模幂的运算。本质上，457 bits的最终模幂的计算就是Frobenius变换和457 bits的全模乘运算。457 bits全模乘运算由457 bits的模乘组成。通过改进的 Montgomery算法和CRT算法将457 bits的模乘运算转化为250 bits以内的模乘运算，从而使得只支持250 bits硬件模乘的平台实现457 bits的Tate对最终模幂。由于多项式参数 N，Ψ，Ψ'，Ψ1，Ψ2，Ψ'1，Ψ'2是固定的，可以预先计算出Ψ2mod N等预计算参数数据，等使用时再调用该参数数据，减少额外的开销。改进方案也可以通过将多项式Ψ分解为3个或者更多的不可约多项式来支持更高的Tate对最终模幂运算。在效率方面，由于运算步骤的增加会降低一定的运算效率。

［1］Miller V S.Short Programs for Functions on Curves［DB/OL］.http://crypto.stanford.edu/miller/miller.pdf，1986-05-06.

［2］Joux A.A one round protocol for tripartite Diffie-Hellman［C］//Proceedings of the 4th International Symposium on Algorithmic Number Theory.2000:385-394.

［3］赵昌安，张方国.双线性对有效计算研究进展［J］.软件学报，2009，20(11):3001-3009.

［4］Boneh D，Franklin M K.Identity-based encryption from the Weil pairing［C］//Proceedings of the 21st Annual International Cryptology Conference on Advances in Cryptology.2001:213-229.

［5］Hess F，Smart N P，Vercauteren F.The Eta pairing revisited［J］.IEEE Transactions on Information Theory，2006，52(10):4595-4602.

［6］Zhao C，Zhang F，Huang J.A note on the Ate pairing［J］.International Journal of Information Security，2008，7(6):379-382.

［7］Alfred Menezes，Scott Vanstone，Tatsuaki Okamoto.Reducing elliptic curve logarithms to logarithms in a finite field［C］//Proceedings of the 23rd Annual ACM Symposi-um on Theory of Computing.1991:80-89.

［8］Jean-Claude Bajard，Laurent Imbert，Graham A Jullien，et al.A CRT-based Montgomery Multiplication for Finite Fields of Small Characteristic［DB/OL］.http://hal.archives-ouvertes.fr/docs/00/10/64/55/pdf/d512.pdf，2005-07-15.

［9］Wang Mao-cai，Hu Han-ping，Dai Guang-ming.An efficient algorithms for Tate pairing computation［J］.Journal of Communication and Computer，2007，4(8):20-23.

［10］Alfred Menezes.An Introduction to Pairing-based Cryptography［DB/OL］.http://cacr.uwaterloo.ca/～ ajmeneze/publications/pairings.pdf，2008-12-20.

［11］Soonhak Kwon.Efficient Tate pairing computation for supersingular elliptic curves over binary fields［C］//Proceedings of the 10th Australasian Conference on Information Security and Privacy.2005:134-145.

［12］McCusker K，O’Connor N，Diamond D.Low-energy finite field arithmetic primitives for implementing security in wireless sensor networks［C］//Proceedings of the 2006 International Conference on Communications，Circuit and Systems.2006，3:1537-1541.

［13］Darrel Hankerson，Alfred Menezes，Scott Vanstone.椭圆曲线密码学导论［M］.张焕国译.北京:电子工业出版社，2005:55.

［14］Henri Cohen，Gerhard Frey，Roberto Avanzi，et al.Handbook of Elliptic and Hyperelliptic Curve Cryptography［M］.Chapman＆ Hall/CRC，2006:196.