版權申明：本文為博主窗戶(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),"\n"
quit
當然，算法中x0,x1是記錄b0的系數，其實對於計算b1的逆元無用，所以可以省略。整個算法的平均時間復雜度為線性（忽略除法的時間復雜度，這里只考慮除法的個數，當然，其實除法的時間復雜度本不可忽略，至少為O(n)，最多為O(n2)，而整體的真正時間復雜度應該是除法時間復雜度在0~n上積分）。
另外，此求逆算法在RSA中的應用不只在於求私鑰的指數，也可用於優化模冪算法。