版权申明：本文为博主窗户(Colin Cai)原创，欢迎转帖。如要转贴，必须注明原文网址 http://www.cnblogs.com/Colin-Cai/p/7354682.html 作者：窗户 QQ：6679072 E-mail：6679072@qq.com
此处所谓求逆运算，是指在模乘群里求逆。
第一节里提到互质的两个定义:
(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),"
"
quit
当然，算法中x0,x1是记录b0的系数，其实对于计算b1的逆元无用，所以可以省略。整个算法的平均时间复杂度为线性（忽略除法的时间复杂度，这里只考虑除法的个数，当然，其实除法的时间复杂度本不可忽略，至少为O(n)，最多为O(n2)，而整体的真正时间复杂度应该是除法时间复杂度在0~n上积分）。
另外，此求逆算法在RSA中的应用不只在于求私钥的指数，也可用于优化模幂算法。