C03: clz

tags: `sysprog2017`

主講人: jserv / 課程討論區: 2017 年系統軟體課程

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

De Bruijn 數列

Stanford 大學博士候選人 Sean Eron Anderson 維護了 Bit Twiddling Hacks 頁面，裡頭提及許多 bit-wise 的密技。在同一個網頁的不同部份有 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.

課程影片：
Image Not Showing Possible Reasons
- The image file may be corrupted
- The server hosting the image is unavailable
- The image path is incorrect
- The image format is not supported
Learn More →

這個數列的精神如下：
「給定 m 種字母，以及字串長度 k, 這個數列 B(m, k) 向左做旋轉後, 取出前 k 個字母。恰好可以湊出這 m 種字母在字串長度 k 下的所有可能組合。」

舉例來說，引用這篇如下：

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

左移 0 次, 取前5個 bit, 得到 0
左移 1 次, 取前5個 bit, 得到 1
左移 2 次, 取前5個 bit, 得到 3
…

這種序列有什麼好處呢？對於0 ~ 31 的任一個平移量, 都可以 map 到一個 0~ 31 的數。但是平移又跟 2 的次方有關。換個說法就是:

「透過這個數列, 可以把像這種資料 :

00000000000000000000000000000001
00000000000000000000000000000010
00000000000000000000000000000100
00000000000000000000000000001000
00000000000000000000000000010000
00000000000000000000000000100000
00000000000000000000000001000000
00000000000000000000000010000000
00000000000000000000000100000000
00000000000000000000001000000000
00000000000000000000010000000000
00000000000000000000100000000000

...

10000000000000000000000000000000

hash 到[0, 31] 當中，而方法就是：

   h(value) =  (value * 0x077CB531U) >> 27

不過 hash 順序不一定是照位元數目。這樣其實也可以說是「把 bit shift hash 進 [0, 31] 當中」。

既然順利 hash 進去, 所以就可以動點手腳，用查表求 2 的次方。

考慮下面這個程式：

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];

(v & -v) ：會把「最低位的1」過濾出來, 除了他以外的 bit 通通都會變成 0。原理是計概學過的「二的補數」的換算方法, 右邊開始往左碰到第一個不是0的數之後, 後面全部變號。所以最低位的1往左的位元全部相反, 右邊又都是0, AND 起來通通消光光
將(v & -v) 乘上前述的數值 (0x07C4ACDD) : 因為這個東西恰好是
$2^{N}$ 所以恰好在做左移。剛剛知道這個數列的特性就是「對給定的平移, 可以在 0~31 中找到唯一的結果」。所以只要建個大小 32 的表，把平移量 map 到的欄位，擺上對應的位數，即可實作一個查表的count trailing zero (ctz) 了。
這支程式就在做這件事。來源是這裡

De Bruijin 數列 - clz 版

把這樣的啟發應用到 clz. 這裡的問題變成這樣：

「能不能找到一個數 M , 使得做了以下的動作：

    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 整數由下面這個步驟

    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 也可以做類似的事，只是需要比較多步：

    value |= value >> 1;
    value |= value >> 2;
    value |= value >> 4;
    value |= value >> 8;
    value |= value >> 16;
    value |= value >> 32;

接下來有兩種選擇：

因為 value 的 MSB 以右都是 1，所以做

value -= value >> 1;

就可以讓數值從本來的這樣：

00000000000000000000000000000001
00000000000000000000000000000011

...

01111111111111111111111111111111
11111111111111111111111111111111

變成這樣：

00000000000000000000000000000001
00000000000000000000000000000010
00000000000000000000000000000100
00000000000000000000000000001000
00000000000000000000000000010000
00000000000000000000000000100000
00000000000000000000000001000000
...

10000000000000000000000000000000

這樣問題就回到上面最原始的 De Bruijn 問題（變成 byte-shift 了嘛）。所以重新理解數值的 log64 就是用這個概念實做的：

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 序列也是可以維持性質的。這就是下面的做法：

直接把 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 的實作:

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

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

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

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 一樣。

作業要求

閱讀重新理解數值裡頭關於 count leading zero (clz) 的描述，設計實驗，比較 32-bit 數值對以下實做的效能差異：
- recursive version
- iteration version
- binary search technique
- byte-shift version
- Harley's algorithm
除了在重新理解數值列出的程式，詳閱劉亮谷的實驗紀錄，予以重現並解釋個別實作的運作原理
- 解說影片: Count Leading Zero [必看!]
走訪全部的 32-bit 數值，用上述演算法帶入計算 clz 值，先驗證正確性，如果演算法不正確，試圖改正
在 GitHub 上 fork clz-tests，以此為基礎進行以下調整 (如在 9 月 23 日就 fork 過，那請重新 fork)
- 用 C99 或 C11 改寫程式碼，其中 C11 的 _Generic 非常好用，詳情可見你所不知道的 C 語言：前置處理器應用篇
- 比照 phonebook 和 compute-pi，設計一套 benchmark suite，得以針對所有的 32-bit 數值進行個別 clz 實做效能的分析，並透過 gnuplot 製圖
- 要附上個別數值實際耗費時間，不能只列出平均值
- 提供落點分析圖，類似 tcp-anaysis (with-code)
- 為了避免編譯器最佳化的影響，務必指定編譯參數 -O0 來抑制最佳化
找至少 3 個案例，說明 clz 的應用場合
- 示範: A Fast Hi Precision Fixed Point Divide
- 提示：透過 Google Books 可以找到一些專門探討 optimization 的書籍，裡頭會有 clz 的應用
將你的觀察、上述要求的解說、應用場合探討，以及各式效能改善過程，善用 gnuplot 製圖，紀錄於「作業區」

截止日期

Oct 8, 2017 (含) 之前
越早在 GitHub 上有動態、越早接受 code review，評分越高

共筆選讀

baimao8437
- 可看到 byte 和 binary 所用 cycles 最少，而 harley 相對平穩
0xff07: 數學之美

參考資料

引用過往作業成果時，務必標註出處 (包含 GitHub 帳號或人名)，連帶相關的超連結

2016 年秋季班課程作業區: 提供錄影解說