==========================================================================
前言：此文來自於www.pediy.com一位Cracker---afanty之手。他建立了一個VC++（MFC）版的大數運算庫。用Delphi的朋友們請到http://ace.ulyssis.student.kuleuven.ac.be/~triade/去下載Pascal的GInt大數運算庫。當然，所有基於大素數分解難題的公開密匙演算法都可以使用這兩個庫。不過請你小心裡面可能存在的Bug...

--- Crossbow
==========================================================================

基於以下原因，俺估計RSA演算法會被越來越多的共享軟體採用：

1、原理簡潔易懂
2、加密強度高
3、專利限制已過期
4、看雪老大三番五次呼籲共享軟體採用成熟的非對稱加密技術

所以，大家應該對RSA演算法進行深入瞭解。

俺寫了一個簡單易懂的大數運算庫，並使用該庫作了一個RSA Demo，大家可以使用這一
Demo生成真正隨機的、各種長度的RSA 金鑰對。其中生成1024位的RSA 金鑰耗時不超過
五分鐘，而對1024位以內的密文進行解密則不超過三秒鐘，應該是可以接受的。

有一點需要說明的是，假如類似於這個Demo的RSA工具被共享軟體作者廣泛用於註冊碼
的生成與驗證，俺認為Cracker們的日子就會過得很無趣了，唉！

RSA 依賴大數運算，目前主流RSA 演算法都建立在512 位到1024位的大數運算之上，所以
我們在現階段首先需要掌握1024位的大數運算原理。

大多數的編譯器只能支援到64位的整數運算，即我們在運算中所使用的整數必須小於等
於64位，即：0xffffffffffffffff，也就是18446744073709551615，這遠遠達不到RSA
的需要，於是需要專門建立大數運算庫來解決這一問題。

最簡單的辦法是將大數當作字串進行處理，也就是將大數用10進位制字元陣列進行表示，
然後模擬人們手工進行“豎式計算”的過程編寫其加減乘除函式。但是這樣做效率很低，
因為1024位的大數其10進位制數字個數就有數百個，對於任何一種運算，都需要在兩個有
數百個元素的陣列空間上做多重迴圈，還需要許多額外的空間存放計算的進位退位標誌
及中間結果。當然其優點是演算法符合人們的日常習慣，易於理解。

另一種思路是將大數當作一個二進位制流進行處理，使用各種移位和邏輯操作來進行加減
乘除運算，但是這樣做程式碼設計非常複雜，可讀性很低，難以理解也難以除錯。
於是俺琢磨了一種介於兩者之間的思路：

將大數看作一個n進位制陣列，對於目前的32位系統而言n可以取值為2的32次方，即0x100
00000，假如將一個1024位的大數轉化成0x10000000 進位制，它就變成了32位，而每一位
的取值範圍就不是0-1或0-9，而是0-0xffffffff。我們正好可以用一個無符號長整數來
表示這一數值。所以1024位的大數就是一個有32個元素的unsigned long陣列。而且0x1
00000000進位制的陣列排列與2 進位制流對於計算機來說，實際上是一回事，但是我們完全
可以針對unsigned long 陣列進行“豎式計算”，而迴圈規模被降低到了32次之內，並
且演算法很容易理解。

例如大數18446744073709551615，等於“ffffffff ffffffff”，它就相當於10
進位制的“99”：有兩位，每位都是ffffffff。而大數18446744073709551616，等於“00000001
00000000 00000000”，它就相當於10進位制的“100”：有三位，第一位是1 ，其它兩位
是0。如果我們要計算18446744073709551616-18446744073709551615，就類似於100-99：

00000001 00000000 00000000
- ffffffff ffffffff
-----------------------------
= 0 0 1
當然，因為在VC裡面存在一個__int64型別可以用來計算進位與借位值，所以將大數當作
0x100000000進位制進行運算是可能的，而在其他編譯系統中如果不存在64位整形，則可以
採用0x40000000進位制，由於在0x40000000進制中，對任何兩個“數字”進行四則運算，
結果都在0x3fffffff*03fffffff之間，小於0xffffffff，都可以用一個32位無符號整數
來表示。

