大數(shù)儲(chǔ)存

RSA 依賴大數(shù)運(yùn)算，目前主流RSA 算法都建立在512 到1024位的大數(shù)運(yùn)算之上。
而大多數(shù)的編譯器只能支持到64位的整數(shù)運(yùn)算，即我們?cè)谶\(yùn)算中所使用的整數(shù)必須小
于等于64位，即：0xffffffffffffffff，也就是18446744073709551615，這遠(yuǎn)遠(yuǎn)達(dá)不
到RSA 的需要，于是需要專門(mén)建立大數(shù)運(yùn)算庫(kù)來(lái)解決這一問(wèn)題。

    最簡(jiǎn)單的辦法是將大數(shù)當(dāng)作數(shù)組進(jìn)行處理，也就是將大數(shù)用0—9這十個(gè)數(shù)字組成
的數(shù)組進(jìn)行表示，然后模擬人們手工進(jìn)行“豎式計(jì)算”的過(guò)程編寫(xiě)其加減乘除函數(shù)。
但是這樣做效率很低，因?yàn)槎M(jìn)制為1024位的大數(shù)其十進(jìn)制也有三百多位，對(duì)于任何
一種運(yùn)算，都需要在兩個(gè)有數(shù)百個(gè)元素的數(shù)組空間上做多重循環(huán)，還需要許多額外的
空間存放計(jì)算的進(jìn)退位標(biāo)志及中間結(jié)果。另外，對(duì)于某些特殊的運(yùn)算而言，采用二進(jìn)
制會(huì)使計(jì)算過(guò)程大大簡(jiǎn)化，這種大數(shù)表示方法轉(zhuǎn)化成二進(jìn)制顯然非常麻煩，所以在某
些實(shí)例中則干脆采用了二進(jìn)制數(shù)組的方法來(lái)記錄大數(shù)，這樣效率就更低了。

    一個(gè)有效的改進(jìn)方法是將大數(shù)表示為一個(gè)n 進(jìn)制數(shù)組，對(duì)于目前的32位系統(tǒng)而言
n 可以取值為2 的32次方，即0x10000000，假如將一個(gè)二進(jìn)制為1024位的大數(shù)轉(zhuǎn)化成
0x10000000進(jìn)制，它就變成了32位，而每一位的取值范圍就不是二進(jìn)制的0—1或十進(jìn)
制的0—9，而是0-0xffffffff，我們正好可以用一個(gè)無(wú)符號(hào)長(zhǎng)整數(shù)來(lái)表示這一數(shù)值。
所以1024位的大數(shù)就是一個(gè)有32個(gè)元素的unsigned long數(shù)組，針對(duì)unsigned long數(shù)
組進(jìn)行各種運(yùn)算所需的循環(huán)規(guī)模至多32次而已。而且0x100000000 進(jìn)制與二進(jìn)制，對(duì)
于計(jì)算機(jī)來(lái)說(shuō)，幾乎是一回事，轉(zhuǎn)換非常容易。

    例如大數(shù)18446744073709551615，等于 ffffffff ffffffff，就相當(dāng)于十進(jìn)制的
99：有兩位，每位都是ffffffff。而18446744073709551616 等于00000001 00000000
00000000，就相當(dāng)于十進(jìn)制的100：有三位，第一位是1 ，其它兩位是0，如此等等。
在實(shí)際應(yīng)用中，“數(shù)字”數(shù)組的排列順序采用低位在前高位在后的方式，這樣，大數(shù)
A 就可以方便地用數(shù)學(xué)表達(dá)式來(lái)表示其值：A=Sum[i=0 to n](A[i]*0x100000000**i)
（其中Sum 表示求和，A[i]表示用以記錄A 的數(shù)組的第i 個(gè)元素，**表示乘方）。

    任何整數(shù)運(yùn)算最終都能分解成數(shù)字與數(shù)字之間的運(yùn)算，在0x100000000 進(jìn)制下其
