本文列舉了對於 演算法 : 第4版 / (美) 塞奇威客 (Sedgewick, R.) , (美) 韋恩 (Wayne, K.) 著 ; 謝路雲譯. -- 北京 : 人民郵電出版社, 2012.10 (2021.5重印)（以下簡稱原書或書）中的練習題 1.1.27 的三種解法（C++ 實現），並對包含原書題中的遞迴方法在內的四種解法的執行時間進行了統計和對比。
◆ 要求
原書中的練習題 1.1.27 要求對如下二項分佈遞迴過程中的值儲存在陣列中，
b(n,k,p) = 1.0 ( n == 0 && k == 0 )
b(n,k,p) = 0.0 ( n < 0 || k < 0 )
b(n,k,p) = (1.0-p) * b(n-1,k,p) + p * b(n-1,k-1,p)
◆ 解一
依然採用遞迴的方式，但使用二維陣列儲存中間結果。
如下程式碼所示，
static long double binomial1(int N, int K, long double p) // #1
{
long double x;
long double** b = new long double*[N+1]; // #2
long double* data = new long double[(N+1)*(K+1)];
...
x = binomial1_impl(b, N, K, p); // #3
...
return x;
}
static long double binomial1_impl(long double** b, int N, int K, long double p)
{
if (N == 0 && K == 0) return 1.0L;
if (N < 0 || K < 0) return 0.0L;
if (b[N][K] == -1) {
... // #4
b[N][K] = (1.0L-p) * binomial1_impl(b, N-1, K, p) + p * binomial1_impl(b, N-1, K-1, p);
}
return b[N][K];
}
保持對外的介面不變（#1），建立一個二維陣列 b[0..N][0..K] 儲存中間計算結果（#2），並將其傳給演算法實現（#3）。演算法雖然還是用遞迴呼叫（#4），但由於中間結果儲存在全域性的二維陣列中，不用頻繁地壓棧和彈棧去獲取中間資料。此解法網路上也見於 [github] reneargento/algorithms-sedgewick-wayne 和 [github] aistrate/AlgorithmsSedgewick 。
◆ 解二
使用二維陣列保持中間結果，但同時將遞迴改進為遞推。若以橫向為 x 軸，縱向為 y 軸，左上角為座標原點，則座標軸上的 (x,y) 點則代表二維陣列的 b[y][x] 單元。
以 N = K = 4 為例，
0 1 2 3 4
0 * * * * * <-- * 代表待計算的單元
1 * * * * *
2 * * * * *
3 * * * * *
4 * * * * ? <-- 最終計算結果的單元 (4,4)
仔細考察遞迴關係式的特點，b(-1,*,p) = 0.0, b(*,-1,p) = 0.0。由
b(0,1,p) = (1.0-p) * b(-1,1,p) + p * b(-1,0,p)
= (1.0-p) * 0.0 + p * 0.0
= 0.0
b(0,2,p) = (1.0-p) * b(-1,1,p) + p * b(-1,1,p)
= (1.0-p) * 0.0 + p * 0.0
= 0.0
...
可推論出，二維陣列中的第 0 行中的所有單元（不含b[0][0]）均為 0.0；由
b(1,0,p) = (1.0-p) * b(0,0,p) + p * b(0,-1,p)
= (1.0-p) * 1.0 + p * 0.0
= 1.0-p
b(2,0,p) = (1.0-p) * b(1,0,p) + p * b(1,-1,p)
= (1.0-p) * (1.0-p) + p * 0.0
= (1.0-p)^2
...
可推論出，二維陣列中的第 0 列的單元為 (1.0-p)^y。
因為每個單元 b[n][k] 結果（n 代表行號，k 代表列號），依賴於 b[n-1][k-1] 和 b[n-1][k] 的結果。為了減少計算量，遞推過程可僅用到二維陣列的部分單元。筆者設定一個 G 點，將待計算單元的區域劃分為 '#' 和 '*' 兩部分，G 點在 '#' 區域中。分為以下三種情況，
第一種情況，N < K：（如 N = 4, K = 6）
0 1 2 3 4 5 6
0 - - G # # # # <-- G 點所在單元為 0.0
1 - - - * * * * <-- '-' 代表不用計算的單元
2 - - - - * * *
3 - - - - - * *
4 - - - - - - ? <-- 最終結果的儲存單元
G 點為 b(0,K-N)。按照遞推關係式容易推匯出，'#' 和 '*' 區域均為 0.0，所以最終結果即 0.0。
第二種情況，N = k：（如 N = 6, K = 6）
0 1 2 3 4 5 6
0 G # # # # # # <-- G 點所在單元為 1.0
1 - * * * * * *
2 - - * * * * *
3 - - - * * * *
4 - - - - * * *
5 - - - - - * *
6 - - - - - - ?
G 點為 b(0,0)。按照遞推關係式容易推匯出，陣列中 n = k 的單元為 p^n。所以最終結果即 p^N。
第三種情況，N > K：（如 N = 6, K = 4）
0 1 2 3 4
0 # # # # #
1 # # # # #
2 G # # # # <-- G 點所在單元為 (1.0-p)^2
3 - * * * *
4 - - * * *
5 - - - * *
6 - - - - ?
G 點為 b(N-K,0)。可先計算 '#' 區域中的單元，再計算 '*' 區域中的單元，得出最終結果。處理'#'區域時，為避免大量的陣列下標越界判斷，可以考慮先計算 0 行和 0 列的所有單元。
如下程式碼所示，
static long double binomial2(int N, int K, long double p)
{
long double x;
if (N < K) { // #1
x = 0.0L;
} else if (N == K) { // #2
x = powl(p, N);
} else { // #3
...
b[0][0] = 1.0L;
// process '#' area // #4
// calcuate [1..N-K][0]
for (i = 1; i <= N-K; ++i)
b[i][0] = powl(1.0L-p, i);
// calcuate [0][1..K]
for (j = 1; j <= K; ++j)
b[0][j] = 0.0L;
// calcuate [1..N-K][1..K]
for (i = 1; i <= N-K; ++i)
for (j = 1; j <= K; ++j)
b[i][j] = (1.0L-p) * b[i-1][j] + p * b[i-1][j-1];
// process '*' area // #5
for (i = N-K+1; i <= N; ++i)
for (j = i-(N-K); j <= K; ++j)
b[i][j] = (1.0L-p) * b[i-1][j] + p * b[i-1][j-1];
x = b[N][K]; // #6
...
}
return x;
}
三條分支（#1、#2、#3）分別對應前述三種情況。在第三種情況下，再先處理 '#' 區域（#4），然後採用遞推求值的方式處理 '*' 區域（#5），最後得到結果（#6）。
◆ 解三
此方法是從遞推解法中引申出來了。進一步探究這個此二項分佈的遞迴式，以 N = 4 且 K = 4 為例，
0 1 2 3 4
0 1.0
1 1.0-p p
2 (1.0-p)^2 2*(1.0-p)*p p^2
3 (1.0-p)^3 3*[(1.0-p)^2]*p 3*(1-p)*(p^2) p^3
4 (1.0-p)^4 4*[(1.0-p)^3]*p 6*[(1.0-p)^2]*(p^2) 4*(1.0-p)*(p^3) p^4
可以發現，從第 0 行到第 N 行的非零單元即“楊輝三角形”，第 n 行中的非零單元之和構成 [(1.0-p) + p]^k 的展開式。因此，解二中的第三種情況，可結合利用通項公式 C(N,K)*[(1.0-p)^(N-K)]*(p^K) 來解決。
如下程式碼所示，
static long double binomial3(int N, int K, long double p)
{
long double x;
if ...
} else {
x = combination(N, K) * powl(1.0L-p, N-K) * powl(p, K);
}
return x;
}
◆ 測試
編譯並執行程式，
$ g++ -std=c++11 -o 27.out 27.cpp
$ ./27.out 10 5 0.25
為易於顯示兩者之間的差異，筆者選擇了硬體配置偏低的測試環境。
- 硬體配置：Raspberry Pi 3 Model B
- Quad Core 1.2GHz 64bit
- 1G RAM
- 16G MicroSD
- 100 Base Ethernet
- 軟體配置：Raspbian Stretch
- g++ (Raspbian 6.3.0-18+rpi1+deb9u1) 6.3.0 20170516
測試並記錄了 (N, K, p) 為 (10, 5, 0.25)， (20, 10, 0.25)， (40, 20, 0.25)， (80, 40, 0.25)， (100, 50, 0.25) 的情況下，原遞迴、解一、解二、解三執行時所消耗的時間。
結果如下圖所示，
對比可以看出，不同的解法在執行時間上的差異隨著計算量的增加而逐步擴大。
◆ 最後
完整示例程式碼和測試結果，請參考 [gitee] cnblogs/17328989 。
寫作過程中，筆者參考了 [github] reneargento/algorithms-sedgewick-wayne 和 [github] aistrate/AlgorithmsSedgewick 的實現方式。致 reneargento 和 aistrate 兩位。