2011年4月10日 星期日

TAOCP閱讀心得: Euclid's algorithm

讀到 TAOCP 1.1 尾聲, Knuth 用數學的形式表示 Euclid's algorithm, 也就是用輾轉相除法求最大公因數。覺得書上的表示方式很妙。原本的英文描述如下:
E1. [Find remainder.] Divide m by n and let r be the remainder. (We will have 0 <= r < n.)

E2. [Is it zero?] If r = 0, the algorithm terminates; n is the answer.

E3. [Reduce.] Set m <-n, n <- r, and go back to step E1.
數學定義摘錄函式定義部份如下:
f((m, n) = (m, n, 0, 1); f((n)) = (n);

f((m, n, r, 1)) = (m, n, remainder of m divided by n, 2);

f((m, n, r, 2)) = (n) if r = 0,      (m, n, r, 3) otherwise;

f((m, n, p, 3)) = (n, p, p, 1)
其中 m, n, p 為 >0 的整數, r 為非負數整數。

在上廁所的時候忽然想通 (廁所真是靈感湧現的好地方), 這和以前讀 Haskell 時看到的東西很像, 將狀態編碼到輸出入裡, 於是可以使用函數定義來達到迴圈的效果。以上例來說, 就是隱藏 while 迴圈, 將 while 內的每個操作拆成連續的數個算式, 用輸入的最後一個數字 (即 1, 2, 3) 來表示和銜接算式的順序。若需要用到 for-loop 的 i 的話 (其實大多數情況用不到它), 可用同樣的方式, 將 i 表現在函數的輸出入上。

於是回頭翻一下 Haskell 的語法, 用生澀的寫法練習表示上面的數學式子, 結果如下:
f(m, n, x, 1) = (m, n, mod m n, 2)
f(m, n, r, 2) | r == 0 = (n, n, n, 4)
f(m, n, r, 2) | r /= 0 = (m, n, r, 3)
f(m, n, p, 3) = (n, p, p, 1)
f(a, b, c, 4) = (a, b, c, 4)

印象中 Haskell 的函數輸出入必須為同一型別, 所以我改寫了一下, 讓程式符合 Haskell 的語法。應該是有其它簡單的寫法, 可以完全表示上面的數學定義, 比方輸入用 list, 自己切參數看數量再決定行為。以後有動力重頭學 Haskell 的話, 再來改寫。

下面是執行 gcd(18, 15) 的過程:
*Main> f(18, 15, 0, 1)
(18,15,3,2)
*Main> f(18, 15, 3, 2)
(18,15,3,3)
*Main> f(18, 15, 3, 3)
(15,3,3,1)
*Main> f(15, 3, 3, 1)
(15,3,0,2)
*Main> f(15, 3, 0, 2)
(3,3,3,4)

一開始我很納悶為什麼 Knuth 要多一步 f((m, n, p, 3)) = (n, p, p 1), 而不是在 r != 0 時就表示成 (n, p, p, 1)。對照原本的數學描述, 才明白這樣才能精確的表示執行成本, 讓每個步驟都是明確且有效的, 這樣寫更符合 algorithm 的定義。

那麼, 以上這串體悟對現實開發軟體有何幫助呢? 我覺得這是一個精簡表示數學和演算法的好例子。過去, 我從學習 imperative programming 到學習 OOP, 還有後來稍微學習 functional programming (相關介紹), 發覺兩者有很大的差異。有些在 imperative programming / OOP 裡很難做的事, 在 FP 裡卻不是什麼困難的事。基本論點不同, 表現想法的方式不同, 帶來的優缺點也不同。一直想找機會從頭練習 FP 的思維。附帶一提, 最近 CMU 決定要用 ML 教大一 FP。FP 似乎愈來愈紅了。

2011-04-10 更新

經曾出國深造的 command 說明, 查了一下定義 data type 的方法, 改寫 Haskell 程式如下:
data I = V1 (Integer)
         | V2 (Integer, Integer)
         | V4 (Integer, Integer, Integer, Integer) deriving Show

g :: I -> I
g(V1 (x)) = V1 (x)
g(V2 (m, n)) = V4 (m, n, 0, 1)
g(V4 (m, n, x, 1)) = V4 (m, n, mod m n, 2)
g(V4 (m, n, r, 2)) | r == 0 = V1 (n)
g(V4 (m, n, r, 2)) | r /= 0 = V4 (m, n, r, 3)
g(V4 (m, n, p, 3)) = V4 (n, p, p, 1)
感覺上好像沒有比較清楚, 或是我沒抓到或不習慣 Haskell 的風格吧。