/****************************************************************/
file://大數運算庫標頭檔案：BigInt.h
file://作者：[email protected]
file://版本：1.0 (2003.4.26)
file://說明：適用於MFC
/****************************************************************/

#define BI_MAXLEN 40
#define DEC 10
#define HEX 16

class CBigInt
{
public:
int m_nSign; file://記錄大數的符號，支援負值運算
int m_nLength; file://記錄0x10000000進位制的位數，0-40之間，相當於2進位制的0-1280位
unsigned long m_ulvalue[BI_MAXLEN]; file://記錄每一位的“數字”

CBigInt();
~CBigInt();

file://將大數賦值為另一個大數
CBigInt& Mov(CBigInt& A);

file://將大數賦值為編譯器能夠理解的任何整形常數或變數
CBigInt& Mov(unsigned __int64 A);

file://比較兩個大數大小
int Cmp(CBigInt& A);

file://計算兩個大數的和
CBigInt Add(CBigInt& A);

file://重載函式以支援大數與普通整數相加
CBigInt Add(long A);

file://計算兩個大數的差
CBigInt Sub(CBigInt& A);

file://重載函式以支援大數與普通整數相減
CBigInt Sub(long A);

file://計算兩個大數的積
CBigInt Mul(CBigInt& A);

file://重載函式以支援大數與普通整數相乘
CBigInt Mul(long A);

file://計算兩個大數的商
CBigInt Div(CBigInt& A);

file://重載函式以支援大數與普通整數相除
CBigInt Div(long A);

file://計算兩個大數相除的餘數
CBigInt Mod(CBigInt& A);

file://重載函式以支援大數與普通整數相除求模
long Mod(long A);

file://將輸入的10進位制或16進位制字串轉換成大數
int InPutFromStr(CString& str, const unsigned int system);

file://將大數按10進位制或16進位制格式輸出到字串
int OutPutToStr(CString& str, const unsigned int system);

file://歐幾裡德演算法求：Y=X.Euc(A)，使滿足：YX mod A = 1
CBigInt Euc(CBigInt& A);

file://蒙哥馬利演算法求：Y=X.Mon(A,B)，使滿足：X^A mod B = Y
CBigInt Mon(CBigInt& A, CBigInt& B);
};

