已知兩整數 M < N
使得 Φ(M,N)
成立。上次我們介紹了 Netty van Gasteren 和 Wim Feijen 的二元搜尋法,能找到某個 l
, M ≤ l < N
,使得 Φ(l,l+1)
成立。給定一個排序好的陣列 a[0..N)
(其元素為 a[0]
, a[1]
... a[N-1]
), 0 ≤ N
。如何用 van Gasteren 與 Feijen 的方法判斷其中是否含有某個關鍵值 K
呢?
第一個念頭是用這個 Φ
:
Φ(i,j) = a[i] ≤ K < a[j]
並令 M = 0
。演算法結束後,Φ(l,l+1) = a[l] ≤ K < a[l+1]
一定成立,只要再看看 a[l]
是否等於 K
就可以了。但問題來了:van Gasteren-Feijen 演算法正確的兩個先決條件之一是 M < N
-- van Gasteren-Feijen 演算法的迴圈設計便假設 [M,N)
這個區段不是空的,因此我們無法處理空陣列。其次是 Φ(0,N)
得要成立,但 a[0] ≤ K < a[N]
並不一定成立呢。
我們可以把上述的例外都分開測試。但 van Gasteren 與 Feijen 建議的作法是在陣列頭尾各補一個想像元素:a[-1]
小於任何數,a[N]
大於任何數。一個等價的說法是用這個Φ
:
Φ(i,j) = (i = -1 ∨ a[i] ≤ K) ∧ (K < a[j] ∨ j = N)
最初令 l, r := -1, N
, 那麼初始條件就滿足了。程式如下:
{ 0 ≤ N ∧ Φ(-1,N) }
l, r := -1, N
{ Inv: -1 ≤ l < r ≤ N ∧ Φ(l,r), Bound: r - l }
; do l+1 ≠ r →
{ l + 2 < r }
m := (l + r) / 2
; if a[m] ≤ K → l := m
[] K < a[m] → r := m
fi
od
{ -1 ≤ l < N ∧ Φ(l,l+1) }
; if l > -1 → found := a[l] = K
[] l = -1 → found := false
fi
如果在陣列頭尾補元素讓你覺得很不妥,別擔心。迴圈恆式保證了 -1 < m < N
,因此從頭到尾,a[-1]
和 a[N]
兩個位置都沒被讀過 -- 陣列 a
並不用真正有什麼改變,兩個想像元素只是為了證明演算法正確而假想出來的。這也讓我們可以處理空陣列了。我覺得這是 van Gasteren-Feijen 演算法漂亮之處。
Bentley 的程式
Jon Bentley 在 Programming Pearls 一書中給的二元搜尋法程式如果翻譯成 guarded command language, 並把陣列索引從 [1..N]
改成 [0..N-1]
, 大約是這樣:
l, r := 0, N-1 ; do l ≤ r → m := (l + r) / 2 ; if a[m] < K → l := m + 1 [] a[m] = K → found := true; break [] K < a[m] → r := m - 1 fi od ; found := false
Bentley 書中也有以口語描述的正確性證明。我原想在課堂上講解這個程式,畢竟這個版本流傳較廣 -- 許多函式庫裡的二元搜尋法程式就是這麼寫的。但困難之一是若要把 K < a[m]
和 r := m - 1
牽上關係,似乎非得提早用上陣列已排序好的性質。這個演算法便比較難在更廣泛的脈絡中解釋了。當然,另一個困難是我不知如何在 Hoare logic 中簡單地解釋 break
。
乍看之下我以為 Bentley 演算法的搜尋範圍縮小得比 van Gasteren-Feijen 快:l
和 r
分別設為 m+1
與 m-1
,並不是 m
. 細看之後又發現並不盡然。變數 l
可設為 m+1
, 因為 m
的位置已由另一個比較 a[m] = K
處理;變數 r
設為 m-1
,則可能僅因為 Bentley 將陣列區段用 a[l..r-1]
表示,而 van Gasteren-Feijen 演算法把同一個區段表達為 a[l..r)
.
Bentley 與 van Gasteren-Feijen 演算法解的問題並不完全相同。當 K
在陣列中出現一次以上,Bentley 演算法不限定傳回哪個,van Gasteren-Feijen 則必定傳回索引最大的那個。當 K
不在陣列中時,van Gasteren-Feijen 演算法似乎較快,因為兩者都得跑完,而 van Gasteren-Feijen 的迴圈中只有一個比較(最後一個比較可省去)。Bentley 演算法若提早發現 K
可隨時跳出迴圈,代價是迴圈中多了一個比較。當陣列中確實找得到 K
, 這樣的交換是否值得呢?Timothy J. Rolfe 的實驗似乎認為只用一個比較的演算法平均上仍快一些。
習題
Van Gasteren 與 Feijen 建議的幾個習題中包括這個:假設陣列 a[0..N)
是一串嚴格遞增的數字緊接著一串嚴格遞減的數字,試用二元搜尋法找到陣列中的最大值。本問題中我覺得合理的假設是遞增和遞減數列均可能為空的,但 a
至少有一個元素:
0 < N ⋀
(∃ M: 0 ≤ M < N :
(∀ i,j : 0 ≤ i < j ≤ M : a[i] < a[j]) ⋀
(∀ i,j : M ≤ i < j < N : a[i] > a[j]))
這剛好是「最大密度區段」問題的演算法之一需要的一個副程式。我記得當時也是先隨便寫寫,幾次寫不對才發現真是難。現在總算是知道一點理論背景,覺得安心多了。你要試試看嗎?
中間索引怎麼取?
關於中間值 m := (l+r)/2
, 後來另有些其他故事。Google 的 Joshua Bloch 發現當陣列夠大時,l
和 r
相加可能造成溢位。他可不是故意挑毛病,這個 "bug" 是 Sun 抓到的。他建議該這麼寫:
int m = l + ((r - l) / 2); /* for Java/C/C++ */
int m = (l + r) >>> 1; /* for Java */
m = ((unsigned int)l + (unsigned int)r)) >> 1; /* for C/C++ */
Bloch 的發現自然引起了不少討論。有人認為這是整數型別的問題,而不是二元搜尋法的 bug. 有人質疑,討論 int
無法索引的陣列是否有意義?接下來可能一路談到組合語言定址法上頭。有人說,要造成索引溢位,陣列得有 4GB 大,你這輩子在 4GB 的陣列中作二元搜尋過嗎?有人便答,人家可是 Google 來的。Google 的陣列比常人的大上幾倍也不是不可能的唷。
如果往抽象的極端走,Roland Backhouse 在 Program Construction: Calculating Implementations from Specifications 書中則建議應該用 m := (l + r - 1) / 2
-- 如此一來不管該語言的整數除法到底是無條件消去、四捨五入、或是無條件進入,算出來的 m
值都會是 ⌊(l + r)/2⌋
.
參考資料
- Roland Backhouse. Program Construction: Calculating Implementations from Specifications. John Wiley and Sons, 2003.
- Jon Bentley. Programming Pearls, Second Edition. Addison-Wesley, 2000.
- Joshua Bloch. Extra, Extra - Read All About It: Nearly All Binary Searches and Mergesorts are Broken. Google Research Blog, June 02, 2006.
- Netty van Gasteren and Wim Feijen. The binary search revisited. AvG127/WF214, 1995.
- Timothy J. Rolfe. Analytic derivation of comparisons in binary search. SIGNUM Newsletter, Vol. 32, No. 4 (Oct 1997), pp. 15-19.
補充一些資料,依照The Art of Computer Programming第6.2.1小節,Single Comparison的二元搜尋法最早是在H. Bottenbruch的論文 “Structure and Use of ALGOL 60,” Journal of the ACM, Vol. 9, No. 2, p.161-221, April 1962中提到的,不過他的不同點是每次都取m=ceil( (l+u)/2 )。按照此節習題23的內容,Knuth分析了Bottenbruch的二元搜尋法的效率為(17.5 lg n + 17)u,而他自己設計的二元搜尋法的效率為(18 lg n – 16)u,u為單位執行時間。