2017q1 Homework1 (clz)

contributed by < 0xff07 >

原始碼在這裡

作業題目再這裡

Reviewed by <andycom12000>

• 原始實驗結果中的圖二沒有圖例，看不出來哪條線是何者
• Commit cbeb2e806366ee294e95fd8edb0a742cef2752f3中一次包含太多東西，應可分成2~3次的commit區分不同功能
• 可以比較Harley's Algorithm還有Recursive的方式

開發環境

Ubuntu 16.04.2
Linux version 4.8.0-36-generic

gcc version 5.4.0 20160609 (Ubuntu 5.4.0-6ubuntu1~16.04.4)

L1d cache:             32K
L1i cache:             32K
L2 cache:              256K
L3 cache:              3072K

看懂演算法

其實只有那個用到 De Bruijn 數列的演算法看不太懂。所以這裡想辦法弄懂它

首先查資料時看到史丹佛網站有一些有趣的東西。

De Bruijn 數列 - 原始定義

然後再同一個網頁的不同部份找到 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.

• 課程影片：
  

其實這個數列的精神是這樣 ：

「給定 m 個字母， 以及字串長度 k, 這個數列向左做 rotation 後, 取出前 k 個字母。 恰好可以湊出所有這 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
0x07C4ACDD is 5bit de Bruin Sequence

可以發現 0x07C4ACDD 就是個長度 5, 單字表是 {0, 1} 的一個 De Bruijn 數列。 因為他在向左旋轉 32 次之後, 前 5 個 bit 恰好湊齊了 0 ~ 31, 而 0~31 也恰代表了枚舉了所有 5 個 {0, 1} 可以組成的字串。

De Bruijn 數列 - byte shift 版

De Bruijn 序列一開始的精神是 「rotation 1, 2 … 31 次之後, 就湊滿 0~31」。 那麼可不可以不要 rortaion 31 次, 而是找「 byte shift 0, 1, …31 次之後， 剛好湊滿 0~31 ?」 答案是可以的。 實際上這跟上面的 De Bruijn 是同一個問題 ： 只要找到的 De Bruijn 序列「前 5 個 bit 是 0」， 那做 rotation 取前 5 個bit, 跟做 byte shift 是一樣的。 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] 當中， 而方法就是：

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

不過 hash 順序不一定是照位元數目。 這樣其實也可以說是「把 byte 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];

1. (v & -v) ： 會把「最低位的1」 過濾出來, 除了他以外的 bit 通通都會變成 0。 原理是計概學過的「二的補數」的換算方法, 右邊開始往左碰到第一個不是0的數之後, 後面全部變號。 所以最低位的1往左的位元全部相反, 右邊又都是0, AND 起來通通消光光

2. 把(v & -v)乘上上面的數(0x07C4ACDD) : 因為這個東西恰好是

   $2^N$所以恰好在做左移。 剛剛知道這個數列的特性就是「對給定的平移, 可以在 0~31 中找到唯一的結果」。 所以只要建個大小 32 的表， 把平移量 map 到的欄位， 擺上對應的位數， 就可以實作一個查表的「count trailing zero」了。

3. 這支程式就在做這件事。 來源是這裡

De Bruijin 數列 - clz 版

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

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

​​​​h(a) = (a * M) >> 27

之後， 把所有形如下面的 a :

00000000000000000000000000000001
00000000000000000000000000000011
00000000000000000000000000000111
00000000000000000000000000001111
00000000000000000000000000011111
00000000000000000000000000111111
00000000000000000000000001111111
00000000000000000000000011111111
00000000000000000000000111111111
00000000000000000000001111111111

...

01111111111111111111111111111111
11111111111111111111111111111111

hash 進 [0, 31] 當中 ？」

先解釋一下為什麼要是左邊一排 0, 右邊一排 1 的形式 , 因為對於任意 32 bit 整數由下面這個步驟

    a |= value >> 1;
    a |= value >> 2;
    a |= value >> 4;
    a |= value >> 8;
    a |= 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;

接下來有兩種選擇：

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

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 ：

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 這種形式的數，也是可以讓性質成立的。

另外，也不一定要把 32 個這種形式的數 hash 進恰好大小為 32 的陣列，只要是合理的 hash 即可。比如說劉亮谷學長的實作

