---
tags: Linux Kernel Internals, 作業系統
---
# 現代硬體架構上的演算法: 以 Binary Search 為例
## 序
時間複雜度是我們在評估演算法效益時經常參考的指標。然而，在真實的應用中，演算法是否能夠「在限定的條件下」，「以最短的時間完成定量的問題」也就是「具有最大的吞吐量(throughput)」才是對我們而言可以決定其優劣的關鍵。而時間複雜度僅僅只是影響該評判標準的一個因子。時間複雜度低並不意味在**任何場景**都是花費最短的時間。以下圖與附文為例，即便我們知道 insertion sort 的時間複雜度($O(n^2)$)高於 merge sort($O(n\lg{n})$)，當陣列很小時可以看到在實務場景上反而 insertion sort 是能在更短時間下完成的。
事實上，在現代的計算機中，舉凡 cache, pipeline, branch prediction, SIMD 等等硬體架構的特性也會對演算法實際運行的時間產生影響。
> [PA4: Sorting Improvements](https://w3.cs.jmu.edu/lam2mo/cs240_2014_08/pa04-sorting.html)
從另一個角度來看，除了 CPU 運算上的成本，演算法真正的執行時間還要加上從記憶體中取得資料的時間。因此，運行時間的瓶頸可能受 CPU 或記憶體限制。如此來看，軟體的開發絕對不能僅僅著眼於程式邏輯的優化，也需要考量硬體特性帶來的效能差異。本文將以 [Optimizing Binary Search](https://youtu.be/1RIPMQQRBWk) 為基礎進行導讀，並延伸探討 Binary Search 在實作上應該如何考慮現代計算機架構產生的影響。
## 定義問題與應用場景
:::warning
讀前須知:
1. 本文將原文的問題定義與實作做了稍許調整，讓所有比較方法都可以基於同樣的函示原型(prototype)實作。藉此我們可以在更為公平的角度來評斷各個演算法的優劣
2. 此外，為方便非 C++ 語言的用戶(我)理解，這裡重新以 C 進行實作。但由於圖表部份目前仍是引用原文。因此這裡要注意縱使 porting 後演算法邏輯相同，但因為編譯器行為差異仍可能在 C 和 C++ 上導出相異的時間結果
:::
```cpp
int *prepare(int *src_arr, int n);
int lower_bound(int *arr, int n, int val);
void clean(int *arr);
```
我們希望實作三個函式:
* `prepare`:
* 在 benchmark 被啟動的最初執行一次
* 輸入長度為 `n` 的**已排序**陣列 `src_arr` 後，演算法可以根據需求決定是否進行額外調整，或者配置額外空間，並將前處理結果輸出到該額外空間中
* 將以回傳的指標將做為後續函式的輸入
* `lower_bound`:
* 功能類似於 c++ 中的 [`lower_bound`](https://cplusplus.com/reference/algorithm/lower_bound/)。輸入的 `arr` 會是 `prepare` 的回傳指標，`n` 則與在 `prepare` 時輸入的 `n` 相同(因此不一定是輸入的 `arr` 之大小)
* `lower_bound` 會隨機的 query 值 `val`，實作上需要找到之前 `prepare` 時給定的 `src_arr` 中最小的 `src_arr[idx]` 以符合 `src_arr[idx] >= val` 並返回之，若符合的答案不存在則回傳 -1
* 注意到 `lower_bound` 是 benchmark 測試的主體，會在 loop 中以不同的隨機 `val` 多次執行。
* `clean`:
* benchmark 結束後被呼叫，可用來回收配置的資源
為了鎖定貼近實務的應用場景，並更為直接的展示演算法與現代硬體的關係，加上額外的限制和說明:
* 假設 int 是 32 bits
* 假設 cache line 大小為 64 bytes
* 假設使用 x86-64 架構的 CPU
* 假設 `prepare` 中若發生記憶體配置，必定成功，且 benchmark 最後會釋放之
* `prepare` 的實作時間複雜度要求為 O(n)，則因為實務上既然已經允許費時取得已排序的陣列(時間複雜度可視為O(nlogn))，且 `prepare` 相較於 `lower_bound` 的執行次數不多，因此 benchmark 中允許可以忽略其執行的時間
* 僅考慮 comparison-based search 的方法，換句話說，不能使用 lookup table 等方式
* 評判的標準是演算法的吞吐量(throughput)
* 以下的測試數據是在 2GHz 時脈頻率的 CPU 平台上獲得(參照原文的 [AMD Zen2](https://www.7-cpu.com/cpu/Zen2.html))
那麼最為直覺的實作方式是搭配 binary search:
```cpp
int *prepare(int *src_arr, int n)
{
// do nothing, we don't need any preparation in this algorithm
return NULL;
}
int lower_bound(int *arr, int n, int val)
{
if (arr[n - 1] < val)
return -1;
int l = 0, r = n - 1;
while (l < r) {
int t = (l + r) / 2;
if (arr[t] >= val)
r = t;
else
l = t + 1;
}
return arr[l];
}
```
由於沒有額外配置記憶體，`clean` 在此可不必實作。
測試這個演算法的執行時間和陣列的大小的關係(虛線從左至右是 L1 - L3 cache 的大小)。我們可以看到以下的數據。
## Branch prediction 對 Binary Search 的影響
### Issue
首先，讓我們關注 < 32K 的部份。這部份的資料量由於是能夠完全被載入到 cache 上，因此可以先忽略 cache 的影響不計。但如果實際去計算 while loop 中每個 iteration 所需要的 CPU cycle，參考被圈起來的座標點:
```
100 ns / 12 iteration * 2 cycle/ns ~= 17 cycle/iteration
```
:::info
* 2 cycle/ns 是因為測試的 CPU 是 2GHz
* 12 iterations 是因為每個 iteration binary serach 都只會拆出其中一半
:::
大約 17 個 CPU cycle 可以完成。這比起實際從 assembly code 計算出來的 cycle 數還要更多一些(參照下方與 [Optimizing Binary Search](https://youtu.be/1RIPMQQRBWk) 影片說明)。如果用 [perf](https://en.wikipedia.org/wiki/Perf_(Linux)) 來觀察，會發現大部分的時間花費在條件分支的 jump instruction，而非 cmp instruction。這其中的原因為何呢? 這邊我們就得從 [CPU pipeline](https://en.wikipedia.org/wiki/Instruction_pipelining) 的原理說起。
```
│35: mov %rax,%rdx
0.52 │ sar %rdx
0.33 │ lea (%rsi,%rdx,4),%rcx
4.30 │ cmp (%rcx),%edi
65.39 │ ↓ jle b0
0.07 │ sub %rdx,%rax
9.32 │ lea 0x4(%rcx),%rsi
0.06 │ dec %rax
1.37 │ test %rax,%rax
1.11 │ ↑ jg 35
```
先簡單為各位複習一下計算機架構 :laughing: 一道 CPU 指令(instruction)的執行通常可以分為多個層級，而最基礎的四階段可以區分為
1. Fetch: CPU 從記憶體中取得下一道指令
2. Decode: 解析指令，例如取得 opcode 決定接下來要執行的邏輯、決定運算的 operand 等
3. Execute: 根據 decode 結果執行指令
4. Write back: 將執行結果寫回對應的記憶體位置(如果需要)
如下圖，CPU 的每個層級可以獨立運作，因此在同一個 clock cycle 中，指令與指令間可以在不同層級平行的運行。藉由這種方式可以讓 CPU 每單位時間能夠處理的指令數量提升。
然而並非所有的指令都是獨立的。假設後面的指令依賴前面指令的輸出，則指令間就必須要存在延遲，這在術語上稱為 [Hazard](https://en.wikipedia.org/wiki/Hazard_(computer_architecture))。
一種案例是: 在程式具有條件分支(branch) 的狀況下，理論而言我們需要等待計算出該條件的結果為 true 或 false，才能決定下一道指令，但為了擅用 pipeline，更理想的方式是 CPU 會透過 [branch prediction](https://en.wikipedia.org/wiki/Branch_predictor) 的機制預測下一道指令並執行，那麼假設預測成功，CPU 就可以最大化 pipeline 的優勢，但反之若預測失敗，就必須承擔已執行非預期指令的成本(flush pipeline, recover from mispredict)，這種狀況稱做 [control hazard](https://en.wikipedia.org/wiki/Hazard_(computer_architecture)#Control_hazards_(branch_hazards_or_instruction_hazards))。
另一種案例是，如果後面的指令所需的資料，需要依賴於前面指令的完成，那麼指令與指令間就沒辦法像上圖這樣排列，可能需要插入適當的延遲來確保資料的正確性，這種狀況則被稱為 [data hazard](https://en.wikipedia.org/wiki/Hazard_(computer_architecture)#Data_hazards)。
而在資料量可容納於 L1 cache 的情況下，影響 binary serach 效率的關鍵就在於 hazard。在具有條件分支的程式碼中，一方面 CPU 會對 branch 進行預測，因此會遇到 control hazard;另一方面，正確的程式邏輯依賴於條件的 true/false 被計算出來，因此也會有 data hazard 的發生。回顧前面 binary serach 中實做上會遇到的 branch:
* `while (l < r)`: 因為這個判斷式大部分時候都是 true，僅在最後要找到正確答案時為 false，得益於現代計算機 branch prediction 的能力，可以預期這部份的預測準確率是很高的
* `if (a[t] >= val)`: 在每個 iteration 中，此判斷式成立與否是沒有明顯的規則的，因此 CPU 就有很高的機會無法準確預測，進而導致額外的成本
### Branchless Binary Search
由於 branch 可能帶來額外的 CPU cycle，因此修改實作為 branchless 方式是有其必要的。
```cpp
int lower_bound(int *arr, int n, int val)
{
if (arr[n - 1] < val)
return -1;
int *base = arr, len = n;
while (len > 1) {
int half = len / 2;
if (base[half - 1] < val) {
base += half;
len = len - half;
} else {
len = half;
}
}
return *base;
}
```
相比於維護當前 iteration 所處理陣列之左和右 index，我們可以調整為維護該陣列的開頭和長度，如上面的實作。
```cpp
int lower_bound(int *arr, int n, int val)
{
if (arr[n - 1] < val)
return -1;
int *base = arr, len = n;
while (len > 1) {
int half = len / 2;
if (base[half - 1] < val)
base += half;
len -= half;
}
return *base;
}
```
觀察可以發現，下一個 iteration 的 len 原則上只會是 $\lceil \text{len/2} \rceil$ 或 $\lfloor \text{len/2} \rfloor$。由於在搜尋的時候，重要是排除不存在正確答案的範圍，因此我們可以一致都取 $\lceil \text{len/2} \rceil$ 就好。
```cpp
int lower_bound(int *arr, int n, int val)
{
if (arr[n - 1] < val)
return -1;
int *base = arr, len = n;
while (len > 1) {
int half = len / 2;
base += (base[half - 1] < val) * half; // will be replaced with a "cmov"
len -= half;
}
return *base;
}
```
如此一來，我們可以進一步改寫為上面的形式，透過將比較的 boolean 結果對應到 0 或 1，能夠排除掉對 if 的依賴。編譯器實際上會將這轉換為 [condition move(cmov)](https://mudongliang.github.io/x86/html/file_module_x86_id_34.html)。另外可以注意到的是，因為我們簡化每次計算 `len` 的方式，因此對於固定的輸入 `n`，while 的迭代次數變成是固定的，從而可以完全消除 branch prediction 錯誤的狀況。通過該實作可以得到如下圖的結果。
## CPU Prefetching 對 Binary Search 的影響
### Issue
回到前述章節最後的實驗結果，我們會發消除了 branch prediction 之後，雖然在 array 較小時的速度上有提昇，但當 array 變大之後，效率反而變差。
背後的原因是由於在有分支的版本中，CPU 擁有 [speculate](https://en.wikipedia.org/wiki/Speculative_execution) 機制，會預測其中一個分支的執行，並在確認預測正確之前，就已經開始獲取下個 iteration 的比較值。而在無分支的版本中，我們失去了 CPU 因為 branch 隱含的預測能力。由於實際上我們可以預測下一個比較值，我們可以引入軟體的方式主動 prefetch。
### Prefetching
修改程式碼如下:
```cpp
int lower_bound(int *arr, int n, int val)
{
if (arr[n - 1] < val)
return -1;
int *base = arr, len = n;
while (len > 1) {
int half = len / 2;
len -= half;
__builtin_prefetch(&base[len / 2 - 1]);
__builtin_prefetch(&base[half + len / 2 - 1]);
base += (base[half - 1] < val) * half;
}
return *base;
}
```
編譯器提供的 builtin [`__builtin_prefetch`](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html#index-_005f_005fbuiltin_005fprefetch) 提供軟體方式的 prefetch。則我們可以在預知下個 iteration 會存取的左 key 和右 key 的情況下，預先將資料預讀到 cache 中，以減少 cache miss 產生的延遲。
觀察結果，在 array size 變大時，branchless 實作方式的效能有明顯的提升。但由於 branch 版本還可以 prefetch 「下下個」、「下下下個」 等 iteration 的 key，因此 array 大到一定程度之後，branchless + prefetching 方式的效能仍然落後。儘管隨著預測越來越不可能是正確的，每個新的推測性讀取的有用程度也會呈指數級下降。
在 branchless 版本，我們也可以提前 prefetch 更後面的 iteration 會操作到的比較值，但我們需要的 prefetch 數量也呈指數增長。與之相對，我們將嘗試另一種不同的方法，以更睿智的方式來優化對 cache 的利用。
## Cache 特性對 Binary Search 的影響
### Issue
對於前面實驗的圖片，由於該圖是用 [Logarithmic scale](https://en.wikipedia.org/wiki/Logarithmic_scale) 繪製的，且我們知道實作的 Binary Serach 時間複雜度是 O(logn)，則理論上繪製出來要是 45 度角的直線。然而實際上可以看到曲線的成長是超出線性速度的。可以發現當陣列大小在超過 cache 大小的時候，讀取的速度有變慢的趨勢。顯然的，整個 array 是否可以完全載入到 cache 上將影響演算法的執行效率。
讓我們簡單回顧一下 cache 的硬體特性以及相關的性質:
* Cache hierarchy: 現代處理器中 cache 經常是多層結構的，從最快/最小到最慢/最大依序是 L1、L2 和 L3 cache。記憶體的訪問會按 L1 -> L2 -> L3 -> 主記憶體的順序
* Cache line: Cache 儲存資料的基本單位，例如比較常見的設計是 64 bytes。這意味著每當我們想載入一個 byte 到 cache 上，包含鄰近的整個 64 bytes 也會被一起被載入;反之，當我們更新 cache 上的一個 byte，同步回主記憶體時，整個 64 bytes 也都要被同步
* Eviction policy: 當新資料要被存在 cache 上，該挑選哪一個 cache line 以對舊資料進行 eviction 的規則
* Cache 常用的 policy 是 least recently used (LRU)，意即離至今最久以前被使用的 cache line 優先能被取代
* Bandwidth: 讀取或存儲數據到 cache 的速率
* Temporal locality: 是一種存取的模式。當一個記憶體位置被讀寫，很高機率在近期內會再次讀寫之
* Spacial locality: 是一種存取的模式。當一個記憶體位置被讀寫，很高機率鄰近的位置近期內會接著被讀寫
Cache 之所以強大，就是通過 temporal locality 和 spacial locality 來決定哪些資料該被保留/預載在 cache，以減少直接從主記憶體存取資料的需求。然而前述 binary serach 的實作之存取方式卻無法符合這些 locality。如下圖展示一個 binary search 的案例，可以觀察到這種存取方式的 locality:
* Spatial locality 僅在最後的 iteration 較佳，而最開始的搜尋總是跨越一大段記憶體位址
* Temporal locality: 在 benchmark 中，雖然我們會不斷的用不同的值 query，但可以預想每一次 query 最初的搜尋順序是很接近的。具例來說，上圖總是以 15 作為第一個 key 與 `val` 比較
我們可以利用實驗來證明上述的第二點對於速度的影響。
```cpp
int lower_bound(int *arr, int n, int val)
{
if (arr[n - 1] < val)
return -1;
int l = 0, r = n - 1;
while (l < r) {
int m = l + rand() % (r - l);
if (arr[m] >= val)
r = m;
else
l = m + 1;
}
return arr[l];
}
```
相較於擇取固定的中間點作為 key，改為隨機選擇在當下搜尋範圍的任一個。理論上(詳見[推導](https://en.algorithmica.org/hpc/data-structures/binary-search/#appendix-random-binary-search))，隨機的方式會讓比較的數量多出 30-40%，然而如下圖展示的結果，在大的 array 上整體運行時間卻可達到足足六倍的成長，且可以看到提升的幅度是指數級成長的，在 L2、L3 cache 的邊界有明顯的提升。這些證據都指出了性能下降不僅僅是因為額外的 `rand()` 呼叫引起的，更多原因來自任何獲取的 key 在這種實作方式中都不太可能被 cache 住。
Cache Associativity 同樣會產生不良影響，想像當 array 的大小是 $2^n$ 且 n 足夠大，則連續的兩個 key 可能映射到同一個 cache line，進而導致後者將前者 evict 的狀況。
而最主要的問題在於 temporal locality 的不良。如下圖所展示，越深色部分表示在 binary serach 下越頻繁存取的部分，可以看到在一般的排序陣列下，hot element 和 cool element 是分散排列的。這就導致 cache 無法利用 temporal locality。舉例來說，在讀取 15 位置時，cache 可能也會把後面的 16, 17, 18 預載起來，但後者顯然只有在相對少數的情況下才被存取。
### Eytzinger Binary Search
為了更有效的利用 locality，我們可以調整陣列的結構。這裡引入的方法是源於 [Michaël Eytzinger](https://en.wikipedia.org/wiki/Micha%C3%ABl_Eytzinger) 著名的族譜編號系統，使得可以將陣列調整為樹狀結構的表示方式。對於陣列(*1-based*)的第 $k$ 個元素，可以在 $2k$ 找到其左節點，並在 $2k + 1$ 找到右節點，利用此種方式，我們可以將已排序的陣列轉變成 binary tree，參考下圖:
從上圖中我們可以看出來為甚麼這種儲存方式對於 locality 更好: 由於 hot element 與 cold element 被更好的組織起來，因此在 `lower_bound` 的每個 iteration 中，cache 就更有機會可以載入接下來高機率會去存取的 key 值。
在這種 layout 下搜尋時，我們總是從第一個 element 開始($k = 1$)，然後每次都根據與 key 值比較的結果，遞進到 $2K$ 或者 $2K + 1$ 的位置。如下圖所示。
### 建立 Eytzinger array
以下的程式碼展示了如何在 O(n) 時間中建立 Eytzinger array。
```cpp
static int max;
int eytzinger(int *src_arr, int *arr, int i, int k, int n) {
if (k <= n) {
i = eytzinger(src_arr, arr, i, 2 * k, n);
arr[k] = src_arr[i++];
i = eytzinger(src_arr, arr, i, 2 * k + 1, n);
}
return i;
}
int *prepare(int *src_arr, int n) {
max = src_arr[n - 1];
int *arr = calloc(n + 1, sizeof(int));
eytzinger(src_arr, arr, 0, 1, n);
arr[0] = -1;
return arr;
}
```
一些需留意的細節:
* 由於 Eytzinger array 會打散原本的排序結構，為了處理 lower bound 不存在的案例，可以在 prepare 時記錄下 array 中的最大值並作為特殊情況處理
* 該函式會按照 index 1 到 n 的順序將對應的 `src_arr` 內容寫到 `arr`
* 當 `arr[k]` 被寫入，屬於其左子樹的對應位置都已經被寫入
* 因此，儘管是 recursive 設計的，但因為記憶體讀取都是順序的，而寫入有 O(logn) 空間會落在相同的 memory block，建立 Eytzinger array 的速度很快
* index 0 被刻意保留下來提供一些額外的優化空間，例如我們可以將 `lower_bound` 找不到答案時的回傳值放在 index 0(雖然為了優化的一致性，我們直接將其作為特例處理)
* 也因此 Eytzinger array 也可以被看作是 1-based 的陣列
### 基於 Eytzinger array 的搜尋法
則搜尋 Eytzinger array 的方式是從 index `k = 1` 開始，將對應 key 值與 `val` 做比較。如果 key 值較大，則往左子樹搜尋 `k := 2*k`;若 key 值較小，則往右子樹搜尋 `k := 2*k+1`。歸納邏輯為以下的程式碼:
```cpp
int k = 1;
while (k <= n)
k = 2 * k + (arr[k] < val);
```
但下面我們會遇到的問題是，迴圈會一直搜尋到 `k` 落在 `n` 之外，但我們想找的 lower bound 應該只是搜尋路徑的其中一個節點。這該怎麼辦呢?
訣竅是: 假設當我們在某一節點找到 lower bound 之後，首先在 lower bound 節點上這數字是 `>= val` 的，因此我們會往左邊(`k := 2*k`)搜尋。之後遇到的數字因為不是 lower bound，數字必然滿足 `< val`，則剩餘的搜尋都會往右(`k := 2*k+1`)，直到搜尋到 leaf 節點為止。那麼我們只要將取消這一系列的右搜尋加上一次的左搜尋，就可以找到正確的 lower bound。
而由於 `k` 計算上的性質，`k := 2*k` 表示將 k 左移 1 bit，相當於是在低位元加上一個 0，而 `k := 2*k + 1` 是將 k 左移 1 bit 並加上 1。因此我們只要通過 k 的二進位表示之低位元有多少個連續的 1，就能夠得知最後有幾次的右搜尋需要取消，然後再加上一次的左搜尋。實作上也就是:
```cpp
int lower_bound(int *arr, int n, int val) {
if (max < val)
return -1;
int k = 1;
while (k <= n)
k = 2 * k + (arr[k] < val);
k >>= __builtin_ffs(~k);
return arr[k];
}
```
得到的結果如下圖的綠線，結果似乎......在 array size 變大時這個作法反而對效能產生負面的影響?
其中的主要原因是 Eytzinger 在大的陣列下並沒有直接獲得 spatial locality 的優勢。觀察搜尋方式我們可以發現，每次的 iteration `k` 會呈 2 倍的成長，因此存取的記憶體會比上一次更遠兩倍，導致比較的最後幾個元素很可能有更高的機會不再位於同一個 cache line 中。而 sorted array 上的查找最後的 iteration 反而會收斂在一個小範圍的空間中，硬體的 prefetch 能預先載入接下來要存取的元素之機會更好。
另一個劣勢是迴圈條件，與之前 branchless 提到的 `while (len > 1)` 不同，這裡的 `while (k <= n)` 編譯器並無法優化為 constant 次數的 iteration 進而排除掉 branch。
```cpp
void clean(int *arr)
{
if (arr) {
free(arr);
}
}
```
最後，由於有額外配置的空間，這裡需要實作 `clean` 來完成回收。
### 透過 Software Prefetch 優化 Eytzinger array 搜尋
既然 hardware 沒辦法自動 prefetch，我們可以效法前面的方式使用 software prefetch。值得留意的是與之前要 prefetch 左右的下個 key 不同，在 Eytzinger array 中下一個可能的 key 會在 `2k` 或 `2k+1`，兩者很可能是落在同一 cacheline 的鄰近關係，因此我們只需要一道 prefetch instruction 即可。甚至觀察再往下一代的關係:
```
4k
/
2k
/ \
4k + 1
k
4k + 2
\ /
2k + 1
\
4k +3
```
可以發現連孫代都可能在一個 prefetch 中包含進來。繼續往下推導的話，根據 64 bytes cache / 4 bytes integer = 16，最佳情況下我們能在一次的 prefetch 中載入到曾曾孫層級 (16k ~ 16K + 15) 所需要所有 key 至 cache。不過以上討論是最佳情形，因為 16k ~ 16K + 15 範圍的元素可能橫跨兩個 cacheline。為了保證可以在單一次 prefetch 就載入 16k ~ 16K + 15 區間中的所有元素，我們可以將陣列對齊 cacheline 的開頭。則因為 cacheline 的 64 bytes 正好可以被 16 整除，上述的目的得以達成。
我們可以修改 prepare 來配置對齊 cache 的空間，並在每個 iteration 都預先 prefetch，參考如下的實作。
```cpp
static int max;
int *prepare(int *src_arr, int n) {
max = src_arr[n - 1];
// assume the multiplication never overflow
int *arr = aligned_alloc(64, (n + 1) * sizeof(int));
eytzinger(src_arr, arr, 0, 1, n);
arr[0] = -1;
return arr;
}
int lower_bound(int *arr, int n, int val) {
if (max < val)
return -1;
int k = 1;
while (k <= n) {
__builtin_prefetch(arr + k * 16);
k = 2 * k + (arr[k] < val);
}
k >>= __builtin_ffs(~k);
return arr[k];
}
```
這讓我們得到如下圖的結果，可以看到 Eytzinger array 所帶來的 locality 優勢終於得以體現。留意到幾個額外的細節:
1. 我們可以 prefetch 更後續 iteration 的 key 值，並且仍然可能只使用 1 道 prefetch instruction 完成。可以直接將 `__builtin_prefetch(arr + k * 16)` 更改成 `__builtin_prefetch(arr + k * 32);`，然後我們依賴硬體的 prefetch 去把鄰近的 cache line 也預先載入。不過這方法也可能無效，取決於硬體實際的行為
2. 結束 while 迴圈前的幾個 prefetch 實際上是無效的，且 prefetch 的範圍會超出 array 配置的合法空間。大多數現代的 CPU 可以將這些非法位置的 prefetch 轉換成 nops，不過也有可能 CPU 無此能力而導致錯誤或是速度變慢。我們也可以選擇微調程式邏輯讓最後的幾個 iteration 不要 prefetch
3. prefetch 並不是毫無代價的，我們實際上是在用多餘的 memory bandwidth 換取延遲的降低。但如果平台上有其他依賴於密集的記憶體存取的任務，它們將顯著影響彼此的性能
### 消除 branch 的成本
我們有提到之前針對 Eytzinger array 實作的另一個問題是迴圈條件無法優化為 constant 次數的 iteration，因此 branch 隱含的預測錯誤將會導致效能的下降。放大之前實驗的效能圖更能看到問題的根源所在:
我們可以發現在 Eytzinger 方法上，latency 的量測結果呈現鋸齒狀，且高處的尖端大約都是在陣列大小 `n` 介於 $2^k$ 和 $2^{k+1}$ 間，且與兩者距離都差不多時。因為該狀況下並不容易預測 loop 會進行 $\lfloor log_2n \rfloor$ 還是 $\lfloor log_2n \rfloor + 1$ 次，因此 branch 預測錯誤的機會提升，進而造成額外的延遲。
一種解決這種問題的方法是刻意將陣列擴展到 2 的次方大小，不過這聽起來有點浪費記憶體空間。比較可接受的方式是，我們可以先執行固定的 $\lfloor log_2n \rfloor$ 次，然後再選擇性地與最後一個元素或某個 dummy 元素(必須保證小於 `val`)進行最後的比較，藉此移除掉最後一個分支。具體的實作方式如下:
```cpp
static int iters;
static int max;
int *prepare(int *src_arr, int n) {
// assume the multiplication never overflow
int *arr = aligned_alloc(64, (n + 1) * sizeof(int));
eytzinger(src_arr, arr, 0, 1, n);
arr[0] = -1;
iters = log2(n + 1);
max = arr[n - 1];
return arr;
}
int lower_bound(int *arr, int n, int val) {
if (max < val)
return -1;
int k = 1;
for (int i = 0; i < iters; i++)
k = 2 * k + (arr[k] < val);
int *loc = (k <= n ? arr + k : arr);
k = 2 * k + (*loc < val);
k >>= __builtin_ffs(~k);
return arr[k];
}
```
可以看到藉此方式得到的實驗結果。
## 基於 B-tree + SIMD 優化 binary search
### Issue
雖然 Eytzinger array + prefetching 的方式看似得到不錯的效能，但此前我們也提到 prefetch 方式需要在 memory bandwidth 上付出的代價。如此看來，減少 prefetch 的次數能夠有效的提昇性能。因為 prefetch cacheline 的次數與樹的高度相關，這代表如果可以降低樹的高度，在使用 cacheline 上的成本也可能相應減少。另一方面，如果可以減少每次 prefetch 的 cacheline 中，實際上不會去比較的值的比例，也許也是助益於效能表現的。
是否有辦法可以降低樹的深度呢? 答案是肯定的。也許你已經想到曾經在資料結構課上學到的 [B-tree](https://en.wikipedia.org/wiki/B-tree)。簡單為不認識 B-tree 的朋友介紹 :wave: : B-tree 可以視為是廣義的 binary tree。對於 order = $k$ 的 B-tree，每個樹的節點中會存在 $k$ 個 key，並存在 $k + 1$ 個 branch。這 k 個 key 在節點上是排序的。對於擁有 k 個 key $(K_0, K_1, ... K_k)$ 的節點，其第 $i$ 個 child 包含的所有 key 值會滿足 $K_{i - 1}$ < key < $K_i$ (如果存在)。 如下圖所示，root 節點的 child 1 所包含的 key 值 9 和 12 皆 < 7 且 > 16。
藉此方式，我們可以將樹的高度從 binary tree 的 $log_2n$ 降低到 $log_kn$，也就是 $log_2k$ 倍。當然這不是毫無代價的，在 B-tree 中搜尋時需要將 key 值與節點上的每個點進行比較，因此實際上該資料結構是將樹的高度和比較次數做 trade-off。那麼用更多的比較來換取樹高度的降低是否是值得的呢? 讓我們慢慢的看下去!
### 建立 B-tree
B-tree 在如題意所述的 binaray search 上可行的事實是因為我們不需新增或者移除新的 key 值，後者操作上(有興趣者可以參閱 [Insertion into a B-tree](https://www.programiz.com/dsa/insertion-into-a-b-tree) 或其他網路資源)有較差的效能以及不友善的記憶體存取方式。反之，我們需要的只是從原始的數據建立出的固定不變(static)之 B-tree，並善用其搜尋上的優勢。可以由如下定義建立符合 Eytzinger array 系統的 B-tree:
* 令 root 節點為 node 0
* 對於 node $k$，其擁有 $B + 1$ 個 child(即樹的 order 為 $B$)，分別的編號是 $k * (B + 1) + i + 1$，$i = 0, 1, ..., B$
簡單來說就是用 level order 來標號，下圖是一個可能建立出的 B-tree 之範例。
從上圖能看到這裡建立的並不是標準的 B-tree(類 B-tree?)! 因為我們並不會去平衡從 root 到 leaf 的深度使其相同。此外，原本我們說 child 包含的 key 值需滿足 $K_{i - 1}$ < key < $K_i$，是因為標準的 B-tree 是沒有重複的 key 出現的。但這裡我們允許重複的情況。對於 binary search 的優化而言，關鍵是最長的深度可以被保證是 $\log_{B+1}n$，這是引入 B-tree 最大原因，而細節上原本就是可以根據使用情境修改的!
為了最佳化對 cache 的利用，這裡定義樹的 order 為 16，因此每個 node 之大小正好符合 `16 * sizeof(int) = 64`，也就是一個 cacheline 的尺寸。受益於此，當我們進行節點上所有 key 值與目標值 `val` 的比較時，對於 key 值的讀取大多將得以在 cache 層級完成。而樹的高度 = prefetch cache 的次數與 binary tree 相比得以降低約 $\log_217 = 4$ 倍。
具體建立 B-tree 的程式碼如下。
```cpp
const int B = 16;
static int nblocks;
static int max;
#define KEY(btree, k, i) btree[k * B + i]
// Get the node index for the child i of node k
static inline int go(int k, int i)
{
return k * (B + 1) + i + 1;
}
static void build(int *src_arr, int *btree, int k, int n)
{
static int t = 0;
if (k < nblocks) {
for (int i = 0; i < B; i++) {
build(src_arr, btree, go(k, i), n);
KEY(btree, k, i) = (t < n ? src_arr[t++] : max);
}
build(src_arr, btree, go(k, B), n);
}
}
int *prepare(int *src_arr, int n)
{
max = src_arr[n - 1];
nblocks = ALIGN_UP(n, B);
/* This is actually a two dimensional array (int btree[nblocks][B]),
* but we allocate as a one dimensional array to make sure that we free
* the space correctly in the bench function. */
int *btree = aligned_alloc(64, nblocks * B * sizeof(int));
build(src_arr, btree, 0, n);
return btree;
}
```
簡單解釋需要留意的細節:
* `go` 用來取得 node k 的第 i 個 child 之 node index
* 我們用 `btree[k][i]` 來表達對應到 node k 的第 i 個 key 值，但實際上 `btree` 是一個一維陣列，這是為了讓本解法與其他方法都具有一致的介面
* 為了實作上的方便，我們會使得 B-tree 中的每個節點都具有 16 個 key，因此當原始 `*src_arr` 中的元素數量不對齊 $B$ 時，會 padding `max` 令其對齊
### 搜尋 B-tree
要在 B-tree 上搜尋 lower bound，我們從 root 開始將目標 `val` 與節點上的所有 key 值 $a_0$, $a1$, ... , $a_{B-1}$ 進行比較，找到最小的 $a_k$ 滿足 $a_k >= val$ (local lower bound)後，前進到其第 i 個 child，並繼續此流程直至搜尋到 leaf 節點。而 global lower bound 可以藉由額外的變數來維護，該值存在於所有 local lower bound 其中深度最深的節點中。實作如下:
```cpp
int lower_bound(int *btree, int n, int val)
{
if (max < val)
return -1;
int k = 0, res = max;
while (k < nblocks) {
int cmp_mask = cmp(btree, k, val);
int i = __builtin_ffs(cmp_mask) - 1;
if (i < B)
res = KEY(btree, k, i);
k = go(k, i);
}
return res;
}
```
* `cmp` 的用途是將節點 `k` 的所有 key 值與 `val` 進行比較，如果 $a_i >= val$，則設置 mask 的 bit i
* 如果 lower bound 不存在的狀況需回傳 `1 << B`
* 通過 `__builtin_ffs` 可以從 mask 推論出最小的 local bound 之 index
而有關 `cmp` 的實作，我們可以選擇用 binary serach 尋找 local lower bound。不過考慮到 B = 16，使用 linear search 事實上足矣。也就是如下的實作:
```cpp
static int cmp(int *btree, int k, int val)
{
int mask = (1 << B);
for (int i = 0; i < B; i++)
mask |= (KEY(btree, k, i) >= val) << i;
return mask;
}
```
選擇這個簡單的實作還有一個重大的好處: 我們可以透過 [SIMD(single instruction, multiple data)](https://en.wikipedia.org/wiki/Single_instruction,_multiple_data) 來加速 `cmp` 的運算。什麼是 SIMD 呢? 顧名思義，該種技術可藉由單一指令(instruction)同時對多個資料平行運算。以此題所設置的實驗場景為例，平台可支援的 SIMD 指令集是 [AVX](https://en.wikipedia.org/wiki/Advanced_Vector_Extensions)，因此允許可以使用 256-bits 的 vector。則我們可以同時平行比較 256 / 32 = 8 個 int。`cmp` 的實作可以調整如下:
```cpp
static int cmp(int *btree, int k, int val)
{
int mask;
__m256i x_vec, y_vec, mask_vec;
x_vec = _mm256_set1_epi32(val);
y_vec = _mm256_load_si256((__m256i *) &KEY(btree, k, 0));
mask_vec = _mm256_cmpgt_epi32(x_vec, y_vec);
int lower = _mm256_movemask_ps((__m256) mask_vec);
y_vec = _mm256_load_si256((__m256i *) &KEY(btree, k, 8));
mask_vec = _mm256_cmpgt_epi32(x_vec, y_vec);
int upper = _mm256_movemask_ps((__m256) mask_vec);
mask = ~(lower | upper << 8);
return mask;
}
```
可以看到效能上又可得到進一步的提昇。
:::info
因為這裡使用的資料結構嚴格來說是特化的 B-tree，原作者以 S-tree 命名之。
:::
### 優化 B-tree 方法
B-tree 方法仍然存在優化的空間。首先，一個可優化的地方是透過 Linux 中的 hugepage。
```cpp
int *prepare(int *src_arr, int n)
{
H = height(n);
max = src_arr[n - 1];
nblocks = ALIGN_UP(n, B) / B;
int sz = ALIGN_UP(nblocks * B * sizeof(int), HUGE_PAGESIZE);
int *btree = aligned_alloc(HUGE_PAGESIZE, sz);
madvise(btree, sz, MADV_HUGEPAGE);
build(src_arr, btree, 0, n);
return btree;
}
```
在現代作業系統中，會使用虛擬記憶體進行記憶體管理。在虛擬記憶體管理的機制下，從 user space 的角度來看每個 process 的記憶體都會是連續的。然而實際上記憶體以 [page](https://en.wikipedia.org/wiki/Page_(computer_memory)) 為單位切分，其通常為 4KB 大小，並由作業系統維護的表將每個 virtual page 各自映射到 physical page，連續的兩個 virtual page 不一定映射到連續的兩段 physical page 空間。且雖然從 user space 看來不同 process 各自擁有很大的一段空間(例如 4GB)，但實體記憶體必然有限，因次存取 virtual address 時其映射的位址並不一定存在正確的資料，可能需要額外從硬碟載入資料的過程。
每次存取 virtual addrss 時，作業系統判斷對應的 physical page 是否已經被載入正確的內容，若否，也就是發生 [page fualt](https://en.wikipedia.org/wiki/Page_fault)，就需要額外從硬碟載入到記憶體上。在此種機制下，可以知道當 page 越小，那麼程式執行過程 page fault 發生的頻率會比較高，需要載入的 page 就會很多，這將對性能產生影響。
對此問題，Hugepage 機制允許將 page 以更大的空間為單位做映射。因此可降低 page fault 發生的頻率，進而提升性能。當然 hugepage 通常只應用在特定的需求上，我們不太會將作業系統上的每個 process 都改成使用 hugepage，因為該機制的 trade-off 是記憶體的[碎片化(segmentation)](https://en.wikipedia.org/wiki/Memory_segmentation)。
在引入 hugepage 後，能看到 B-tree 方法在 array 大小提升的狀況下有更佳的 latency。
:::info
對於前面的方法，其實我們也有機會用 hugepage 去優化。不過在已經涵蓋 prefetch 方式優化的前提下，hugepage 的效益可能不高，因此原文也沒有對此進行額外的實驗。
:::
```cpp
#define HUGE_PAGESIZE (1 << 21)
static int H;
static int height(int n)
{
// grow the tree until its size exceeds n elements
int s = 0, // total size so far
l = B, // size of the next layer
h = 0; // height so far
while (s + l - B < n) {
s += l;
l *= (B + 1);
h++;
}
return h;
}
// Add "H = height(n);" in the prepare function
```
單純引入 hugepage 只能帶來細微的效能提昇，而針對整體的優化，一個可以考量的地方是改寫程式碼以允許編譯器主動的優化。首先，雖然在本文的命題下可能沒有幫助，但如果 binary search 的使用的場景可以假設輸入陣列的大小 `n` 在編譯時期已知，那麼我們可以搭配 C++ 的 `constexpr` 或是 C macro 來實作如上的 `height` 算法，以在編譯時期就獲得 B-tree 的高度。則這將有助於編譯器可以進行 loop unrolling 等等的優化。
```cpp
// __m256i x_vec = _mm256_set1_epi32(val - 1);
static int rank(int *btree, int k, __m256i x_vec)
{
__m256i a = _mm256_load_si256((__m256i *) &KEY(btree, k, 0));
__m256i b = _mm256_load_si256((__m256i *) &KEY(btree, k, 8));
__m256i mask_a = _mm256_cmpgt_epi32(a, x_vec);
__m256i mask_b = _mm256_cmpgt_epi32(b, x_vec);
/* The result of comparison mask_a([0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 ]) and
* mask_b([8 | 9 | 10 | 11 | 12 | 13 | 14 | 15]) will be packed to
* [0 | 1 | 2 | 3 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15] */
__m256i mask_vec = _mm256_packs_epi32(mask_a, mask_b);
int mask = _mm256_movemask_epi8(mask_vec);
// we need to divide the result by two because we call movemask_epi8 on
// 16-bit masks:
return __builtin_ctz(mask) >> 1;
}
```
接著，我們可以調整 SIMD 的使用進一步優化 local lower bound 搜尋的效率。相較於優化以前是將 node 上的 16 個 int 拆半，並各自由 SIMD 計算出 mask 再整合起來，更好的方式是直接利用 SIMD 的 `pack` 在 SIMD register 上組合後再取出 mask，進而減少一次的 `movemask`，後者的運行是緩慢的。實作的細節請看上面展示的 `rank`。
留意到 `rank` 設計與原先 `cmp` 實作上的區別:
* 比較的次序作調換(原本是 `_mm256_cmpgt_epi32(x, y)` 變成 `_mm256_cmpgt_epi32(y, x)`)，因此我們就不需要在最後對 mask 做 bit 反轉，但比較的 `x_vec` 要改成用 `val - 1` 而非 `val` 去設置
* `_mm256_packs_epi32` 會將原本用 32bits 表示的一次比較結果變成 16bits 表示。但 `_mm256_movemask_epi8` 是以 8bit 獲取最後的 mask，因此最後 ctz 的結果需除 2 才正確
要注意的地方是 `_mm256_packs_epi32` 是以交錯的方式組合兩個 mask 的，參見 code 中註解的說明，mask `a1 a2` 與 mask `b1 b2` pack 的結果會是 `a1 b1 a2 b2` 而非 `a1 a2 b1 b2`。因此正確的 mask 按理來說需要透過額外的置換獲得。
但相對於每次搜尋都置換一次 mask，取而代之，可以在 `prepare` 階段直接預先調整好 node 中 key 的順序，則 pack 之後就可以直接得到理想的 mask 了，參見以下的 `permute` 實作，我們在 build 完 B-tree 的每個節點時透過之來改變 key 的順序。
```cpp
static void permute(int *node)
{
const __m256i perm = _mm256_setr_epi32(4, 5, 6, 7, 0, 1, 2, 3);
// Move value of btree[k][4]/btree[k][8 + 4] to a destination vector
__m256i *middle = (__m256i *) (node + 4);
__m256i x = _mm256_loadu_si256(middle);
/* This will let the node
* [0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15]
* being [0 | 1 | 2 | 3 | 8 | 9 | 10 | 11 | 4 | 5 | 6 | 7 | 12 | 13 | 14 |
* 15] */
x = _mm256_permutevar8x32_epi32(x, perm);
_mm256_storeu_si256(middle, x);
}
```
由於節點的 key 被置換過順序，我們需要額外的函式 `update` 來把 mask 對應到正確的位置上。
```cpp
const int translate[16] = {0, 1, 2, 3, 8, 9, 10, 11,
4, 5, 6, 7, 12, 13, 14, 15};
static int update(int res, int *btree, int k, int i)
{
if (i >= B) {
return res;
}
return KEY(btree, k, translate[i]);
}
```
將上述的函式介面總結在一起，則 `prepare` 和 `lower_bound` 的實作如下:
```cpp
int *prepare(int *src_arr, int n)
{
H = height(n);
max = src_arr[n - 1];
nblocks = ALIGN_UP(n, B) / B;
int sz = ALIGN_UP(nblocks * B * sizeof(int), HUGE_PAGESIZE);
int *btree = aligned_alloc(HUGE_PAGESIZE, sz);
madvise(btree, sz, MADV_HUGEPAGE);
build(src_arr, btree, 0, n);
return btree;
}
int lower_bound(int *btree, int n, int val)
{
if (max < val)
return -1;
__m256i x_vec = _mm256_set1_epi32(val - 1);
int k = 0, res = INT_MAX;
for (int h = 0; h < H; h++) {
int i = rank(btree, k, x_vec);
res = update(res, btree, k, i);
k = go(k, i);
}
// the last branch:
if (k < nblocks) {
int i = rank(btree, k, x_vec);
res = update(res, btree, k, i);
}
return res;
}
```
得到的實驗數據則如圖所展示。經過各種技巧所優化的 B-tree 方法最終將比起原先的實作帶來約 15 - 20 % 的效能提昇!
## B+ tree
### Issue
來到這裡，或許你已經眼冒金星了 :dizzy:，不過這可還不是終點。讓我們先喝口水休息一會XD，回顧一下目前在優化的 B-tree 中存在的主要問題:
* 由於正解的 lower bound 可能出現在搜尋路徑上的任一個節點，每個 iteration 都需要進行 `update` 來存取節點上的內容，這造成了一定的成本
* 最後會有額外的 branch，因此存在 branch prediction 的成本
為了解決這些問題，我們將會把儲存的資料結構調整為 [B+ tree](https://en.wikipedia.org/wiki/B%2B_tree)。該種資料結構與 B-tree 最大的區別在於:
* 節點可以區分為 internal node 和 data node(leaf) 兩種，internal node 只有搜尋索引(key)，data node 才存在資料
* 對於 internal node，parent 節點的 key 值會重複出現在 child，節點第 i 個 key 值會與第 i + 1 個 child 的最小 key 值相同
* 對於 data node(leaf)，會儲存 key 對應的資料(key-value map)，並維護額外的指標使得可以找到下一個 leaf node，如下圖的紅點
> 細節暫不贅述，對 B+ tree 有興趣可以參考 [30-11 之資料庫層的核心 - 索引結構演化論 B+樹](https://ithelp.ithome.com.tw/articles/10221111)，寫得蠻清楚的~
此種設計的優勢是搜尋速度上可能更快，因為 internal 只存 key;且因為 leaf 呈現 linked list 而自帶排序，在範圍搜尋上效率更好。當然天下沒有白吃的午餐，在 B+ tree 中我們付出的額外的成本是 internal nodes 中重複的 key 所需的空間。
則在本題的狀況下，B+ tree 對於前面的問題可以良好的解決:
* lower bound 將會存在搜尋路徑上的最後一個節點中，或是下一個 leaf 節點的第一個 key，因此不再需要頻繁的 `update`
* 因為 B+tree 成長的方向是由下而上，其從 root 到任一個 leaf 的深度將會是固定的，因此沒有額外的 branch
### 建立 B+tree
與之前的做法相同，我們透過陣列來表示被建立的 B+tree。樹上的 key 從最深的 layer 0 (也就是 leaf 所在的 layer)開始至 root 依序儲存。
```cpp
static int max;
static int H;
static int S;
static int *offset;
#define DIV_UP(n, d) (((n) + (d) -1) / (d))
#define BLOCKS(n) DIV_UP(n, B)
/* For n key values, count the number of keys on the previous
* layer. The calculation comes from the fact that for a group
* of B + 1 child node, one node with B keys will be needed. */
static int prev_keys(int n)
{
return DIV_UP(BLOCKS(n), B + 1) * B;
}
static int height(int n)
{
return (n <= B ? 1 : height(prev_keys(n)) + 1);
}
/* This must be initialized first because we'll use
* it for prepare later. */
// H = height(n);
```
`prev_keys` 函式可以計算在有 `n` 個 key 的情形下，它的前一層需要多包含其中的多少個 key。這計算的根據是我們每有 B + 1 個 child node，其就會有共同一個 parent node，而一個 parent node 將包含 B 個 key。因此 $\text{prev_key} = \lceil \frac{\lceil\text{key / B}\rceil}{B + 1} \rceil * B$。
則一直遞迴到 `key <= B` 的時候，表示該層只存在 1 個 node，也就是 root 了，此時計算已經遞迴了多少次則可以推算出 B+tree 的高度為何。
```cpp
/* For the tree with n keys, get the offset set of the starts
* oflayer h starts (layer 0 is the largest and the leaf node) */
static int get_offset(int n, int h)
{
int k = 0;
while (h--) {
k += BLOCKS(n) * B;
n = prev_keys(n);
}
return k;
}
// The size of tree will equal to the start offset of (non-existent) layer H
// S = get_offset(n, H);
```
透過前述的特性則可以實作 `get_offset` 函式得到指定高度的起始 offset。如果想知道整個樹共有多少個 key 以得知需要動態配置的空間大小，可以用 `get_offset(n, H)` 來獲取。
注意到 `get_offset` 的計算根據 n 和 h 會是固定的輸出結果，因此我們之後能夠計算一次後直接建立表格，再透過查表來避免重新計算 offset 的成本。
```cpp
int *prepare(int *src_arr, int n)
{
max = src_arr[n - 1];
/* This must be initialized first because we'll use
* it for prepare later. */
H = height(n);
/* Construct the table to get the offset of each height faster */
offset = malloc((H + 1) * sizeof(int));
for (int h = 0; h <= H; h++) {
offset[h] = get_offset(n, h);
}
// The size of tree will equal to the start offset of (non-existent) layer H
S = offset[H];
int sz = ALIGN_UP(sizeof(int) * S, HUGE_PAGESIZE);
int *btree = aligned_alloc(HUGE_PAGESIZE, sz);
madvise(btree, sz, MADV_HUGEPAGE);
memcpy(btree, src_arr, sizeof(int) * n);
for (int i = n; i < S; i++)
btree[i] = max;
```
因此 prepare 的實作上如前述，我們可以預先建立 table 來存放每個高度對應到陣列中的 offset。之後配置出 B+tree 所需要的陣列空間後，接著就可以填上每個 element。layer 0 的部分顯而易見，就等同於 src_arr 本身。
```cpp
/* Construct the tree layer by layer */
for (int h = 1; h < H; h++) {
int offset_cur = offset[h];
int offset_range = offset[h + 1] - offset_cur;
for (int i = 0; i < offset_range; i++) {
/* k is the index of node */
int k = i / B;
/* j will range from 0 - 15, which is the sub-index of
* the corresponding node */
int j = i - k * B;
/* Descend to the offset of the right node */
k = k * (B + 1) + j + 1;
/* and then go left until reaching the leaf */
for (int l = 0; l < h - 1; l++)
k *= (B + 1);
// The smallest key on the leaf will the key
btree[offset_cur + i] = (k * B < n ? btree[k * B] : max);
}
}
```
剩餘的 internal node 部分可以一層一層進行。對於每個位置的值，找法是先找到對應 key 存在的 leaf node 之 node index，然後再取該 node 中最小的 key。
```cpp
// Don't permute layer 0 to avoid index translation on it
for (int i = offset[1]; i < S; i += B)
permute(btree + i);
return btree;
}
```
和之前的 B+tree 類似，為了待會能使用 `pack` 來加速，我們先透過 `permute` 調整 key 的 index 順序。注意到 `offset(0)` 也就是 leaf 層的 key 並未被 permute，原因是為了避免最後需要之前在 `update` 的實作用到的 index 的轉換，後者造成額外的時間成本。因此當搜尋至 leaf 層時，我們會回到最早先將 node 分成前後兩半比較，再將 mask 組合起來的方式。
### 搜尋 B+tree
```cpp
int lower_bound(int *btree, int n, int val)
{
if (max < val)
return -1;
int k = 0;
__m256i x_vec = _mm256_set1_epi32(val - 1);
for (int h = H - 1; h > 0; h--) {
int i = permuted_rank(btree, offset[h] + k, x_vec);
k = k * (B + 1) + i * B;
}
int i = direct_rank(btree, k, x_vec);
return btree[k + i];
}
```
則搜尋上很簡單，從 root 節點開始，我們在節點上找到 local lower bound 後，前進到下個對應分支的節點，直到 leaf node，最後再從 leaf node 上獲取正確的 global lower bound 即可。這讓我們避免了之前在 B-tree 實作下，每一次找到新的 local lower bound 都要 `update` 的成本，也消除了最後的 branch prediction。
`permuted_rank` 就類似於之前的 `rank`，而 `direct_rank` 則類似於 `cmp`，細節不再贅述，可以參考 [binary-search-playground](https://github.com/RinHizakura/binary-search-playground)。從實驗數據上，我們能看到 B+tree 在效能上相較於 B-tree 具有更理想的效率!
```cpp
void clean(int *btree)
{
if (btree) {
free(btree);
}
if (offset) {
free(offset);
offset = NULL;
}
}
```
最後 `clean` 的部分，除了 `btree` 之外也需要釋放用來存 `offset` 的空間。
## Bonus: Shar's algorithm
TODO
https://probablydance.com/2023/04/27/beautiful-branchless-binary-search/
https://orlp.net/blog/bitwise-binary-search/
## 小結
以上就是對 [Optimizing Binary Search](https://youtu.be/1RIPMQQRBWk) 這部影片內容的所有整理和討論。但是即便我們在過程中認識許多特殊的資料結構，也學習到不少妥善運用現代硬體架構特性的技巧，進步的可能性卻是沒有極限的 :smile: ! 有興趣的讀者可以在原作者書中的 [Data Structures Case Studies](https://en.algorithmica.org/hpc/data-structures/) 一章深入了解更多的優化 binary search 的方式。
本文中的所有實作都可以在 [binary-search-playground](https://github.com/RinHizakura/binary-search-playground) 找到。值得注意的是，這裡是將原作者的實作用 C 重寫，並修正其中的錯誤(原文在 B-tree 和 B+tree 的實作上存在正確性的問題，大家可以仔細找找發生在哪XD)，但並未真正測量運行時間。考慮到編譯器優化的差異，以及為了解釋方便演算法而移除掉的 `const expr` 等等因素，binary-search-playground 中的範例實際執行效率可能和原作者的實驗數據有出入。這部份目前就先留給感興趣的朋友們親自去研究了~
:::danger
註: [binary-search-playground](https://github.com/RinHizakura/binary-search-playground) 目前已經嘗試整合測試用的 [python script](https://github.com/RinHizakura/binary-search-playground/blob/main/scripts/perf_plot.py)，可以用來得到類似原文效果的圖表，然而執行後就可以發現本文中提供的 code 所得的實驗結果與原文存在差異。可想而知無論是實驗環境(需小心對待 CPU isolation, cache 的影響)和移植到 C 的實作手法都需要再調整，請讀者務必留意。
:::
## Reference
* [Fast(er) binary search in Rust](https://www.bazhenov.me/posts/faster-binary-search-in-rust/)
* [binary-search-cppcon](https://github.com/CppCon/CppCon2022/blob/main/Presentations/binary-search-cppcon.pdf)
* [Eytzinger Binary Search](https://algorithmica.org/en/eytzinger)
* [Implicit Static B-trees](https://algorithmica.org/en/b-tree)
* [ARRAY LAYOUTS FOR COMPARISON-BASED SEARCHING](https://arxiv.org/ftp/arxiv/papers/1509/1509.05053.pdf)
* [Algorithms for Modern Hardware - Binary Search](https://en.algorithmica.org/hpc/data-structures/binary-search/)
* [sslotin/amh-code - binsearch](https://github.com/sslotin/amh-code/tree/main/binsearch)
Sign in with Wallet
Connect another wallet