注意以上函式的宣告格式，完全遵循普通整數運算的習慣，例如大數
Y=X+Z 相當於 Y.Mov(X.(Add(Z))，這樣似乎沒有Mov(Y,Add(X,Z))
看起來舒服，但是一旦我們過載運算子“=”為“Mov”，“+”為“Add”，
則Y.Mov(X.(Add(Z))的形式就等價於 Y=X+Z。

俺不知道其他程式語言裡是否支援運算浮過載，至少這樣定義函式格式
在C++裡可以帶來很大的方便。

下面讓我們來實現大數類的主要成員函式：
/****************************************************************/
file://大數運算庫原始檔：BigInt.cpp
file://作者：[email protected]
file://版本：1.0 (2003.4.26)
file://說明：適用於MFC
/****************************************************************/

#include "stdafx.h"
#include "BigInt.h"

file://初始化大數為0
CBigInt::CBigInt()
{
m_nSign=1;
m_nLength=1;
for(int i=0;i<BI_MAXLEN;i++)m_ulvalue=0;
}

file://採用預設的解構函式
CBigInt::~CBigInt()
{
}

file://大數比較，如果大數A位數比大數B多，當然A>B
file://如果位數相同，則從高位開始比較，直到分出大小
int CBigInt::Cmp(CBigInt& A)
{
if(m_nLength>A.m_nLength)return 1;
if(m_nLength<A.m_nLength)return -1;
for(int i=m_nLength-1;i>=0;i--)
{
if(m_ulvalue>A.m_ulvalue)return 1;
if(m_ulvalue<A.m_ulvalue)return -1;
}
return 0;
}

file://照搬引數的各屬性
CBigInt& CBigInt::Mov(CBigInt& A)
{
m_nLength=A.m_nLength;
for(int i=0;i<BI_MAXLEN;i++)m_ulvalue=A.m_ulvalue;
return *this;
}

file://大數相加
file://調用形式：N.Add(A)，返回值：N+A
file://若兩大數符號相同，其值相加，否則改變引數符號再呼叫大數相減函式
/******************************************************************/
例如：
A B C
+ D E
--------------
= S F G H

其中，若C+E<=0xffffffff，則H=C+E，carry(進位標誌)=0
若C+E>0xffffffff，則H=C+E-0x100000000，carry=1

若B+D+carry<=0xfffffff，則G=B+D，carry=0
若B+D+carry>0xfffffff，則G=B+D+carry-0x10000000，carry=1

若carry=0，則F=A，S=0
若carry=1，A<0xfffffff，則F=A+1，S=0
若carry=1，A=0xfffffff，則F=0，S=1
/*****************************************************************/
CBigInt CBigInt::Add(CBigInt& A)
{
CBigInt X;
if(X.m_nSign==A.m_nSign)
{
X.Mov(*this);
int carry=0;
unsigned __int64 sum=0;
if(X.m_nLength<A.m_nLength)X.m_nLength=A.m_nLength;
for(int i=0;i<X.m_nLength;i++)
{
sum=A.m_ulvalue;
sum=sum+X.m_ulvalue+carry;
X.m_ulvalue=(unsigned long)sum;
if(sum>0xffffffff)carry=1;
else carry=0;
}
if(X.m_nLength<BI_MAXLEN)
{
X.m_ulvalue[X.m_nLength]=carry;
X.m_nLength+=carry;
}
return X;
}
else{X.Mov(A);X.m_nSign=1-X.m_nSign;return Sub(X);}
}

file://大數相減
file://調用形式：N.Sub(A)，返回值：N-A
file://若兩大數符號相同，其值相減，否則改變引數符號再呼叫大數相加函式
/******************************************************************/
例如：
A B C
- D E
--------------
= F G H

其中，若C>=E，則H=C-E，carry(借位標誌)=0
若C<E，則H=C-E+0x100000000，carry=1

若B-carry>=D，則G=B-carry-D，carry=0
若B-carry<D，則G=B-carry-D+0x10000000，carry=1

若carry=0，則F=A
若carry=1，A>1，則F=A-1
若carry=1，A=1，則F=0
/*****************************************************************/
CBigInt CBigInt::Sub(CBigInt& A)
{
CBigInt X;
if(m_nSign==A.m_nSign)
{
X.Mov(*this);
int cmp=X.Cmp(A);
if(cmp==0){X.Mov(0);return X;}
int len,carry=0;
unsigned __int64 num;
unsigned long *s,*d;
if(cmp>0)
{
s=X.m_ulvalue;
d=A.m_ulvalue;
len=X.m_nLength;
}
if(cmp<0)
{
s=A.m_ulvalue;
d=X.m_ulvalue;
len=A.m_nLength;
X.m_nSign=1-X.m_nSign;
}
for(int i=0;i<len;i++)
{
if((s-carry)>=d)
{
X.m_ulvalue=s-carry-d;
carry=0;
}
else
{
num=0x100000000+s;
X.m_ulvalue=(unsigned long)(num-carry-d);
carry=1;
}
}
while(X.m_ulvalue[len-1]==0)len--;
X.m_nLength=len;
return X;
}
else{X.Mov(A);X.m_nSign=1-X.m_nSign;return Add(X);}
}

file://大數相乘
file://調用形式：N.Mul(A)，返回值：N*A
/******************************************************************/
例如：
A B C
* D E
----------------
= S F G H
+ T I J K
----------------
= U V L M N

其中，SFGH=ABC*E，TIJK=ABC*D

而對於：
A B C
* E
-------------
= S F G H

其中，若C*E<=0xffffffff，則H=C*E，carry(進位標誌)=0
若C*E>0xffffffff，則H=(C*E)&0xffffffff
carry=(C*E)/0xffffffff
若B*E+carry<=0xffffffff，則G=B*E+carry，carry=0
若B*E+carry>0xffffffff，則G=(B*E+carry)&0xffffffff
carry=(B*E+carry)/0xffffffff
若A*E+carry<=0xffffffff，則F=A*E+carry，carry=0
若A*E+carry>0xffffffff，則F=(A*E+carry)&0xffffffff
carry=(A*E+carry)/0xffffffff
S=carry
/*****************************************************************/
CBigInt CBigInt::Mul(CBigInt& A)
{
CBigInt X,Y;
unsigned __int64 mul;
unsigned long carry;
for(int i=0;i<A.m_nLength;i++)
{
Y.m_nLength=m_nLength;
carry=0;
for(int j=0;j<m_nLength;j++)
{
mul=m_ulvalue[j];
mul=mul*A.m_ulvalue+carry;
Y.m_ulvalue[j]=(unsigned long)mul;
carry=(unsigned long)(mul>>32);
}
if(carry&&(Y.m_nLength<BI_MAXLEN))
{
Y.m_nLength++;
Y.m_ulvalue[Y.m_nLength-1]=carry;
}
if(Y.m_nLength<BI_MAXLEN-i)
{
Y.m_nLength+=i;
for(int k=Y.m_nLength-1;k>=i;k--)Y.m_ulvalue[k]=Y.m_ulvalue[k-i];
for(k=0;k<i;k++)Y.m_ulvalue[k]=0;
}
X.Mov(X.Add(Y));
}
if(m_nSign+A.m_nSign==1)X.m_nSign=0;
else X.m_nSign=1;
return X;
}

file://大數相除
file://調用形式：N.Div(A)，返回值：N/A
file://除法的關鍵在於“試商”，然後就變成了乘法和減法
file://這裡將被除數與除數的試商轉化成了被除數最高位與除數最高位的試商
CBigInt CBigInt:iv(CBigInt& A)
{
CBigInt X,Y,Z;
int len;
unsigned __int64 num,div;
unsigned long carry=0;
Y.Mov(*this);
while(Y.Cmp(A)>0)
{
if(Y.m_ulvalue[Y.m_nLength-1]>A.m_ulvalue[A.m_nLength-1])
{
len=Y.m_nLength-A.m_nLength;
div=Y.m_ulvalue[Y.m_nLength-1]/(A.m_ulvalue[A.m_nLength-1]+1);
}
else if(Y.m_nLength>A.m_nLength)
{
len=Y.m_nLength-A.m_nLength-1;
num=Y.m_ulvalue[Y.m_nLength-1];
num=(num<<32)+Y.m_ulvalue[Y.m_nLength-2];
if(A.m_ulvalue[A.m_nLength-1]==0xffffffff)div=(num>>32);
else div=num/(A.m_ulvalue[A.m_nLength-1]+1);
}
else
{
X.Mov(X.Add(1));
break;
}
Z.Mov(div);
Z.m_nLength+=len;
for(int i=Z.m_nLength-1;i>=len;i--)Z.m_ulvalue=Z.m_ulvalue[i-len];
for(i=0;i<len;i++)Z.m_ulvalue=0;
X.Mov(X.Add(Z));
Z.Mov(Z.Mul(A));
Y.Mov(Y.Sub(Z));
}
if(Y.Cmp(A)==0)X.Mov(X.Add(1));
if(m_nSign+A.m_nSign==1)X.m_nSign=0;
else X.m_nSign=1;
return X;
}

file://大數求模
file://調用形式：N.Mod(A)，返回值：N%A
file://求模與求商原理相同
CBigInt CBigInt::Mod(CBigInt& A)
{
CBigInt X,Y;
int len;
unsigned __int64 num,div;
unsigned long carry=0;
X.Mov(*this);
while(X.Cmp(A)>0)
{
if(X.m_ulvalue[X.m_nLength-1]>A.m_ulvalue[A.m_nLength-1])
{
len=X.m_nLength-A.m_nLength;
div=X.m_ulvalue[X.m_nLength-1]/(A.m_ulvalue[A.m_nLength-1]+1);
}
else if(X.m_nLength>A.m_nLength)
{
len=X.m_nLength-A.m_nLength-1;
num=X.m_ulvalue[X.m_nLength-1];
num=(num<<32)+X.m_ulvalue[X.m_nLength-2];
if(A.m_ulvalue[A.m_nLength-1]==0xffffffff)div=(num>>32);
else div=num/(A.m_ulvalue[A.m_nLength-1]+1);
}
else
{
X.Mov(X.Sub(A));
break;
}
Y.Mov(div);
Y.Mov(Y.Mul(A));
Y.m_nLength+=len;
for(int i=Y.m_nLength-1;i>=len;i--)Y.m_ulvalue=Y.m_ulvalue[i-len];
for(i=0;i<len;i++)Y.m_ulvalue=0;
X.Mov(X.Sub(Y));
}
if(X.Cmp(A)==0)X.Mov(0);
return X;
}

file://暫時只給出了十進位制字串的轉化
int CBigInt::InPutFromStr(CString& str, const unsigned int system=DEC)
{
int len=str.GetLength();
Mov(0);
for(int i=0;i<len;i++)
{
Mov(Mul(system));
int k=str-48;
Mov(Add(k));
}
return 0;
}

file://暫時只給出了十進位制字串的轉化
int CBigInt::OutPutToStr(CString& str, const unsigned int system=DEC)
{
str="";
char ch;
CBigInt X;
X.Mov(*this);
while(X.m_ulvalue[X.m_nLength-1]>0)
{
ch=X.Mod(system)+48;
str.Insert(0,ch);
X.Mov(X.Div(system));
}
return 0;
}

file://歐幾裡德演算法求：Y=X.Euc(A)，使滿足：YX mod A=1
file://相當於對不定方程ax-by=1求最小整數解
file://實際上就是初中學過的輾轉相除法
/********************************************************************/
例如：11x-49y=1，求x

11 x - 49 y = 1 a)
49%11=5 -> 11 x - 5 y = 1 b)
11%5 =1 -> x - 5 y = 1 c)