unsigned clz(uint32_t x)
{
    static uint8_t const Table[] = { 
      0xFF, 0, 0xFF, 15, 0xFF, 1, 28, 0xFF,
      16, 0xFF, 0xFF, 0xFF, 2, 21, 29, 0xFF,
      0xFF, 0xFF, 19, 17, 10, 0xFF, 12, 0xFF,
      0xFF, 3, 0xFF, 6, 0xFF, 22, 30, 0xFF,
      14, 0xFF, 27, 0xFF, 0xFF, 0xFF, 20, 0xFF,
      18, 9, 11, 0xFF, 5, 0xFF, 0xFF, 13, 
      26, 0xFF, 0xFF, 8, 0xFF, 4, 0xFF, 25, 
      0xFF, 7, 24, 0xFF, 23, 0xFF, 31, 0xFF,
    };  

    /* Propagate leftmost 1-bit to the right */
    x = x | (x >> 1); 
    x = x | (x >> 2); 
    x = x | (x >> 4); 
    x = x | (x >> 8); 
    x = x | (x >> 16);
 
    /* x = x * 0x6EB14F9 */
    x = (x << 3) - x;   /* Multiply by 7. */
    x = (x << 8) - x;   /* Multiply by 255. */
    x = (x << 8) - x;   /* Again. */
    x = (x << 8) - x;   /* Again. */

    return Table[x >> 26];
}

就是用 0x6EB14F9 把 64 個 a 值 hash 進一個大小為 64 的陣列。不過這樣就會有些位置沒有使用到。

所以就把計算 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 一樣。

clz的應用

雖然我還沒深入了解clz, 但實際上第一次看到clz是用在這裡：

「Efficient and easy segment trees」

利用陣列模擬線段樹，並且利用 clz 計算樹高。不過當初並不知道為什麼這個實作是正確的，所以並沒有使用它。

常用的實作是老老實實的把每一個節點刻出來，而且用遞迴版的深度優先搜尋遍歷。很常超時。

檢查程式正確性

看了一下裡面 code, 發現有好的東西 ：

#if defined(correct)

    for (int try = 0; try < 20; try++) {
                    timec = 0;
                    get_cycles(&timec_high1, &timec_low1);
                    for (uint32_t i = 0; i < 31; i++) {
                        printf("%u:%d \n",1<<i,clz(1<<i));
                        for (uint32_t j = (1 << i); j < (1 << (i + 1)); j++) {
                            assert( __builtin_clz (j) == clz(j));
                        }
                    }
                    get_cycles_end(&timec_high2, &timec_low2);
                    timec = diff_in_cycles(timec_high1, timec_low1, timec_high2, timec_low2);
                    printf("executiom time : %lu cycles\n", timec);
                }
#else

所以就把 makefile 參數改成：

CFLAGS_common ?= -Wall -std=gnu99 -g -Dcorrect -DDEBUG

然後放著等 assert 爆掉

後來他跑了很久。 再看了一眼發現可以把 for 迴圈弄掉, 不然每個都會跑 25 次。

跑完之後發現沒有因為 assert 失敗而中止。 然後再看了一眼， 發現當 i = 30 時, j 的範圍為 [1 << 30, 1 << 31), 但是這樣一來， (1 << 31) + 1 這個數字就沒有檢驗到了。 不知道是不是因為這些數字 clz 都是 0 的關係？

原始實驗結果

本來以為最快的 harley 居然最慢？感覺是大小 64 的陣列惹的禍。也許可以看一下 cache miss 的情況。

_Generic 修改

稍微看了一下， 發現裡面有 .hpp 與 .cc ， 打開一看發現有不同型別的實作。 嘗試把它用 _Generetic 改寫 ：

#define clz(T) _Generic((T), \
    uint8_t : clz8,\
    uint16_t : clz16, \
    default : clz32, \
)(T)

然後發現 cc 裡面的實作似乎有點問題 ：

unsigned clz8(uint8_t x)
{
    uint8_t upper = (uint8_t)(x >> 4);
    uint8_t lower = (uint8_t)(x & 0xF);
    return upper ? magic[upper] : 8 + magic[lower];
}

如果帶 0x1 進去的話看起來不太對。 先塞個陽春的版本：

unsigned clz8(uint8_t x)
{
    if (x == 0) return 8;
    int n = 0;
    if (x <= 0x0F) { n += 4; x <<= 4; }
    if (x <= 0x3F) { n += 2; x <<= 2; }
    if (x <= 0x7F) { n += 1; x <<= 1; }
    return n;
}

然後仿照前面的 code 去跑測試。 assert 均有通過。

覺得最有趣的地方是這個 _Generic 寫成的函數(?)遞迴下去的時候會在不同型別的實作間跳來跳去。