“數(shù)字”最大達(dá)到0xffffffff，其數(shù)字與數(shù)字之間的運(yùn)算，結(jié)果也必然超出了目前32
系統(tǒng)的字長(zhǎng)。在VC++中，存在一個(gè)__int64 類型可以處理64位的整數(shù)，所以不用擔(dān)心
這一問(wèn)題，而在其它編譯系統(tǒng)中如果不存在64位整形，就需要采用更小的進(jìn)制方式來(lái)
存儲(chǔ)大數(shù)，例如WORD類型（16位）可以用來(lái)表示0x10000 進(jìn)制，但效率更高的辦法還
是采用32位的DWORD 類型，只不過(guò)將0x100000000 進(jìn)制改成0x40000000進(jìn)制，這樣兩
個(gè)數(shù)字進(jìn)行四則運(yùn)算的最大結(jié)果為 0x3fffffff * 0x3fffffff，小于0xffffffff，只
是不能簡(jiǎn)單地用高位低位來(lái)將運(yùn)算結(jié)果拆分成兩個(gè)“數(shù)字”。

加法

設(shè)：
  A=Sum[i=0 to p](A[i]*0x100000000**i)
  B=Sum[i=0 to q](B[i]*0x100000000**i)，p>=q
  C=Sum[i=0 to n](C[i]*0x100000000**i)=A+B

顯然：
  C[i]不是簡(jiǎn)單地等于A[i]+B[i]，因?yàn)槿绻鸆[i]>0xffffffff就需要進(jìn)位，當(dāng)然計(jì)算
  C[i-1]時(shí)也可能產(chǎn)生了進(jìn)位，所以計(jì)算C[i]時(shí)還要加上上次的進(jìn)位值。

如果用carry[i]記錄每次的進(jìn)位則有：
  C[i]=A[i]+B[i]+carry[i-1]-carry[i]*0x100000000
  其中carry[-1]=0
  若A[i]+B[i]+carry[i-1]>0xffffffff，則carry[i]=1；反之則carry[i]=0
  若carry[p]=0，則n=p；反之則n=p+1

減法

設(shè)：
  A=Sum[i=0 to p](A[i]*0x100000000**i)
  B=Sum[i=0 to q](B[i]*0x100000000**i)，p>=q
  C=Sum[i=0 to n](C[i]*0x100000000**i)=A-B

顯然：
  C[i]不是簡(jiǎn)單地等于A[i]-B[i]，因?yàn)槿绻鸄[i]<B[i]就需要借位，當(dāng)然計(jì)算
  C[i-1]時(shí)也可能產(chǎn)生了借位，所以計(jì)算C[i]時(shí)還要減去上次的借位值。

如果用carry[i]記錄每次的借位則有：
  C[i]=A[i]+carry[i]*0x100000000-B[i]-carry[i-1]
  其中carry[-1]=0
  若A[i]>B[i]則carry[i]=0；反之則carry[i]=1
  若C[p]=0，則n=p-1；反之則n=p

乘法

設(shè)：
  A=Sum[i=0 to p](A[i]*0x100000000**i)
  B=Sum[i=0 to q](B[i]*0x100000000**i)，p>=q
  C=Sum[i=0 to n](C[i]*0x100000000**i)=A*B

顯然：
  C=Sum[i=0 to q](A*B[i]*0x100000000**i)
  而(A*B[i]*100000000**i)=Sum[j=0 to p](A[j]*B[i]*0x100000000**(i+j))
  所以C=Sum[i=0 to q](Sum[j=0 to p](A[j]*B[i]*0x100000000**(i+j)))

因此：
  C[i]=Sum[j=0 to q](A[i-j]*B[j])+carry[i-1]-carry[i]*0x100000000
  其中carry[-1]=0
  carry[i]=(Sum[j=0 to q](A[i-j]*B[j])+carry[i-1])/0x100000000
  n=p+q-1，若carry[n]>0，則n=n+1，C[n]=carry

