# [2020q3](http://wiki.csie.ncku.edu.tw/sysprog/schedule) 第 3 週測驗題
###### tags: `sysprog2020`

:::info
目的: 檢驗學員對 **[數值系統](https://hackmd.io/@sysprog/c-numerics)** 及 [bitwise 操作](https://hackmd.io/@sysprog/c-bitwise) 的認知
:::

==[作答表單](https://docs.google.com/forms/d/e/1FAIpQLSfMPyH-05j5BBCgqkmeOMZ5Xzw-usNeBcv050xMtfegQ9SciA/viewform)==

### 測驗 `1`

依據 [ISO/IEC 9899:TC2](http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1124.pdf) (即 C99) 標準的 6.5.7 章節 (第 84 到第 85 頁):
> "The result of `E1 >> E2` is E1 right-shifted E2 bit positions. ... If E1 has a signed type and a negative value, the resulting value is ==implementation-defined==."

在 C 語言中，對有號型態的物件進行算術右移 (arithmetic right shift，以下簡稱 ASR 或 asr) 歸屬於 [unspecified behavior](https://hackmd.io/@sysprog/c-undefined-behavior)，包含 [Cppcheck](http://cppcheck.sourceforge.net/) 在內的靜態分析器會嘗試找出 C 程式中這項 [unspecified behavior](https://hackmd.io/@sysprog/c-undefined-behavior)。考慮到程式碼的相容性，我們想實作一套可跨越平台的 ASR，假設在[二補數系統](https://hackmd.io/@sysprog/binary-representation)，無論標的環境的 shift 究竟輸出為何，我們的 ASR 總可做到:
$-5 \gg^{arith} 1 = -3$
上述 $-3$ 正是 $\frac{-5}{2}$ 的輸出，並向 $-\infty$ 進位 (rounding)。

針對有號整數的跨平台 ASR 實作如下:

int asr_i(signed int m, unsigned int n)
{
    const int logical = (((int) -1) OP1) > 0;
    unsigned int fixu = -(logical & (OP2));
    int fix = *(int *) &fixu;
    return (m >> n) | (fix ^ (fix >> n));
}

請補完程式碼。

==作答區==

OP1 = ?
* `(a)` `<< 1`
* `(b)` `>> 1`
* `(c)` `* m + n`
* `(d)` `+ m`
* `(e)` `+ n`
* `(f)` `- m`
* `(g)` `- n`

OP2 = ?
* `(a)` `m`
* `(b)` `n`
* `(c)` `m < 0`
* `(d)` `m == 0`
* `(e)` `m > 0`
* `(f)` `n < 0`
* `(g)` `n == 0`
* `(h)` `n > 0`

提示: 透過 gcc/clang 編譯程式碼時，可加上 `-fsanitize=undefined` 編譯參數，在執行時期若遇到[未定義行為](https://hackmd.io/@sysprog/c-undefined-behavior)，會得到以下錯誤訊息:
> runtime error: load of misaligned address 0x7ffd8a89be8f for type 'int', which requires 4 byte alignment

補充資訊:
* [Arm Assembly Language Programming: Arithmetic Shift Operations](http://www-mdp.eng.cam.ac.uk/web/library/enginfo/mdp_micro/lecture4/lecture4-3-3.html)

:::success
延伸問題:
1. 解釋上述程式運作原理，應提供數學證明和 Graphviz 製作的圖解;
2. 練習實作其他資料寬度的 ASR，可參照 [C 語言：前置處理器應用篇](https://hackmd.io/@sysprog/c-preprocessor) 撰寫通用的巨集以強化程式碼的共用程度;
:::

---

### 測驗 `2`

在 [C 語言：數值系統篇](https://hackmd.io/@sysprog/c-numerics) 中，我們介紹過 power of 2 (漢字可寫作「二的冪」)。
> 女星楊冪的名字裡，也有一個「冪」字。且，她取這個名字，就有「次方」的意義，因為她一家三口都姓楊，所以是「楊的三次方」。

[GCC extension](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html) 其中兩個是 ctz 和 clz:
> Built-in Function: `int __builtin_ctz (unsigned int x)`
> * Returns the number of trailing 0-bits in x, starting at the least significant bit position. If x is 0, the result is undefined.
>
> Built-in Function: `int __builtin_clz (unsigned int x)`
> * Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.

我們可用 `__builtin_ctz` 來實作 LeetCode [342. Power of Four](https://leetcode.com/problems/power-of-four/)，假設 `int` 為 32 位元寬度。實作程式碼如下:

bool isPowerOfFour(int num)
{
    return num > 0 && (num & (num - 1))==0 &&
           !(__builtin_ctz(num) OPQ);  
}

請補完程式碼。

==作答區==

OPQ = ?
* `(a)` `> 0`
* `(b)` `> 31`
* `(c)` `>> 0x2`
* `(d)` `>> 0x1`
* `(e)` `| 0x1`
* `(f)` `& 0x1`
* `(g)` `^ 0x1`

:::success
延伸問題:
1. 解釋上述程式運作原理;
2. 改寫上述程式碼，提供等價功能的不同實作，儘量降低 branch 的數量;
3. 練習 LeetCode [1009. Complement of Base 10 Integer](https://leetcode.com/problems/complement-of-base-10-integer/) 和 [41. First Missing Positive](https://leetcode.com/problems/first-missing-positive/)，應善用 clz;
4. 研讀 [2017 年修課學生報告](https://hackmd.io/@3xOSPTI6QMGdj6jgMMe08w/Bk-uxCYxz)，理解 `clz` 的實作方式，並舉出 [Exponential Golomb coding](https://en.wikipedia.org/wiki/Exponential-Golomb_coding) 的案例說明，需要有對應的 C 原始程式碼和測試報告;
> [x-compressor](https://github.com/jserv/x-compressor) 是個以 [Golomb-Rice coding](https://en.wikipedia.org/wiki/Golomb_coding) 為基礎的資料壓縮器，實作中也用到 clz/ctz 指令，可參見 [Selecting the Golomb Parameter in Rice Coding](https://ipnpr.jpl.nasa.gov/progress_report/42-159/159E.pdf)。
:::

---

### 測驗 `3`

[population count](https://en.wikichip.org/wiki/population_count) 簡稱 popcount 或叫 sideways sum，是計算數值的二進位表示中，有多少位元是 `1`，在一些場合下很有用，例如計算 0-1 稀疏矩陣 (sparse matrix)或 bit array 中非 `0` 元素個數、計算兩個字串的 [Hamming distance](https://en.wikipedia.org/wiki/Hamming_weight)。Intel 在 2008 年 11 月 Nehalem 架構的處理器 Core i7 引入 SSE4.2 指令集，其中就有 `CRC32` 和 `POPCNT` 指令，`POPCNT` 可處理 16-bit, 32-bit, 64-bit 整數。

對應到 C 程式的實作:

unsigned popcnt_naive(unsigned v) {
    unsigned n = 0;
    while (v)
        v &= (v - 1), n = -(~n);
    return n;
}

函式 `popcnt_naive()` 利用不斷清除 LSB 直到輸入的數值 `v` 為 0。

`v &= (v - 1)` 的作用是將 LSB 設為 0 (**reset LSB**) 舉例來說，假設輸入數值為 `20`:

0001 0100  # 20          ; LSB in bit position 2
0001 0011  # 20 - 1
0001 0000  # 20 & (20 - 1)

> 類似的操作還有 `x & -x`，將 `x` 的 LSB 取出 (**isolate LSB**)

`n = -(~n)` 等同 `n++`，因為在二補數系統中，
$-n =\ \sim n + 1$
$-(\sim n) = n + 1$

因此 `popcnt_naive()` 的執行時間取決於輸入數值中 1 (set bit) 的個數。可改寫為以下常數時間複雜度的實作:

unsigned popcnt_branchless(unsigned v)
{
    unsigned n;
    n = (v >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;

    v = (v + (v >> 4)) & 0x0F0F0F0F;
    v *= 0x01010101;                                     

    return v >> 24;
}

對於一個 32 bit 的無號整數，popcount 可以寫成以下數學式:

$popcnt(x) = x - \left \lfloor{{\dfrac{x}{2}}}\right \rfloor - \left \lfloor{{\dfrac{x}{4}}}\right \rfloor - ... - \left \lfloor{{\dfrac{x}{2^{{31}}}}}\right \rfloor$

假設 $x = b_{31}...b_3b_2b_1b_0$，先看看 $x[3:0]$ 4 個位元，用以上公式可以計算得:

$(2^3b_3 + 2^2b_2 + 2^1b_1 + 2^0b_0) -
(2^2b_3 + 2^1b_2 + 2^0b_1) - (2^1b_3 + 2^0b_2) - 2^0b_3$
> $\left \lfloor{{\dfrac{x}{2}}}\right \rfloor$ 相當於 C 表達式中 `x >> 1`

稍微改寫可得到:
$(2^3 - 2^2 - 2^1 - 2^0)b_3 + (2^2 - 2^1 - 2^0)b_2 + (2^1 - 2^0)b_1 + 2^0b_0$

因此 popcnt 的一般式可改寫為:
$popcnt(x) = \sum\limits_{n=0}^{31} {}(2^n - \sum\limits_{i=0}^{n-1} 2^{i})b_n = \sum\limits_{n=0}^{31}b_n$

因為 $2^n - \sum\limits_{i=0}^{n-1} 2^{i} = 1$，只要對應的 $b_n$ 為 1，這個 bit 就會在 popcnt 的總和中加一，剛好對應 `popcnt_naive()`，因此映證上述數學式確實可計算出 population count。

且一個 32 位元無號整數最多有 32 個 1 (set bit)，剛好可用一個 byte 表示，所以可分成幾個區塊平行計算，最後再全部加總到一個 byte 中，進行避免檢查 32 次。

`popcnt_branchless` 實作一開始以**每 4 個位元 (nibble) 為一個單位**計算 1 的個數，利用最初的公式計算 $x - \left \lfloor{{\dfrac{x}{2}}}\right \rfloor - \left \lfloor{{\dfrac{x}{4}}}\right \rfloor - \left \lfloor{{\dfrac{x}{8}}}\right \rfloor$

關鍵的程式碼，逐行解釋:

n = (v >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;

1. `n = (v >> 1) & 0x77777777` : 將輸入數值 `v` 除以 2，得到 $\left \lfloor{{\dfrac{v}{2}}}\right \rfloor$

    b_31 b_30 b_29 b_28 ... b7 b6 b5 b4 b3 b2 b1 b0  // v
       0 b_31 b_30 b_29 ...  0 b7 b6 b4  0 b3 b2 b1  // (v >> 1) & 0x77777777

2. `v -= n` : 計算結果相當於 $v - \left \lfloor{{\dfrac{v}{2}}}\right \rfloor$
3. `n = (n >> 1) & 0x77777777` : 再對 `n` 除以 2，得到 $\left \lfloor{{\dfrac{v}{4}}}\right \rfloor$
4. `v -= n` : 計算出 $v - \left \lfloor{{\dfrac{v}{2}}}\right \rfloor - \left \lfloor{{\dfrac{v}{4}}}\right \rfloor$
5. 和 6. 重複同樣動作

最後這段結束後計算出 $v - \left \lfloor{{\dfrac{v}{2}}}\right \rfloor - \left \lfloor{{\dfrac{v}{4}}}\right \rfloor - \left \lfloor{{\dfrac{v}{8}}}\right \rfloor$，得到每 4 個位元為一個單位中 set bit 的個數

**`v = (v + (v >> 4)) & 0x0F0F0F0F`** : 將每 4 個位元中 set bit 的個數加到 byte 中:
1. 假設 $B_n$ 代表第 n 個 nibble (4 位元) 中的數值

    B7 B6 B5 B4 B3 B2 B1 B0  // v
     0 B7 B6 B5 B4 B3 B2 B1  // (v >> 4)

2. 加總可得到:

    // (v + (v >> 4))
    B7 (B7+B6) (B6+B5) (B5+B4) (B4+B3) (B3+B2) (B2+B1) (B1+B0)

3. 最後使用 `0x0F0F0F0F` 做 mask 可得

    // (v + (v >> 4)) & 0x0F0F0F0F
    0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)

**`v *= 0x01010101`** : 在最後一道敘述中，將 `v` 乘上 `0x01010101`。寫成直式如下

                         0 A6  0 A4  0 A2  0 A0
                      x  0  1  0  1  0  1  0  1
---------------------------------------------------
                         0 A6  0 A4  0 A2  0 A0
                   0 A6  0 A4  0 A2  0 A0  0
             0 A6  0 A4  0 A2  0 A0  0
       0 A6  0 A4  0 A2  0 A0  0
---------------------------------------------------
                            ↑_______________________A6+A4+A2+A0

我們可發現期望輸出就在原本 $A_6$ 的位置 ($2^7$)，因此將 `v` 右移 24 bits 即為所求，剩下的位數會 overflow ，右移後不影響結果。
* 假設 $A = B_7 + B_6$, $B = B_5 + B_4$, $C = B_3 + B_2$, $D = B_1 + B_0$, 根據分配律可得:

    v * 0x01010101 = 
    (A + B + C + D)  (B + C + D)      (C + D)          (D)
    |<-- 1 byte -->|<-- 1 byte -->|<-- 1 byte -->|<-- 1 byte -->|

**`return v >> 24`** : 最後得到的結果會放在 Most significant byte 中，因此向右位移 24 位元，即為所求的 popcount 值。

GCC 提供對應的內建函式:
> `__builtin_popcount(x)`: 計算 x 的二進位表示中，總共有幾個 `1`

使用示範：

int x = 5328; // 00000000000000000001010011010000
printf("%d\n",  __builtin_popcount(x)); // 5

針對 LeetCode [1342. Number of Steps to Reduce a Number to Zero](https://leetcode.com/problems/number-of-steps-to-reduce-a-number-to-zero/)，我們可運用 gcc 內建函式實作，假設 `int` 寬度為 32 位元，程式碼列表如下:

int numberOfSteps (int num)
{
    return num ? AAA + 31 - BBB : 0;
}

請補完程式碼。

==作答區==

AAA = ?
* `(a)` `__builtin_popcount(num)`
* `(b)` `__builtin_clz(num)`

BBB = ?
* `(a)` `__builtin_popcount(num)`
* `(b)` `__builtin_clz(num)`

:::success
延伸問題:
1. 解釋上述程式運作原理;
2. 改寫上述程式碼，提供等價功能的不同實作並解說;
> 提示: [Bit Twiddling Hacks](http://graphics.stanford.edu/~seander/bithacks.html) 中提及許多 bitwise operation 的使用，如 [bit inverse](https://hackmd.io/@sysprog/ByzoiggIb)、abs 等等，是極好的參考材料。
3. 避免用到編譯器的擴充功能，只用 C99 特徵及 bitwise 運算改寫 LeetCode [1342. Number of Steps to Reduce a Number to Zero](https://leetcode.com/problems/number-of-steps-to-reduce-a-number-to-zero/)，實作出 branchless 解法，儘量縮減程式碼行數和分支數量;
:::

---

### 測驗 `4`

考慮以下 64-bit GCD (greatest common divisor, 最大公因數) 求值函式:

#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v) {
    if (!u || !v) return u | v;
    while (v) {                               
        uint64_t t = v;
        v = u % v;
        u = t;
    }
    return u;
}

改寫為以下等價實作:

#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v) {
    if (!u || !v) return u | v;
    int shift;
    for (shift = 0; !((u | v) & 1); shift++) {
        u /= 2, v /= 2;
    }
    while (!(u & 1))
        u /= 2;
    do {
        while (!(v & 1))
            v /= 2;
        if (u < v) {
            v -= u;
        } else {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while (XXX);
    return YYY;
}

請補完程式碼。

==作答區==

XXX = ?
* `(a)` `u`
* `(b)` `v`
* `(c)` `u + v`
* `(d)` `u - v`
* `(e)` `u >> 1`
* `(f)` `v >> 1`

YYY = ?
* `(a)` `u`
* `(b)` `v`
* `(c)` `u << v`
* `(d)` `v << u`
* `(e)` `u << shift`
* `(f)` `v << shift`

:::success
延伸問題:
1. 解釋上述程式運作原理;
2. 在 x86_64 上透過 `__builtin_ctz` 改寫 GCD，分析對效能的提升;
:::

---

### 測驗 `5`

在影像處理中，[bit array](https://en.wikipedia.org/wiki/Bit_array) (也稱 `bitset`) 廣泛使用，考慮以下程式碼:

#include <stddef.h>
size_t naive(uint64_t *bitmap, size_t bitmapsize, uint32_t *out)
{
    size_t pos = 0;
    for (size_t k = 0; k < bitmapsize; ++k) {
        uint64_t bitset = bitmap[k];
        size_t p = k * 64;
        for (int i = 0; i < 64; i++) {
            if ((bitset >> i) & 0x1)
                out[pos++] = p + i;
        }
    }
    return pos;
}

可用 clz 改寫為效率更高且等價的程式碼:

#include <stddef.h>
size_t improved(uint64_t *bitmap, size_t bitmapsize, uint32_t *out)
{
    size_t pos = 0;
    uint64_t bitset;
    for (size_t k = 0; k < bitmapsize; ++k) {
        bitset = bitmap[k];
        while (bitset != 0) {
            uint64_t t = KKK;
            int r = __builtin_ctzll(bitset);
            out[pos++] = k * 64 + r;
            bitset ^= t;                           
        }
    }
    return pos;
}

請補完程式碼。

==作答區==

KKK = ?
* `(a)` `bitset`
* `(b)` `bitset & 0x1`
* `(c)` `bitset | 0x1`
* `(d)` `bitset ^ 0x1`
* `(e)` `bitset & -bitset`

:::success
延伸問題:
1. 解釋上述程式運作原理，並舉出這樣的程式碼用在哪些真實案例中;
2. 設計實驗，檢驗 clz 改寫的程式碼相較原本的實作有多少改進？應考慮到不同的 [bitmap density](http://www.vki.com/2013/Support/docs/vgltools-3.html);
3. 思考進一步的改進空間;
:::

---