Diffie-Hellman金鑰交換演算法的原理及程式示範

在http://en.wikipedia.org/wiki/Diffie-Hellman上面給出了這個金鑰交換協議的曆史，原理，重要文獻的連結，以及示範代碼。它的數學基礎就是離散對數這個數學難題。用它進行金鑰交換的過程簡述如下：
選取兩個大數p和g並公開，其中p是一個素數，g是p的一個模p本原單位根(primitive root module p)，所謂本原單位根就是指在模p乘法運算下，g的1次方，2次方……(p-1)次方這p-1個數互不相同，並且取遍1到p-1；
對於Alice(其中的一個通訊者)，隨機產生一個整數a，a對外保密，計算Ka = g^a mod p，將Ka發送給Bob；
對於Bob(另一個通訊者)，隨機產生一個整數b，b對外保密，計算Kb = g^b mod p，將Kb發送給Alice；
在Alice方面，收到Bob送來的Kb後，計算出密鑰為：key = Kb^a mod p=(g^b)^a=g^(b*a) mod p；
對於Bob，收到Alice送來的Ka後，計算出密鑰為：key = Ka ^ b mod p=(g^a)^b=g^(a*b) mod p。
攻擊者知道p和g，並且截獲了Ka和Kb，但是當它們都是非常大的數的時候，依靠這四個數來計算a和b非常困難，這就是離散對數數學難題。
要實現Diffie-Hellman金鑰交換協議，需要能夠快速計算大數模冪，在模冪演算法中，仍需計算大數的乘法和模運算，所以整個過程需要三個演算法：高精度乘法，高精度除法(用來同時求出一個大數除以另一個大數的商和餘數)，快速模冪演算法。
高精度的乘法和除法可以程式類比手算。快速模冪演算法也是從手算中總結出規律來，例如：
5^8 = (5^2)^4 = (25)^4 = (25^2)^2 = (625)^2，這樣，原來計算5^8需要做8次乘法，而現在則只需要三次乘法，分別是：5^2, 25^2, 625^2。這就是快速模冪演算法的基礎。將演算法描述出來，那就是：
演算法M:輸入整數a,b,p，計算a^b mod p：
M1.初始化c = 1
M2.如果b為0，　則c就是所要計算的結果。返回c的值。演算法結束。
M3.如果b為奇數，則令c = c *a mod p，　令b = b-1，轉到M2。
M4.如果b為偶數，則令a = a * a mod p，　令b = b / 2，轉到M2。
高精度試除法原理簡單，但是代碼實現起來需要仔細考慮一些細節。
我的示範代碼如下：
高精度運算類：

