這學期在學離散數學的代數與集合論,剛好上到了 Stirling 數這一章,老師便叫我們寫個程式來計算 Stirling 數並撰寫實驗報告,要用遞迴和非遞迴方式實現並比較兩者差異,計算公式本來就是遞推,所以用遞迴很好實現,主要是非遞迴的實現思路較為難想,當時其實已經把這篇文章的原型寫好了,最後的時間比較我想要用圖表來呈現,不過當時事物繁忙沒時間研究怎麼繪製,擱到了寒假才繼續完成。

概念

Stirling 數是用來計算一個集合中可以分組的組合數量,公式如下,n 是集合中元素數量,k 是組合數,利用遞迴關係可以求出 Stirling 數。

$$S(n, k)=S(n-1, k-1)+k \cdot S(n-1, k)$$

表格遞推法

假設欄數代表公式中的 n,列則是 k,已知當 $n=1$ 和 $n=k$ 時的組合數為 1(n 個元素分 1 組的組合數為 1,且 n 個元素分 n 組的組合數為 1),由此可遞推更高階的 Stirling 數。

根據上面的遞迴關係式,要求的那格等於欄數乘上一列格子的值再加上左上方格子的數字,假設我要求第 6 列第 3 欄的值,已知要求 $S(6, 3)$,左上角($S(5, 2)$)數字是 15,在加上上面那格 $S(6, 2)$,就可以得到 $15+3\times25=90$。

1 2 3 4 5 6
1 1
2 1 1
3 1 3 1
4 1 7 6 1
5 1 15 25 10 1
6 1 31 90 65 15 1

程式撰寫

假設

我認為遞迴是最簡單快速的方法,因為不用考慮其他複雜邏輯,而非遞迴還需要一個時間複雜度將所有堆疊彈出求和,所以我先假設遞迴的時間應該會比用非遞迴的時間時間短,不過根據我改良的方法時間複雜度應該會大幅下降。

輸入

輸入兩個數字,第一個數字是欄數 $n$,第二個數字是列數 $k$。

使用遞迴

用遞迴來寫是最直覺簡單的方式了,程式碼也較為簡單,遞迴中止的條件設 $k=1$ 和 $n=k$ 時返回 1,核心程式碼如下。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
int s(n, k)
{
    if(n == 1 || n == k)
    {
        return 1;
    }
    else
    {
        return s(n-1, k-1)+k*s(n-1, k);
    }
}

非遞迴

一般人可能直覺想到用陣列來畫表格推算,但是這樣時間、空間複雜度會很高,經過我一夜的思考,想到用堆疊(stack)來存儲遞迴的過程,實驗多次總算寫出成品了,先寫個結構來記錄 Stirling 數的參數,n, k 前面解釋過了,而 times 的目的是用來記錄這個 Stirling 數前面的乘數,還記得前面 Stirling 數的遞迴關係式嗎?

$$S(n, k)=S(n-1, k-1)+k \cdot S(n-1, k)$$

times 就是代表 $S(n-1, k)$ 前面的 $k$,預設值為 1。

1
2
3
4
5
6
struct Stirling
{
    int n;
    int k;
    int times;
};

定義兩個堆疊,一個存結構,一個存算出來的值。

1
2
stack<Stirling> s;
stack<int> ans;

先輸入 n, k,並寫入到 Stirling 的結構裡推入推疊中。

1
2
3
4
int n, k;
scanf("%d %d", &n, &k);
Stirling stirling = {n, k, 1};
s.push(stirling);

接下來就是核心部份了,會先從 s 取堆疊頂端,看 n, k 是否符合條件(n==1 || n==k),如果條件為真就在 ans 推入 times 的值並彈出 s,否則將 s 彈出,聲明兩個 Stirling 型態的結構 new1new2,分別代表公式中的 $S(n-1, k-1)$ 和 $S(n-1, k)$ 並推入堆疊中,這步驟是為了將原本的 $S(n, k)$ 拆分成子問題,當條件符合中就會彈出 s,直到 s 為空迴圈終止,而拆分出來的子問題解答會存在 ans 堆疊中。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
while (!s.empty())
{
    Stirling top, new1, new2;
    top = s.top();
    if (top.n == top.k || top.n == 1)
    {
        ans.push(top.times);
        s.pop();
    }
    else
    {
        new1 = {top.n - 1, top.k - 1, top.times * 1}, new2 = {top.n - 1, top.k, top.times * top.k};
        s.pop();
        s.push(new1);
        s.push(new2);
    }
}

最後將 ans 中的值一個一個彈出求和就是答案了。

1
2
3
4
5
6
7
int sum = 0;
while (!ans.empty())
{
    int a = ans.top();
    sum += a;
    ans.pop();
}

時間比較

我在計算完成後將計算時間寫入到檔案裡,並用 Python + Matplotlib 繪製三維散佈圖表,可視化兩種實現方法的時間差異,初學 Matplotlib 製圖,下面圖片都是直接在 GUI 界面儲存的,如果有更好的方法歡迎告訴我。

下面是 n 從 1~10 的時間比較,橘色三角形是非遞迴的實現時間,可以明顯看出非遞迴在 $k$ 值愈高的時候所花費的計算時間比遞迴多愈多。

n 從 1~10

接下來把 $n$ 的範圍加大到 25,其實我本來是想跑到 100 的,但是我發現這樣計算量太大,所以只算到 $n = 25$,這樣差不多有兩千多筆資料,可以看見非遞迴所用時間隨著 k 的變大呈指數型成長。

n 從 1~25

結論

即使我的非遞迴實現採用堆疊去減少時間複雜度和空間複雜度,但還是比不過用遞迴實現的速度,這種遞推公式貌似遞迴做是最快的。