# C03: clz
###### tags: `sysprog2017`
:::info
主講人: [jserv](http://wiki.csie.ncku.edu.tw/User/jserv) / 課程討論區: [2017 年系統軟體課程](https://www.facebook.com/groups/system.software2017/)
:mega: 返回「[進階電腦系統理論與實作](http://wiki.csie.ncku.edu.tw/sysprog/schedule)」課程進度表
:::
## [De Bruijn 數列](https://en.wikipedia.org/wiki/De_Bruijn_sequence)
* Stanford 大學博士候選人 Sean Eron Anderson 維護了 [ Bit Twiddling Hacks](http://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2) 頁面，裡頭提及許多 bit-wise 的密技。在同一個網頁的[不同部份](http://graphics.stanford.edu/~seander/bithacks.html#ZerosOnRightMultLookup) 有 De Bruijin 數列的用途:
> The constant 0x077CB531UL is a de Bruijn sequence, which produces a unique pattern of bits into the high 5 bits for each possible bit position that it is multiplied against. When there are no bits set, it returns 0. More informaton can be found by reading the paper Using de Bruijn Sequences to Index 1 in a Computer Word by Charles E. Leiserson, Harald Prokof, and Keith H. Randall.
* 課程影片：
{%youtube 6hR1X5d_1VU%}
這個數列的精神如下：
「給定 m 種字母，以及字串長度 k, 這個數列 B(m, k) 向左做旋轉後, 取出前 k 個字母。 恰好可以湊出這 m 種字母在字串長度 k 下的所有可能組合。」
舉例來說，引用[這篇](http://www.stmintz.com/ccc/index.php?id=306404)如下 ：
```
My small De Bruijn Test of 130329821:
0x07C4ACDD
00000111110001001010110011011101B
bit 31 - bit 27 00000 0
bit 30 - bit 26 00001 1
bit 29 - bit 25 00011 3
bit 28 - bit 24 00111 7
bit 27 - bit 23 01111 15
bit 26 - bit 22 11111 31
bit 25 - bit 21 11110 30
bit 24 - bit 20 11100 28
bit 23 - bit 19 11000 24
bit 22 - bit 18 10001 17
bit 21 - bit 17 00010 2
bit 20 - bit 16 00100 4
bit 19 - bit 15 01001 9
bit 18 - bit 14 10010 18
bit 17 - bit 13 00101 5
bit 16 - bit 12 01010 10
bit 15 - bit 11 10101 21
bit 14 - bit 10 01011 11
bit 13 - bit 9 10110 22
bit 12 - bit 8 01100 12
bit 11 - bit 7 11001 25
bit 10 - bit 6 10011 19
bit 9 - bit 5 00110 6
bit 8 - bit 4 01101 13
bit 7 - bit 3 11011 27
bit 6 - bit 2 10111 23
bit 5 - bit 1 01110 14
bit 4 - bit 0 11101 29
bit 3 - bit 31 11010 26
bit 2 - bit 30 10100 20
bit 1 - bit 29 01000 8
bit 0 - bit 28 10000 16
0x07C4ACDD is 5bit de Bruin Sequence
```
可以發現 `0x07C4ACDD` 就是個長度 5, 單字表是 {0, 1} 的一個 De Bruijn 數列。 因為他在向左旋轉 32 次之後, 前 5 個 bit 恰好湊齊了 0 ~ 31, 而 0~31 也恰代表了枚舉了所有 5 個 {0, 1} 可以組成的字串。
## De Bruijn 數列: bit shift 版
De Bruijn 序列一開始的精神是 「rotation 1, 2 ... 31 次之後, 就湊滿 0~31」。那可否省略 rotation 31 次的操作, 而是找「 left bit shift 0, 1, ...31 次之後， 剛好湊滿 0~31 ?」 答案是可以的。 實際上這跟上面的 De Bruijn 是同一個問題 ： 只要找到的 De Bruijn 序列「前 5 個 bit 是 0」， 那做 rotation 取前 5 個bit, 跟做 left bit shift 取前 5 個 bit 是一樣的。 0 都會從右方補進來。
比如 `0x077CB531U`
```
00000111011111001011010100110001
l-shamnt binary decimal
0 00000 0
1 00001 1
2 00011 3
3 00111 7
4 01110 14
5 11101 29
6 11011 27
7 10111 23
8 01111 15
9 11111 31
10 11110 30
11 11100 28
12 11001 25
13 10010 18
14 00101 5
15 01011 11
16 10110 22
17 01101 13
18 11010 26
19 10101 21
20 01010 10
21 10100 20
22 01001 9
23 10011 19
24 00110 6
25 01100 12
26 11000 24
27 10001 17
28 00010 2
29 00100 4
30 01000 8
31 10000 16
```
1. 左移 0 次, 取前5個 bit, 得到 0
2. 左移 1 次, 取前5個 bit, 得到 1
3. 左移 2 次, 取前5個 bit, 得到 3
4. ...
這種序列有什麼好處呢 ？ 對於0 ~ 31 的任一個平移量, 都可以 map 到一個 0~ 31 的數。 但是平移又跟 2 的次方有關。 換個說法就是:
「 透過這個數列, 可以把像這種資料 :
```
00000000000000000000000000000001
00000000000000000000000000000010
00000000000000000000000000000100
00000000000000000000000000001000
00000000000000000000000000010000
00000000000000000000000000100000
00000000000000000000000001000000
00000000000000000000000010000000
00000000000000000000000100000000
00000000000000000000001000000000
00000000000000000000010000000000
00000000000000000000100000000000
...
10000000000000000000000000000000
```
++hash++ 到[0, 31] 當中， 而方法就是：
```C
h(value) = (value * 0x077CB531U) >> 27
```
不過 hash 順序不一定是照位元數目。 這樣其實也可以說是「把 bit shift hash 進 [0, 31] 當中」。
既然順利 hash 進去, 所以就可以動點手腳， 用查表求 2 的次方。
考慮下面這個程式 ：
```C
unsigned int v; // find the number of trailing zeros in 32-bit v
int r; // result goes here
static const int MultiplyDeBruijnBitPosition[32] =
{
0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8,
31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
};
r = MultiplyDeBruijnBitPosition[((uint32_t)((v & -v) * 0x077CB531U)) >> 27];
```
1. `(v & -v)` ： 會把「最低位的1」 過濾出來, 除了他以外的 bit 通通都會變成 0。 原理是計概學過的「二的補數」的換算方法, 右邊開始往左碰到第一個不是0的數之後, 後面全部變號。 所以最低位的1往左的位元全部相反, 右邊又都是0, AND 起來通通消光光
2. 將`(v & -v)` 乘上前述的數值 (0x07C4ACDD) : 因為這個東西恰好是 $2^N$ 所以恰好在做左移。 剛剛知道這個數列的特性就是「對給定的平移, 可以在 0~31 中找到唯一的結果」。 所以只要建個大小 32 的表， 把平移量 map 到的欄位， 擺上對應的位數，即可實作一個查表的count trailing zero (ctz) 了。
3. 這支程式就在做這件事。 來源是[這裡](http://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2)
## De Bruijin 數列 - clz 版
把這樣的啟發應用到 clz. 這裡的問題變成這樣：
「能不能找到一個數 M , 使得做了以下的動作 ：
```C
h(value) = (value * M) >> 27
```
之後， 把所有形如下面的 a :
```
00000000000000000000000000000001
00000000000000000000000000000011
00000000000000000000000000000111
00000000000000000000000000001111
00000000000000000000000000011111
00000000000000000000000000111111
00000000000000000000000001111111
00000000000000000000000011111111
00000000000000000000000111111111
00000000000000000000001111111111
...
01111111111111111111111111111111
11111111111111111111111111111111
```
hash 進 [0, 31] 當中 ？」
先解釋一下為什麼要是左邊一排 0, 右邊一排 1 的形式 , 因為對於任意 32 bit 整數由下面這個步驟
```C
value |= value >> 1;
value |= value >> 2;
value |= value >> 4;
value |= value >> 8;
value |= value >> 16;
```
把 Most Significant Bit 以左的位元全部設成 0, 以右全部設為 1。
這個正名很容易， 因為在這個過程當中, 最高位的1 被向左複製了 1 + 2 + 4 + 8 + 16 = 31 次, 大概像這樣 ：
```
...001001101.... // 某個 32 bit 整數中, MSB 附近的數。 最左邊的 1 是 MSB
11111111111111111111111111111111 // 剛剛複製出來的 1
// 兩者做位元或
```
但是 MSB 右邊的 bit 數目一定小於 31 （因為他是 32 bit 整數嘛~）
對於 64 bit 也可以做類似的事，只是需要比較多步：
```C
value |= value >> 1;
value |= value >> 2;
value |= value >> 4;
value |= value >> 8;
value |= value >> 16;
value |= value >> 32;
```
接下來有兩種選擇：
1. 因為 value 的 MSB 以右都是 1，所以做
```C
value -= value >> 1;
```
就可以讓數值從本來的這樣：
```
00000000000000000000000000000001
00000000000000000000000000000011
...
01111111111111111111111111111111
11111111111111111111111111111111
```
變成這樣：
```
00000000000000000000000000000001
00000000000000000000000000000010
00000000000000000000000000000100
00000000000000000000000000001000
00000000000000000000000000010000
00000000000000000000000000100000
00000000000000000000000001000000
...
10000000000000000000000000000000
```
這樣問題就回到上面最原始的 De Bruijn 問題（變成 byte-shift 了嘛）。所以重新理解數值的 `log64` 就是用這個概念實做的：
```C
const int tab64[64] = {
63, 0, 58, 1, 59, 47, 53, 2,
60, 39, 48, 27, 54, 33, 42, 3,
61, 51, 37, 40, 49, 18, 28, 20,
55, 30, 34, 11, 43, 14, 22, 4,
62, 57, 46, 52, 38, 26, 32, 41,
50, 36, 17, 19, 29, 10, 13, 21,
56, 45, 25, 31, 35, 16, 9, 12,
44, 24, 15, 8, 23, 7, 6, 5
};
int log2_64 (uint64_t value) {
value |= value >> 1;
value |= value >> 2;
value |= value >> 4;
value |= value >> 8;
value |= value >> 16;
value |= value >> 32;
return tab64[((uint64_t)((value - (value >> 1 ))*0x07EDD5E59A4E28C2)) >> 58];
}
```
注意最下面的 `value - (value >> 1)` ，就把它變回 對 De Bruijn 序列做 byte-shift。
然後，因為：
(value - value >> 1) * 0x100000
展開來：
(value) * 0x100000 - value * 0x10000
所以這也等於：
(value) * 0x011111
所以實際上，變成 0x01111 這種形式後，直接乘 De Bruijn 序列也是可以維持性質的。這就是下面的做法：
2. 直接把 32 個形如 0000...0011111...1 的數字(我就叫他 a 好了) hash 進 [0:31]。假定我們用之前「一個乘法 + 取高 5 位」的作法，也就是我們要找到一個 M ，並且 ：
1. 用上面的作法, 把一個 32 bit 整數弄成如上的形式
2. 把處理過的整數, 帶進 `h(a) = a*M >> 5` 找查對應的值
3. 因為一個 a 只會對應 [0, 31] 唯一一個數, 而一個 a 也恰好只有唯一的一個 leading zero 數目。 所以只要建個 table[32] 陣列, 把 hash 算出來數值的位置, 擺上對應的 leading zero 數目, 就可以完成一個查表的快速 clz 了。
但是還有一個問題 ： 這個 M 存不存在呢 ? 答案是存在的。 其中一個 M 就是
0x07C4ACDD
```
num of 1 byte after mult decimal
1 00000 0
2 00010 2
3 00110 6
4 01110 14
5 11110 30
6 11101 29
7 11011 27
8 10111 23
9 10000 16
10 00001 1
11 00011 3
12 01000 8
13 10001 17
14 00100 4
15 01001 9
16 10100 20
17 01010 10
18 10101 21
19 01011 11
20 11000 24
21 10010 18
22 00101 5
23 01100 12
24 11010 26
25 10110 22
26 01101 13
27 11100 28
28 11001 25
29 10011 19
30 00111 7
31 01111 15
32 11111 31
```
可以發現再做這個 hash 方法恰好把上述所有的 a 成功對應到 0~31 當中了！ 所以就可以建一個由右邊查到左邊的表， 然後很快的找出 1 的數目。
比如說以下 log2 的實作:
```C
const int tab32[32] = {
0, 9, 1, 10, 13, 21, 2, 29,
11, 14, 16, 18, 22, 25, 3, 30,
8, 12, 20, 28, 15, 17, 24, 7,
19, 27, 23, 6, 26, 5, 4, 31
};
int log2_32 (uint32_t value)
{
value |= value >> 1;
value |= value >> 2;
value |= value >> 4;
value |= value >> 8;
value |= value >> 16;
return tab32[(uint32_t)(value * 0x07C4ACDD) >> 27];
}
```
可以注意到 `table32[32]` 跟上面的表根本一樣， 只是把表左邊的數字通通減 1 , 因為 「右邊 1 的數目 - 1」就是取 log2 的結果。
但是這引起了一個問題：當初 De Bruijin 數列定義的是「旋轉」，byte shift 之所以成立，是因為剛好有前幾個 byte-shift 跟「旋轉」 等價，所以仍然成立。但這裡似乎完全不一樣 -- 乘上形如 0x00FF 這種數看起來一點也不像是旋轉的特例。
這個原理其實就跟上面提到的一樣
因為：
(value - value >> 1) * 0x100000
展開來：
(value) * 0x100000 - value * 0x10000
所以這也等於：
(value) * 0x011111
所以仍然如果乘上形如 00..01111...1 這種形式的數，也是可以讓性質成立的。上述表格能變形，把計算 log 轉變成設計 minimal hash 的問題了。
## iterative
```C
int clz(uint32_t x) {
int n = 32, c = 16;
do {
uint32_t y = x >> c;
if (y) { n -= c; x = y; }
c >>= 1;
} while (c);
return (n - x);
}
```
計算要左移幾次才會變 0。算完之後再用這個數有幾位去減。
每次都把數字折半，對位數做 binary search。
## binary search
```C
int clz(uint32_t x) {
if (x == 0) return 32;
int n = 0;
if (x <= 0x0000FFFF) { n += 16; x <<= 16; }
if (x <= 0x00FFFFFF) { n += 8; x <<= 8; }
if (x <= 0x0FFFFFFF) { n += 4; x <<= 4; }
if (x <= 0x3FFFFFFF) { n += 2; x <<= 2; }
if (x <= 0x7FFFFFFF) { n += 1; x <<= 1; }
return n;
}
```
一樣是 binary search，只是用 bit mask 做。
更進一步的解釋是這樣：一個二進位的數可以表示成：
$$d_{1} \cdot 16 + d_{2} \cdot 8 + d_{3} \cdot 4 + d_{4} \cdot 2 + d_{5} \cdot 1$$
所以一堆 `if` 只是在判斷 $d_{1}$ ~ $d_{5}$ 哪些為 0, 哪些不為 0
## byte shift
```C
int clz(uint32_t x) {
if (x == 0) return 32;
int n = 1;
if ((x >> 16) == 0) { n += 16; x <<= 16; }
if ((x >> 24) == 0) { n += 8; x <<= 8; }
if ((x >> 28) == 0) { n += 4; x <<= 4; }
if ((x >> 30) == 0) { n += 2; x <<= 2; }
n = n - (x >> 31);
return n;
}
```
精神仍然跟 binary search 一樣。
## 作業要求
* 閱讀 [重新理解數值](https://hackmd.io/s/BkRKhQGae) 裡頭關於 count leading zero ([clz](https://en.wikipedia.org/wiki/Find_first_set#CLZ)) 的描述，設計實驗，比較 32-bit 數值對以下實做的效能差異：
* recursive version
* iteration version
* binary search technique
* byte-shift version
* Harley's algorithm
* 除了在 [重新理解數值](https://hackmd.io/s/BkRKhQGae) 列出的程式，詳閱[劉亮谷的實驗紀錄](/s/BJBZn6Q6)，予以重現並解釋個別實作的運作原理
* 解說影片: [Count Leading Zero](https://www.youtube.com/watch?v=DRkXFjLfaVg) [必看!]
* 走訪全部的 32-bit 數值，用上述演算法帶入計算 clz 值，先驗證正確性，如果演算法不正確，試圖改正
* 在 GitHub 上 fork [clz-tests](https://github.com/sysprog21/clz-tests)，以此為基礎進行以下調整 (如在 9 月 23 日就 fork 過，那請重新 fork)
* 用 C99 或 C11 改寫程式碼，其中 C11 的 [_Generic](http://www.robertgamble.net/2012/01/c11-generic-selections.html) 非常好用，詳情可見 [你所不知道的 C 語言：前置處理器應用篇](https://hackmd.io/s/S1maxCXMl)
* 比照 [phonebook](/s/rJYD4UPKe) 和 [compute-pi](/s/HkiJpDPtx)，設計一套 benchmark suite，得以針對所有的 32-bit 數值進行個別 clz 實做效能的分析，並透過 gnuplot 製圖
* 要附上個別數值實際耗費時間，不能只列出平均值
* 提供落點分析圖，類似 [tcp-anaysis](https://upload.wikimedia.org/wikipedia/commons/6/65/Gnuplot_tcp_analysis.png) ([with-code](https://commons.wikimedia.org/wiki/File:Gnuplot_tcp_analysis.png))
* 為了避免編譯器最佳化的影響，務必指定編譯參數 `-O0` 來抑制最佳化
* 找至少 3 個案例，說明 clz 的應用場合
* 示範: [A Fast Hi Precision Fixed Point Divide](http://me.henri.net/fp-div.html)
* 提示：透過 [Google Books](https://books.google.com/) 可以找到一些專門探討 optimization 的書籍，裡頭會有 clz 的應用
* 將你的觀察、上述要求的解說、應用場合探討，以及各式效能改善過程，善用 gnuplot 製圖，紀錄於「[作業區](https://hackmd.io/s/HyxQTaZj-)」
## 截止日期
* Oct 8, 2017 (含) 之前
* 越早在 GitHub 上有動態、越早接受 code review，評分越高
## 共筆選讀
- [ ] [baimao8437](https://hackmd.io/s/BkZOVSIcx)
* 可看到 byte 和 binary 所用 cycles 最少，而 harley 相對平穩
- [ ] [0xff07](https://hackmd.io/s/rkTYU0vKx#): 數學之美
## 參考資料
:::warning
引用過往作業成果時，務必標註出處 (包含 GitHub 帳號或人名)，連帶相關的超連結
:::
* [2016 年秋季班課程作業區](https://hackmd.io/s/H1B7-hGp): 提供錄影解說
Sign in with Wallet
Connect another wallet