原理介紹
RSA 原理:
選取兩個不同的大素數p、q,並計算N=p*q,選取小素數d,並計算e,使d*e % (p-1)(q-1)=1,
對於任意A<N:
若B=A**d % N
則A=B**e % N
可見d、e形成了非對稱秘鑰關系,加密者用公鑰d加密,解密者可用私鑰e解密,第三者即使攔截了密文B、公鑰d和N,在不知道p、q的前提下,無法推算出e,從而無法獲得明文A。當N取非常大的值時,將其因式分解成p、q是非常困難的,例如當N為1024 bit時,據分析,需動用價值數千萬美金的大型計算機系統並耗費一年的時間。
RSA 密鑰的選取和加解密過程都非常簡潔,在算法上主要要實現四個問題:
1、如何處理大數運算
2、如何求解同余方程 XY % M = 1
3、如何快速進行模冪運算
4、如何獲取大素數
實際上,在實現RSA 算法的過程中大家會發現后三個問題不是各自獨立的,它們互有關聯,環環相套,相信屆時你會意識到:RSA算法是一種“優美”的算法!
大數存儲:
RSA 依賴大數運算,目前主流RSA 算法都建立在1024位的大數運算之上。而大多數的編譯器只能支持到64位的整數運算,即我們在運算中所使用的整數必須小於等於64位,即:0xffffffffffffffff,也就是18446744073709551615,這遠遠達不到RSA 的需要,於是需要專門建立大數運算庫來解決這一問題。
最簡單的辦法是將大數當作數組進行處理,數組的各元素也就是大數每一位上的數字,通常采用最容易理解的十進制數字0—9。然后對“數字數組”編寫加減乘除函數。但是這樣做效率很低,因為二進制為1024位的大數在十進制下也有三百多位,對於任何一種運算,都需要在兩個有數百個元素的數組空間上多次重循環,還要許多額外的空間存放計算的進退位標志及中間結果。另外,對於某些特殊的運算而言,采用二進制會使計算過程大大簡化,而這種大數表示方法轉化成二進制顯然非常麻煩,所以在某些實例中則干脆采用了二進制數組的方法來記錄大數,當然這樣效率就更低了。
一個有效的改進方法是將大數表示為一個n 進制數組,對於目前的32位系統而言n 可以取值為2 的32次方,即 0x100000000,假如將一個二進制為1024位的大數
轉化成0x10000000進制,就變成了32位,而每一位的取值范圍不再是二進制的0—1或十進制的0—9,而是0-0xffffffff,我們正好可以用一個32位的DWORD (如:無符號長整數,unsigned long) 類型來表示該值。所以1024位的大數就變成一個含有32個元素的 DWORD數組,而針對 DWORD數組進行各種運算所需的循環規模至多32次而已。而且0x100000000 進制與二進制,對於計算機來說,幾乎是一回事,轉換非常容易。
例如大數18446744073709551615,等於 0xffffffff ffffffff,就相當於十進制的99:有兩位,每位都是0xffffffff。而18446744073709551616等於0x00000001 00000000 00000000,就相當於十進制的100:有三位,第一位是1 ,其它兩位都是0 ,如此等等。在實際應用中,“數字數組”的排列順序采用低位在前高位在后的方式,這樣,大數A 就可以方便地用數學表達式來表示其值:
A=Sum[i=0 to n](A[i]*r**i),r=0x100000000,0<=A<r
其中Sum 表示求和,A[i]表示用以記錄A 的數組的第i 個元素,**表示乘方。
任何整數運算最終都能分解成數字與數字之間的運算,在0x100000000 進制下其“數字”最大達到0xffffffff,其數字與數字之間的運算,結果也必然超出了目前32位系統的字長。在VC++中,存在一個__int64 類型可以處理64位的整數,所以不用擔心這一問題,而在其它編譯系統中如果不存在64位整形,就需要采用更小的
進制方式來存儲大數,例如16位的WORD類型可以用來表示0x10000 進制。但效率更高的辦法還是采用32位的 DWORD類型,只不過將0x100000000 進制改成0x40000000進制,這樣兩個數字進行四則運算的最大結果為 0x3fffffff * 0x3fffffff,小於0xffffffffffffff,可以用一個雙精度浮點類型(double,52位有效數字)來儲存這一中間結果,只是不能簡單地用高位低位來將中間結果拆分成兩個“數字”。
大數加減乘除:
設有大數A、B和C,其中A>=B:
A=Sum[i=0 to p](A[i]*r**i)
B=Sum[i=0 to q](B[i]*r**i)
C=Sum[i=0 to n](C[i]*r**i)
r=0x100000000(32位),0<=A[i],B[i],C[i]<r,p>=q(A[i],B[i],C[i]都是32位的數)
則當C=A+B、C=A-B、C=A*B時,我們都可以通過計算出C來獲得C:
1.加法
C=A+B,顯然C[i]不總是等於A[i]+B[i],因為A[i]+B[i]可能>0xffffffff,而C[i]必須<=0xffffffff,這時就需要進位,當然在計算C[i-1]時也可能產生了進位,所以計算C[i]時還要加上上次的進位值。如果用一個64位變量result來記錄和(64位是為乘法准備的,實際加減法只要33位即可),另一個32位變量carry來記錄進位(為什么要32位?為乘法准備的,實際加減法進位只有1),則有:
carry=0;
for(i=0;i<=p;i++) { //i從0到p 因為A>B
result=A[i]+B[i]+carry;
C[i]=result%0x100000000 ; //從這里看result應該大於64位,至少65位
carry=result/0x100000000;
}
if(carry=0) n=p;
else n=p+1;
2.減法
C=A-B,同理C[i]不總是等於A[i]-B[i],因為A[i]-B[i]可能<0,而C[i]必須>=0,這時就需要借位,同樣在計算C[i-1]時也可能產生了借位,所以計算C時還要減去上次的借位值:
carry=0
for(i=0;i<=p;i++) { //i從0到p 因為A>B
if((A[i]-B[i]-carry)>=0){
C[i]=A[i]-B[i]-carry;
carry=0;
}
else{
C[i]=0x100000000+A[i]-B[i]-carry;
carry=1;
}
}
n=p;
while (C[n]==0) n=n-1; //將前邊的0去掉
3.乘法
C=A*B,首先我們需要觀察日常做乘法所用的“豎式計算”過程:
A3 A2 A1 A0
* B2 B1 B0
------------------------------------------
= A3B0 A2B0 A1B0 A0B0
+ A3B1 A2B1 A1B1 A0B1
+ A3B2 A2B2 A1B2 A0B2
------------------------------------------
= C5 C4 C3 C2 C1 C0
可以歸納出:C[i]=Sum[j=0 to q](A[i-j]*B[j]) (注意是C[i]),其中i-j必須>=0且<=p。
當然這一結論沒有考慮進位,雖然計算A[i-j]*B[j]和Sum的時候都可能發生進位,顯然這兩種原因產生的進位可以累加成一個進位值。最終可用如下算法完成乘法:
C = Sum[i= 0 to n](C[i]*r**i) = Sum[i= 0 to n] ( Sum[j=0 to q](A[i-j]*B[j]) *r**i).(這里的n=p+q-1,但當第n位的運算有進位時n應加1)
C也可以表示成 C= Sum[i= 0 to q](A*B[i] *r**i)
n=p+q-1
carry=0
for(i=0;i<=n;i++){
result=carry;
for(j=0;j<=q;j++){
if (0<=i-j<=p ){
result=result+A[i-j]*B[j];
C[i]=result%0x100000000;
carry=result/0x100000000;
}
}
}
if(carry!=0) {
n=n+1;
C[n]=carry
}
4.除法
對於C=A/B,由於無法將B 對A“試商”,我們只能轉換成B[q]對A[p]的試商來得到一個近似值,所以無法直接通過計算C來獲得C,只能一步步逼近C。由於:
B*(A[p]/B[q]-1)*0x100000000**(p-q)<C
令:X=0,重復A=A-X*B,X=X+(A[p]/B[q]-1)*0x100000000**(p-q),(為什么?) 直到A<B
則:X=A/B,且此時的A=A%B
注意對於任意大數A*0x100000000**k,都等價於將A 的數組中的各元素左移k 位,
不必計算;同樣,A/0x100000000**k則等價於右移。
歐幾里得方程:
在RSA 算法中,往往要在已知A、N的情況下,求 B,使得 (A*B)%N=1。即相當於求解B、M都是未知數的二元一次不定方程 A*B-N*M=1 的最小整數解。
而針對不定方程ax-by=c 的最小整數解,古今中外都進行過詳盡的研究,西方有著名的歐幾里德算法,即輾轉相除法,中國有秦九韶的“大衍求一術”。事實上二元一次不定方程有整數解的充要條件是c為a、b的最大公約數。即當c=1時,a、b必須互質。而在RSA算法里由於公鑰d為素數,素數與任何正整數互質,所以可以通
過歐幾里德方程來求解私鑰e。
歐幾里德算法是一種遞歸算法,比較容易理解:
例如:11x-49y=1,求x
(a) 11 x - 49 y = 1 49%11=5 ->
(b) 11 x - 5 y = 1 11%5 =1 ->
(c) x - 5 y = 1
令y=0 代入(c)得x=1
令x=1 代入(b)得y=2
令y=2 代入(a)得x=9
同理可使用遞歸算法求得任意 ax-by=1(a、b互質)的解。實際上通過分析歸納將遞歸算法轉換成非遞歸算法就變成了大衍求一術:
x=0,y=1
WHILE a!=0
i=y
y=x-y*a/b
x=i
i=a
a=b%a
b=i
IF x<0 x=x+b
RETURN x
模冪運算
模冪運算是RSA 的核心算法,最直接地決定了RSA 算法的性能。針對快速模冪運算這一課題,西方現代數學家提出了大量的解決方案,通常都是先將冪模運算轉化為乘模運算。
例如求D=C**15 % N,由於:a*b % n = (a % n)*(b % n) % n,所以:
C1 =C*C % N =C**2 % N
C2 =C1*C % N =C**3 % N
C3 =C2*C2 % N =C**6 % N
C4 =C3*C % N =C**7 % N
C5 =C4*C4 % N =C**14 % N
C6 =C5*C % N =C**15 % N
即:對於E=15的冪模運算可分解為6 個乘模運算,歸納分析以上方法可以發現對於任意E,都可采用以下算法計算D=C**E % N:
D=1
WHILE E>=0
IF E%2=0
C=C*C % N
E=E/2
ELSE
D=D*C % N
E=E-1
RETURN D
繼續分析會發現,要知道E 何時能整除 2,並不需要反復進行減一或除二的操作,只需驗證E 的二進制各位是0 還是1 就可以了,從左至右或從右至左驗證都可以,從左至右會更簡潔,設E=Sum[i=0 to n](E*2**i),0<=E<=1,則:
D=1
FOR i=n TO 0
D=D*D % N
IF E[i]=1 D=D*C % N
RETURN D
這樣,模冪運算就轉化成了一系列的模乘運算。
模乘運算
對於乘模運算 A*B%N,如果A、B都是1024位的大數,先計算A*B,再% N,就會產生2048位的中間結果,如果不采用動態內存分配技術就必須將大數定義中的數組空間增加一倍,這樣會造成大量的浪費,因為在絕大多數情況下不會用到那額外的一倍空間,而采用動態內存分配技術會使大數存儲失去連續性而使運算過程中的循
環操作變得非常繁瑣。所以模乘運算的首要原則就是要避免直接計算A*B。
設A=Sum[i=0 to k](A[i]*r**i),r=0x10000000,0<=A[i]<r,則:
可以用一個循環來處理:
C=0;
FOR i=n to 0
C=C*r
C=C+A[i]*B
RETURN C
這樣將一個多位乘法轉換成了一系列單位乘法和加法,由於:
a*b %n = (a%n * b%n) %n
a+b %n = (a%n + b%n) %n
所以,對於求C=A*B %N,我們完全可以在求C的過程中完成:
C=0;
FOR i=n to 0
C=C*r %N
C=C+A[i]*B %N RETURN C
這樣產生的最大中間結果是A*B 或C*r,都不超過1056(1024+32)位,空間代價會小得多,但是時間代價卻加大了,因為求模的過程由一次變成了多次。對於孤立的乘模運算而言這種時間換空間的交易還是值得的,但是對於反復循環的乘模運算,這種代價就無法承受,必須另尋出路。
蒙哥馬利模乘要解決的問題:
{注:蒙哥馬利模乘實際上是解決了這樣一個問題,即不使用除法(用移位操作)而求得模乘運算的結果。時刻注意是模運算而且是大數運算的基礎上的。}
{例如:假設進制 R=10 一個數(大數表示)23 =2*10^1+3*10^0 。欲求23 mod 5,先求 23 *10**-K mod 5的值 ,不使用乘法,我們可以采取下面的辦法 就是將 23+5*q 這時不影響模運算的結果 當23+5*q 是10的倍數時 就可以用移位操作除以10,一般k的值取23(大數)的位數(在進制R的基礎上)一直移位,最后剩下一個大於模5 小於2*5的數,在減去一個5 就是最后的結果了,這樣求出的是23 * 10^-k mod 5 的結果。我們想求23 mod 5的結果只需要先求出23*10^-k mod 5 = Z即可 在去求 z*10^k mod 5 即可 ,因為 23 % 5 = (23 * 10**K * 10 ** (-K))%5 = ((23 * 10**(-K)%5 *10 ** (-K)%5)%5。
這個例子舉得不好,因為23+5q 永遠也加不出10的倍數,但是由於RSA算法的特殊性,在蒙哥馬利模乘時是可以加出進制倍數,而實現移位的,這個例子只是領會精神,也可以不看}
蒙哥馬利模乘(下列表述不是很容易理解,可以參考上面的注)
由於RSA 的核心算法是模冪運算,模冪運算又相當於模乘運算的循環,要提高RSA 算法的效率,首要問題在於提高模乘運算的效率。不難發現,模乘過程中復雜度最高的環節是求模運算,因為一次除法實際上包含了多次加法、減法和乘法,如果在算法中能夠盡量減少除法甚至避免除法,則算法的效率會大大提高。
設A=Sum[i=0 to k](A[i]*2**i),0<=A[i]<=1,則:
C= A*B = Sum[i=0 to k](A*B*2**i) 可用循環處理為:
C=0
FOR i FROM k TO 0
C=C*2
C=C+A[i]*B
RETURN C
若令 C'= A*B *2**(-k)則
C'= Sum[i=0 to k](A[i]*B*2**(i-k))
用循環處理即:
C'=0
FOR i FROM 0 TO k
C'=C'+A[i]*B
C'=C'/2
RETURN C'
通過這一算法求A*B*2**(-k)是不精確的,因為在循環中每次除以2都可能有余數被舍棄了,但是可以通過這一算法求A*B*2**(-k) %N的精確值,方法是在對C'除
2之前,讓C'加上C'[0]*N。由於在RSA中N是兩個素數的積,總是奇數,所以當C'是奇數時,C'[0]=1,C'+C'[0]*N 就是偶數,而當C'為偶數時C'[0]=0,C'+C'[0]*N還是偶數,這樣C'/2 就不會有余數被舍棄。又因為C'+N %N = C' %N,所以在計算過程中加若干次N,並不會影響結果的正確性。可以將算法整理如下:
{A*B*2**(-k)%N 的含義是找到一個與A*B%N 同余的一個數H*2**K,再計算H%N}
C'=0
FOR i FROM 0 TO k
C'=C'+A*B
C'=C'+C'[0]*N
C'=C'/2
IF C'>=N C'=C'-N
RETURN C'
由於在RSA中A、B總是小於N,又0<=A,C'[0]<=1,所以:
C' = (C'+A*B+C'[0]*N)/2
C' < (C'+2N)/2
2C' < C'+2N
C' < 2N
既然C'總是小於2N,所以求C' %N 就可以很簡單地在結束循環后用一次減法來完成,即在求A*B*2**(-k) %N的過程中不用反復求模,達到了我們避免做除法的目
的。當然,這一算法求得的是A*B*2**(-k) %N,而不是我們最初需要的A*B %N。但是利用A*B*2**(-k)我們同樣可以求得A**E %N。
設R=2**k %N,R'=2**(-k) %N,E=Sum[i=0 to n](E[i]*2**i):
A'=A*R %N //這一步是怎么求的?
X=A'
FOR i FROM n TO 0
X=X*X*R' %N
IF E[i]=1 X=X*A'*R' %N
X=X*1*R' %N
RETURN X
最初:
X = A*R %N,
開始循環時:
X = X*X*R' %N
= A*R*A*R*R' %N
= A**2*R %N
反復循環之后:
X = A**E*R %N
最后:
X = X*1*R' %N
= A**E*R*R' %N
= A**E %N
如此,我們最終實現了不含除法的模冪算法,這就是著名的蒙哥馬利算法,而X*Y*R' %N 則被稱為“蒙哥馬利模乘”。以上討論的是蒙哥馬利模乘最簡單,最容
易理解的二進制形式。蒙哥馬利算法的核心思想在於將求A*B %N轉化為不需要反復取模的A*B*R' %N,(移位即可,因為R是2^K,總之R是與進制相關的數),但是利用二進制算法求1024位的A*B*R' %N,需要循環1024次之多,我么必然希望找到更有效的計算A*B*R' %N的算法。
{A*B%N
=(A*B*2**(-K) * 2**K )%N
=(A*B*2**(-K)%N * 2**K%N)%N
=((H*2**K * 2**(-K))%N * 2**K%N)%N
=(H%N * 2**K%N)%N //這里把A*B*2**(-K) 通過移位轉化成了 H%N 而H是一個小於2N的數,做H-N即可。
}
考慮將A表示為任意的r進制:
A = Sum[i=0 to k](A*r**i) 0<=A<=r
我們需要得到的蒙哥馬利乘積為:
C'= A*B*R' %N R'=r**(-k)
則以下算法只能得到C'的近似值
C'=0
FOR i FROM 0 TO k
C'=C'+A*B
C'=C'/r
IF C'>=N C'=C'-N
RETURN C'
因為在循環中每次C'=C'/r 時,都可能有余數被舍棄。假如我們能夠找到一個系數 q,使得(C' + A*B + q*N) %r =0,並將算法修改為:
C'=0
FOR i FROM 0 TO k
C'=C'+A*B+q*N
C'=C'/r
IF C'>=N C'=C'-N
RETURN C'
則C'的最終返回值就是A*B*R' %N的精確值,所以關鍵在於求q。由於:
(C' + A*B + q*N) %r =0
==> (C' %r + A*B %r + q*N %r) %r =0
==> (C'[0] + A*B[0] + q*N[0]) %r =0
若令N[0]*N[0]' %r =1,q=(C'[0]+A*B[0])*(r-N[0]') %r,則:
(C'[0] + A*B[0] + q*N[0]) %r
= (C'[0]+A*B[0] - (C'[0]+A*B[0])*N[0]'*N[0]) %r) %r
= 0
於是我們可以得出r為任何值的蒙哥馬利算法:
m=r-N[0]'
C'=0
FOR i FROM 0 TO k
q=(C'[0]+A*B[0])*m %r
C'=(C'+A*B+q*N)/r
IF C'>=N C'=C'-N
RETURN C'
如果令 r=0x100000000,則 %r 和 /r 運算都會變得非常容易,在1024位的運算中,循環次數k 不大於32,整個運算過程中最大的中間變量C'=(C'+A*B+q*N)
< 2*r*N < 1057位,算法效率就相當高了。唯一的額外負擔是需要計算 N[0]',使N[0]*N[0]' %r =1,而這一問題前面已經用歐幾里德算法解決過了,而且在模冪運算轉化成反復模乘運算時,N是固定值,所以N[0]'只需要計算一次,負擔並不大。
素數測試方法
數論學家利用費馬小定理研究出了多種素數測試方法,目前最快的算法是拉賓米勒測試算法,其過程如下:
(1)計算奇數M,使得N=(2**r)*M+1
(2)選擇隨機數A<N
(3)對於任意i<r,若A**((2**i)*M) % N = N-1,則N通過隨機數A的測試
(4)或者,若A**M % N = 1,則N通過隨機數A的測試
(5)讓A取不同的值對N進行5次測試,若全部通過則判定N為素數
若N 通過一次測試,則N 不是素數的概率為 25%,若N 通過t 次測試,則N 不是素數的概率為1/4**t。事實上取t 為5 時,N 不是素數的概率為 1/128,N 為素數的概率已經大於99.99%。
在實際應用中,可首先用300—500個小素數對N 進行測試,以提高拉賓米勒測試通過的概率,從而提高整體測試速度。而在生成隨機素數時,選取的隨機數最好讓r=0,則可省去步驟(3) 的測試,進一步提高測試速度。
素數測試是RSA 選取秘鑰的第一步,奇妙的是其核心運算與加解密時所需的運算完全一致:都是模冪運算。而模冪運算過程中中所需求解的歐幾里德方程又恰恰
正是選取密鑰第二步所用的運算。可見整個RSA 算法具有一種整體的和諧。