1. class SuperNumber {
2. public:
3.     SuperNumber() {
4.         memset(data, 0, MAX_SIZE);
5.         high = 0;
6.     }
7.     // 一般整型到SuperNumber的轉換，該版本中不支援負數
8.     SuperNumber(unsigned long l) {
9.         memset(data, 0, MAX_SIZE);
10.         high = 0;
11.         while(l) {
12.             data[++high] = l % 10;
13.             l /= 10;
14.         }
15.     }
16.     // str為字串形式表示的十進位數
17.     SuperNumber(const char* str) {
18.         assert(str != NULL);
19.         high = strlen(str);
20.         for(int i = high, j = 0; i >= 1; i--, j++) {
21.             data[i] = str[j] - '0';
22.         }
23.     }
24.     SuperNumber(const SuperNumber& s) {
25.         memcpy(data, s.data, MAX_SIZE);
26.         high = s.high;
27.     }
28.     operator const char*() const {
29.         return toString(10);
30.     }
31.     SuperNumber& operator=(const SuperNumber& s) {
32.         if(this != &s) {
33.             memcpy(data, s.data, MAX_SIZE);
34.             high = s.high;
35.         }
36.         return *this;
37.     }
39.     // 將資料置為0
40.     void reset() {
41.         memset(data, 0, MAX_SIZE);
42.         high = 0;
43.     }
45.     // str為字串形式表示的十進位數
46.     void setToStr(const char* str) {
47.         assert(str != NULL);
48.         high = strlen(str);
49.         for(int i = high, j = 0; i >= 1; i--, j++) {
50.             data[i] = str[j] - '0';
51.         }
52.     }
53.     // 將資料轉換成以base指定的進位的字串形式，預設為十進位
54.     const char* toString(int base = 10) const {
55.         static char buf[MAX_SIZE];
56.         const char table[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
58.         if(high == 0) return "0";
60.         assert(base >= 2);      // 指定的進位應該不小於2
62.         // 進位轉換
63.         buf[MAX_SIZE-1] = '/0';
64.         int begin = MAX_SIZE-1;
65.         char temp[MAX_SIZE];
66.         memcpy(temp, data, MAX_SIZE);
67.         while(1) {
68.             // 找最高位的起始位置
69.             int h = high;
70.             while(temp[h] == 0 && h >= 1) h--;
71.             if(h == 0) break;
73.             // 除基取餘
74.             int t = 0;
75.             while(h >= 1) {
76.                 t = t * 10 + temp[h];
77.                 temp[h] = t / base;
78.                 t = t % base;
79.                 h--;
80.             }
81.             buf[--begin] = table[t];
82.         }
84.         return buf+begin;
85.     }
87.     // 乘法
88.     SuperNumber operator*(const SuperNumber& s) const {
89.         SuperNumber result;     // default set to 0
90.         int i, j;
91.         
92.         // 相乘
93.         for(i = 1; i <= high; i++) {
94.             for(j = 1; j <= s.high; j++) {
95.                 int k = data[i] * s.data[j] + result.data[i+j-1];
96.                 result.data[i+j-1] = k % 10;
97.                 result.data[i+j] += k / 10;
98.             }
99.         }
100.         // 進位
101.         for(i = 1; i < MAX_SIZE - 1; i++) {
102.             if(result.data[i] >= 10) {
103.                 result.data[i+1] += result.data[i] / 10;
104.                 result.data[i] %= 10;
105.             }
106.         }
107.         // 確定最高位
108.         for(i = MAX_SIZE-1; i >= 1 && result.data[i] == 0; i--);
109.         result.high = i;
111.         return result;
112.     }
114.     // 除法，內部調用doDivide來實現
115.     SuperNumber operator/(const SuperNumber& s) const {
116.         SuperNumber q, r;
117.         doDivide(s, q, r);
118.         return q;
119.     }
120.     // 模運算，內部調用doDivide來實現
121.     SuperNumber operator%(const SuperNumber& s) const {
122.         SuperNumber q, r;
123.         doDivide(s, q, r);
124.         return r;
125.     }
127.     // 除法運算，一次除法運算中同時得到商和餘數，運算子/和%的重載
128.     // 內部會調用這個函數，dest為除數，Q為商，R為餘數，演算法使用試除法
129.     void doDivide(const SuperNumber& dest, SuperNumber& Q, SuperNumber& R) const {
130.         int i, j, t;
132.         Q.reset();
133.         Q.high = high - dest.high + 1; // 商的初始位元
134.         R = *this;                     // 餘數初始實為被除數
135.         t = dest.high;
136.         // 判斷除法是否結束
137.         while(R >= dest) {
138.             // 迴圈相減進行試除
139.             while(dest >= R.sub(1, t)) {
140.                 Q.data[Q.high--] = 0;
141.                 ++t;
142.             }
143.             while(R.sub(1, t) >= dest) {
144.                 // i為相減時被除數最低下標，j為除數最低下標
145.                 for(i=R.high-t+1,j=1; j<=dest.high; i++,j++) {
146.                     R.data[i] -= dest.data[j];
147.                     if(R.data[i] < 0) {
148.                         R.data[i] += 10;
149.                         R.data[i+1] -= 1;
150.                     }
151.                 }
152.                 while(R.data[i] < 0 && i <= R.high) {
153.                     R.data[i] += 10;
154.                     R.data[i+1] -= 1;
155.                     ++i;
156.                 }
157.                 Q.data[Q.high] += 1;
158.             }
159.             // 一次試除結束，更新商的最高位下標
160.             Q.high -= 1;
161.             // 更新被除數的最高位下標
162.             while(R.data[R.high] == 0) {
163.                 R.high--;
164.                 t--;
165.             }
166.             t+=1;               // 下一位被除數
167.         }
169.         Q.high = high - dest.high + 1; 
170.         while(Q.data[Q.high] == 0) Q.high -= 1;
171.         R.high = high;
172.         while(R.data[R.high] == 0) R.high -= 1;
173.     }
175.     // 大數模冪演算法，很簡單的自然演算法，即將指數分解為二進位，換句
176.     // 更簡單的話來說，就是不斷地找平方模冪，而不是全部乘方後再
177.     // 做一次最終的模運算
178.     // a.power_mod(p, n)計算a^p mod n
179.     SuperNumber power_mod(int power, SuperNumber n) const {
180.         SuperNumber c("1"), t(*this);
182.         while(power) {
183.             if(power % 2) {
184.                 c = c * t % n;
185.                 power -= 1;
186.             } else {
187.                 t = t * t % n;
188.                 power /= 2;
189.             }
190.         }
192.         return c%n;
193.     }
195.     bool operator>=(const SuperNumber& s) const {
196.         if(high == s.high) {
197.             int k = high;
198.             while(data[k] == s.data[k] && k >= 1)k--;
199.             if(k < 1) return true; // equal
200.             return data[k] > s.data[k];
201.         } else if(high > s.high) return true;
202.         return false;        
203.     }
205.     bool operator<(const SuperNumber& s) const {
206.         return !(*this >= s);
207.     }
209.     // 從十進位表示的最高位開始數起，數到第i位，從第i位開始截取連續
210.     // 的c位元字出來組成一個新的數。例如：資料是12345678925698，則
211.     // sub(3, 5)返回數字34567，如果數字不夠取，例如在34567上運行
212.     // sub(3, 5)，因為34567從高位元起第3個數字是5，剩下的數字是567，
213.     // 最多隻有三個，不夠取要求的5個，這個時候返回567，不報錯。
214.     SuperNumber sub(int i, int c) const {
215.         SuperNumber ret;
216.         assert(high >= i);   // 保證可截取
218.         i = high - i + 1;    // 從高位元起的第i個數位的下標
219.         if(i >= c) {
220.             ret.high = c;
221.             while(c >= 1) ret.data[c--] = data[i--];
222.         } else {
223.             ret.high = i;
224.             while(i >= 1) {
225.                 ret.data[i] = data[i];
226.                 i--;
227.             }
228.         } 
229.         // 過濾前置0
230.         while(ret.data[ret.high] == 0) ret.high--;
231.         return ret;
232.     }
234.     // I/O
235.     friend istream& operator>>(istream& in, SuperNumber& s) {
236.         char t[256];
237.         in >> t;
238.         s.setToStr(t);
239.         return in;
240.     }
241.     friend ostream& operator<<(ostream& out, const SuperNumber& s) {
242.         return out << s.toString(10);
243.     }
244. private:
245.     enum {MAX_SIZE=256};        // 最大十進位位元
246.     // 須注意，使用data[0]儲存最高位所在下標是自己的一點小聰明，後來在
247.     // 調試的時候發現這是一個極大的錯誤，不過對於此題目來說可以應付
248.     char data[MAX_SIZE];        // 資料的內部表示，字串形式的十進位
249.                                 // 其中data[0]儲存最高位所在下標，data[1]
250.                                 // 儲存資料的最低位，也就是個位
251.     int high;
252. };

主函數：

1. int main(int argc, char **argv) {
2.     freopen("in.txt", "r", stdin);
4.     SuperNumberTest st;
6. //    st.run();
8.     // g和n都是超過2^127的素數。它們在DH演算法中公開
9.     SuperNumber g, n;
10.     int a, b;
11.     SuperNumber ka, kb, key;
13.     srand(time(0));
15.     cin >> g >> n;
17.     cout << "g = " << g << endl
18.          << "n = " << n << endl;
20.     cout << "/nThis is Alice:/n";
21.     a = rand();
22.     cout << "Alice get a random integer a = " << a << endl;
23.     cout << "Alice computer g^a mod n:/n";
24.     ka = g.power_mod(a, n);
25.     cout << "Alice compute out ka = " << ka << endl;
27.     cout << "/nThis is Bob:/n";
28.     b = rand();
29.     cout << "Bob get a random integer b = " << b << endl;
30.     cout << "Bob compute g^b mod n:/n";
31.     kb = g.power_mod(b, n);
32.     cout << "Bob compute out kb = " << kb << endl;
34.     cout << "/nAlice get kb from Bob, she compute out key is:/n";
35.     cout << kb.power_mod(a, n) << endl;
37.     cout << "/nBob get ka from Alice, he compute out key is:/n";
38.     cout << ka.power_mod(b, n) << endl;
41.     return 0;
42. }

運行結果：
g = 170141183460469231731687303715884105757
n = 170141183460469231731687303715884106309
 
This is Alice:
Alice get a random integer a = 20276
Alice computer g^a mod n:
Alice compute out ka = 102075421398841759242347870420481896337
 
This is Bob:
Bob get a random integer b = 28664
Bob compute g^b mod n:
Bob compute out kb = 62348451302684698452476840835428450852
 
Alice get kb from Bob, she compute out key is:
80402514625208456390620786920929643017
 
Bob get ka from Alice, he compute out key is:
80402514625208456390620786920929643017
在http://oldpiewiki.yoonkn.com/cgi-bin/moin.cgi/DiffieHellmanKeyExchange 上面給出了使用GNU gmp庫實現的這個演算法的示範代碼