除法

設(shè)：
  A=Sum[i=0 to p](A[i]*0x100000000**i)
  B=Sum[i=0 to q](B[i]*0x100000000**i)，p>=q
  C=Sum[i=0 to n](C[i]*0x100000000**i)=A/B

由于無(wú)法將B 對(duì)A “試商”，我們只能轉(zhuǎn)換成B[q]對(duì)A[p]的試商來(lái)得到一個(gè)近似值，
所以我們不能夠直接計(jì)算C。但是，我們可以一步一步地逼近C

顯然：
  (A[p]/B[q]-1)*0x100000000**(p-q)<C

令：
  X=0

重復(fù)：
  A=A-X*B，X=X+(A[p]/B[q]-1)*0x100000000**(p-q)，直到A<B

則有：
  X=C

注意：
  由于大數(shù)可理解為0x100000000進(jìn)制，所以對(duì)于任意大數(shù)A*0x100000000**k
  都等價(jià)于將A 的數(shù)組中的各元素左移k 位，不必計(jì)算；同樣，除法則等價(jià)于右移

取模

設(shè)：
  A=Sum[i=0 to p](A[i]*0x100000000**i)
  B=Sum[i=0 to q](B[i]*0x100000000**i)，p>=q
  C=Sum[i=0 to n](C[i]*0x100000000**i)=A%B

求模與求商的過(guò)程一致，只是由于不需要記錄商而更加簡(jiǎn)單：

重復(fù)：
  A=A-(A[p]/B[q]-1)*0x100000000**(p-q)*B，直到A<B

則有：
  A=C

二元一次方程

在RSA 算法中，往往要在已知A、M的情況下，求 B，使得 (A*B)%M=1。即相當(dāng)于
求解B、N都是未知數(shù)的二元一次不定方程 A*B-M*N=1，的最小整數(shù)解。

    而針對(duì)不定方程ax-by=1 的最小整數(shù)解，古今中外都進(jìn)行過(guò)詳盡的研究，西方有
著名的歐幾里德算法，即輾轉(zhuǎn)相除法，中國(guó)有秦九韶的“大衍求一術(shù)”。歐幾里德算
法是一種遞歸算法，比較容易理解：

    例如：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互質(zhì)）的解，實(shí)際上通過(guò)分析歸納
將遞歸算法轉(zhuǎn)換成非遞歸算法就變成了大衍求一術(shù)。

冪模運(yùn)算

冪模運(yùn)算是RSA 核心算法，最直接地決定了RSA 算法的性能，針對(duì)快速冪模運(yùn)算
這一課題，許多西方現(xiàn)代數(shù)學(xué)家提出了大量的解決方案。通常都是先將冪模運(yùn)算化簡(jiǎn)
為乘模運(yùn)算。

    例如求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

    即：
        對(duì)于E=15的冪模運(yùn)算可分解為6個(gè)乘模運(yùn)算

    歸納分析以上方法可以發(fā)現(xiàn)對(duì)于任意E，可采用以下算法計(jì)算D=C**E % N：
    D=1
    WHILE E>=0
      IF E為奇數(shù)
        D=D*C % N
        D=D*D % N
        E=E-1
      IF E為偶數(shù)
        D=D*D % N
        E=E/2
    RETURN D

    再加以分析會(huì)發(fā)現(xiàn)，要知道D 何時(shí)需乘 C，不需要反復(fù)對(duì)E 進(jìn)行減一或除二的操
作，只需要驗(yàn)證E 的二進(jìn)制各位是0 還是1 就可以了，而且從左至右驗(yàn)證和從右至左
驗(yàn)證都行，反而從左至右驗(yàn)證更簡(jiǎn)單：

    若E=Sum[i=0 to n](E[i]*2**i)，0<=E[i]<=1
    D=1
    FOR i=n TO 0
      D=D*D % N
      IF E[i]=1
        D=D*C % N
    RETURN D

    剩下的問(wèn)題就是乘模運(yùn)算了，對(duì)于A*B % N，如果A、B 都是1024位的大數(shù)，先計(jì)
