此处所谓求逆运算,是指在模乘群里求逆。
第一节里提到互质的两个定义:
(1)p,q两整数互质指p,q的最大公约数为1。
(2)p.q两整数互质指存在整数a,b,使得ap+bq=1。
只要明白了欧几里得算法,很容易就可以求出两整数的最大公约数,而这是一个小学时候就学习到的算法。这个算法有个可能让我们更熟悉的名字,叫辗转相除法。
我经常搞不清楚被除数和除数,不知道会不会有人和我一样。所以我要先在这里写明一下,防止混淆,一个除法,除号前的叫被除数,除号后的脚除数。
单次除法,X=m*Y+n,X为被除数,Y为除数,m为商,n为余数,X和Y的最大公约数等于Y和n的最大公约数。辗转相除法的每一轮除法,求最大公约数都是由求被除数、除数的最大公约数转变为被除数和玉树的最大公约数,最大公约数不变,数变小了。直到余数为0,求得最大公约数就是最一个除法下的除数。
顺便说一下,整数环具有这种相除法的结构,但不是所有的环都具有此种结构,可以做相除法的环叫欧几里得整环(Euclidean domain),给个其他的例子,比如复系数多项式环、实系数多项式环、整数系数多项式环……跑题了,就此打住。
互质的第二个定义里,如果对于互质的两个正整数p,q,p<q,我再加一个条件,要求0<a<q,那么a和b存在且唯一。这个时候,a就是q的以p为模的模乘逆元了。
可以用辗转相除法伴随逆元的生成,原理大致如下:
如果b0、b1开始做辗转相除法,步骤如下:
b0 = a0*b1 + b2
b1 = a1*b2 + b3
b2 = a2*b3 + b4
...
bn-2 = an-2*bn-1 + bn
bn-1 = an-1*bn + bn+1
bn = an*bn+1
最后一步余数为0,也就是最大公约数是bn+1,除了最后一个式子其他式子移象,余数写在左边
(1) b2 = b0 - a0*b1
(2) b3 = b1 - a1*b2
(3) b4 = b2 - a2*b3
...
(n-1) bn=bn-2 - an-2*bn-1
(n) bn+1 = bn-1 - an-1*bn
我们开始分析, (1)式可以看成是把b2表示为b0和b1的线性组合,
如果把(1)式带入(2)式,则得到把b3表示为b0和b1的线性组合,姑且称为(2.1),为了方便,把(1)给个一样的表示(1.1),
把(1.1)和(2.1)带入(3),则得到把b4表示为b0和b1的线性组合,称为(3.1),
把(2,1)和(3.1)带入(4),则得到把b5表示为b0和b1的线性组合,称为(4.1),
...
直到把bn+1表示为b0和b1的线性组合
我们这里是求逆元,如果b0和b1互质,那么bn+1应为1。
bn+1表示为b0和b1的线性组合,b1前的系数就是b1在b0模乘下的逆元了,当然该系数还要除以b0取个余数。
同样,还是写个bc程序来表示一下这个算法。
#!/usr/bin/bc -q define inv(b0, b1) { m=b0; x0 = 1; y0 = 0; x1 = 0; y1 = 1; while(1) { a = b0/b1; b = b0%b1; if(b==0) { if(b1==1) { y1 = y1 % m; if(y1<0) { y1+=m; } return y1; } else { return -1; } } /* tmp <= (x1,y1) (x1,y1) <= (x0,y0)-a(x1,y1) (x0,y0) <= tmp */ tmpx = x1; tmpy = y1; x1 = x0-a*x1; y1 = y0-a*y1; x0 = tmpx; y0 = tmpy; b0 = b1; b1 = b; } } b0 = read(); b1 = read(); c1 = inv(b0,b1) print "c1 = ",inv(b0,b1),"\n" quit
当然,算法中x0,x1是记录b0的系数,其实对于计算b1的逆元无用,所以可以省略。整个算法的平均时间复杂度为线性。
另外,此求逆算法在RSA中的应用不只在于求私钥的指数,也可用于优化模幂算法。