令y=1 代入c)式 得x=6
令x=6 代入b)式 得y=13
令y=13 代入a)式 得x=58
/********************************************************************/
CBigInt CBigInt::Euc(CBigInt& A)
{
CBigInt X,Y;
X.Mov(*this);
Y.Mov(A);
if((X.m_nLength==1)&&(X.m_ulvalue[0]==1))return X;
if((Y.m_nLength==1)&&(Y.m_ulvalue[0]==1)){X.Mov(X.Sub(1));return X;}
if(X.Cmp(Y)==1)X.Mov(X.Mod(Y));
else Y.Mov(Y.Mod(X));
X.Mov(X.Euc(Y));
Y.Mov(*this);
if(Y.Cmp(A)==1)
{
X.Mov(X.Mul(Y));
X.Mov(X.Sub(1));
X.Mov(X.Div(A));
}
else
{
X.Mov(X.Mul(A));
X.Mov(X.Add(1));
X.Mov(X.Div(Y));
}
return X;
}

file://蒙哥馬利演算法求：Y=X.Mon(A,B)，使滿足：X^A mod B=Y
file://俺估計就是高中學過的反覆平方法
CBigInt CBigInt::Mon(CBigInt& A, CBigInt& B)
{
CBigInt X,Y,Z;
X.Mov(1);
Y.Mov(*this);
Z.Mov(A);
while((Z.m_nLength!=1)||Z.m_ulvalue[0])
{
if(Z.m_ulvalue[0]&1)
{
Z.Mov(Z.Sub(1));
X.Mov(X.Mul(Y));
X.Mov(X.Mod(B));
}
else
{
Z.Mov(Z.Div(2));
Y.Mov(Y.Mul(Y));
Y.Mov(Y.Mod(B));
}
}
return X;
}

最後需要說明的是因為在VC裡面存在一個__int64型別可以
用來計算進位與借位值，所以將大數當作0x100000000進位制
進行運算是可能的，而在其他編譯系統中如果不存在64位
整形，則可以採用0x40000000進位制，由於在0x40000000
進制中，對任何兩個“數字”進行四則運算，結果都在
0x3fffffff*03fffffff之間，小於0xffffffff，都可以用
一個32位無符號整數來表示。事實上《楚漢棋緣》採用的
freelip大數庫正是運用了0x40000000進位制來表示大數的，
所以其反彙編後大數的值在記憶體中表現出來有些“奇怪”。