算A*B，再% N，就會(huì)產(chǎn)生2048位的中間結(jié)果，如果不采用動(dòng)態(tài)內(nèi)存分配技術(shù)就必須將
大數(shù)定義中的數(shù)組空間增加一倍，這樣會(huì)造成大量的浪費(fèi)，因?yàn)樵诮^大多數(shù)情況下不
會(huì)用到那額外的一倍空間，而采用動(dòng)態(tài)內(nèi)存分配技術(shù)會(huì)使大數(shù)存儲(chǔ)失去連續(xù)性而使運(yùn)
算過(guò)程中的循環(huán)操作變得非常繁瑣。所以計(jì)算的首要原則就是要避免計(jì)算A*B。

    由于：
        A*B=A*(Sum[i=0 to n](B[i]*0x100000000**i))

    所以：
        A*B % N = (Sum[i=0 to n]((A*B[i])*0x100000000**i)) % N

    可以用一個(gè)循環(huán)求得：
    C=0;
    FOR i=0 to n
      C=C+A*B[i]*0x100000000 % N
    RETURN C

    事實(shí)上，有一種蒙哥馬利算法能夠更快地完成多次循環(huán)的乘模運(yùn)算，但是其原理
涉及較多的數(shù)論知識(shí)，且實(shí)現(xiàn)起來(lái)比較麻煩，對(duì)速度雖有提高，經(jīng)測(cè)試也不會(huì)超過(guò)一
個(gè)數(shù)量級(jí)，所以暫且不予考慮。

素?cái)?shù)測(cè)試

    數(shù)論學(xué)家利用費(fèi)馬小定理研究出了多種素?cái)?shù)測(cè)試方法，目前最快的算法是拉賓米
勒測(cè)試算法，其過(guò)程如下：

（1）計(jì)算奇數(shù)M，使得N=(2**r)*M+1
（2）選擇隨機(jī)數(shù)A<N
（3）對(duì)于任意i<r，若A**((2**i)*M) MOD N = N-1，則N通過(guò)隨機(jī)數(shù)A的測(cè)試
（4）或者，若A**M MOD N = 1，則N通過(guò)隨機(jī)數(shù)A的測(cè)試
（5）讓A取不同的值對(duì)N進(jìn)行5次測(cè)試，若全部通過(guò)則判定N為素?cái)?shù)

    若N 通過(guò)一次測(cè)試，則N 不是素?cái)?shù)的概率為 25%，若N 通過(guò)t 次測(cè)試，則N 不是
素?cái)?shù)的概率為1/4**t。事實(shí)上取t 為5 時(shí)，N 不是素?cái)?shù)的概率為 1/128，N 為素?cái)?shù)的
概率已經(jīng)大于99.99%。

    在實(shí)際應(yīng)用中，可首先用300—500個(gè)小素?cái)?shù)對(duì)N 進(jìn)行測(cè)試，以提高拉賓米勒測(cè)試
通過(guò)的概率，從而提高測(cè)試速度。而在生成隨機(jī)素?cái)?shù)時(shí)，選取的隨機(jī)數(shù)最好讓 r=0，
則可省去步驟（3） 的測(cè)試，進(jìn)一步提高測(cè)試速度。

輸入輸出

    大數(shù)的輸入輸出是通過(guò)字符串來(lái)完成的，事實(shí)上很容易實(shí)現(xiàn)，例如按照十進(jìn)制格
式進(jìn)行處理，則：

    輸入：
    X=0
    FOR i=0 TO n
      X=X*10
      X=X+(int)(str[n]-48)
    RETURN X