RSA简介(四)——求逆算法

　　此处所谓求逆运算，是指在模乘群里求逆。

　　第一节里提到互质的两个定义:

　　(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的逆元无用，所以可以省略。整个算法的平均时间复杂度为线性（忽略除法的时间复杂度，这里只考虑除法的个数，当然，其实除法的时间复杂度本不可忽略，至少为O(n)，最多为O(n2)，而整体的真正时间复杂度应该是除法时间复杂度在0~n上积分）。　

　　另外，此求逆算法在RSA中的应用不只在于求私钥的指数，也可用于优化模幂算法。