# 2020q3 Homework (quiz4)
contributed by < `hankluo6` >

### 測驗 1

int hammingDistance(int x, int y) {
    return __builtin_popcount(x OP y);
}

#### `OP`

- [x] `(c)` `^`

使用 XOR 會將兩個有差異的 bit 設為 1，相同的設為 0，所以回傳的值就是 hamming distance。

#### 延伸問題
* 舉出不使用 gcc extension 的 C99 實作程式碼

    int hammingDistance(int x, int y) {
        int v = x ^ y;
        v = v - ((v >> 1) & 0x55555555);
        v = (v & 0x33333333) + ((v >> 2) & 0x33333333);
        return ((v + (v >> 4) & 0xF0F0F0F) * 0x1010101) >> 24;
    }

參考 [Bit Twiddling Hacks](http://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel)，模擬 `popcount()` 的行為。
* 練習 LeetCode [477. Total Hamming Distance](https://leetcode.com/problems/total-hamming-distance/)，你應當提交 (submit) 你的實作 (C 語言實作)，並確保執行結果正確且不超時

一開始使用兩個 for 迴圈，但是會有 TLE 的問題，改成對每個數字的 bit 著手，這樣可以確保經過 32 次迴圈後一定能結束。觀察在 `numsSize` 數字內的同一個 bit 造成的 hamming distance 為 **1 的數量 * 0 的數量**，這是因為每個 1 與每個 0 都可以互相配對，此種方法為將原本需要搜尋所有可能的 pair，移除兩個 bit 相同的情形，進而減少程式碼。

    int totalHammingDistance(int* nums, int numsSize){
        int total = 0;
        int n;
        int count;
        for (;;) {
            n = 0;
            count = 0;
            for (int j = 0; j < numsSize; ++j) {
                if (nums[j] & 1)
                    ++n;
                if (!nums[j])
                    ++count;
                nums[j] >>= 1;
            }
            if (count == numsSize)
                break;
            total += (numsSize - n) * n;
        }
        return total;
    }

這邊使用額外一個變數來讓迴圈提早結束，此程式可過 leetcode 測試：
>Runtime: **44 ms**, faster than **71.56%** of C online submissions for Total Hamming Distance.
Memory Usage: **6.7 MB**, less than **37.62%** of C online submissions for Total Hamming Distance.

看起來還能再改進。
* 研讀[錯誤更正碼簡介](https://w3.math.sinica.edu.tw/math_media/d184/18404.pdf)及[抽象代數的實務應用](https://www.math.sinica.edu.tw/www/file_upload/summer/crypt2017/data/2017/%E4%B8%8A%E8%AA%B2%E8%AC%9B%E7%BE%A9/[20170710][%E6%8A%BD%E8%B1%A1%E4%BB%A3%E6%95%B8%E7%9A%84%E5%AF%A6%E5%8B%99%E6%87%89%E7%94%A8].pdf)，摘錄其中 Hamming Distance 和 Hamming Weight 的描述並重新闡述，舉出實際的 C 程式碼來解說;
* 搭配閱讀 [Hamming codes, Part I]((https://www.youtube.com/watch?v=X8jsijhllIA&feature=youtu.be&ab_channel=3Blue1Brown)) 及 [Hamming codes, Part II](https://www.youtube.com/watch?v=b3NxrZOu_CE&feature=youtu.be)，你可適度摘錄影片中的素材，但應當清楚標註來源
#### (7,4) Hamming Code
假設我們要傳入的為 4 bit 之資料，而我們會將原始資料再加上 3 bit 的 parity bit，形成 7 bit 的資料來傳輸。

前 4 個 bit 可以對應到下圖 set 當中交集的 1, 2, 3, 4 部分，而 剩餘的 3 個 bit 則是 set 外圍的部分。
parity bit 的選擇是依據每個圓圈裡面 1 的個數，只要圓圈內 1 的個數為奇數，則將該圓對應的 parity bit 設為 1，反之則設為 0，目的是讓每個圓圈裡面 1 的個數都為偶數。

接著我們就能將上述規則寫成代數形式，這邊加法為不含 carry 的二進位加法，也可以當成 XOR 運算子：
$$
\left\{
\begin{array}{c}
x1 + x2 + x3 + x5 = 0 \\
x1 + x2 + x4 + x6 = 0 \\
x1 + x3 + x4 + x7 = 0 \\
\end{array}
\right.
$$
將此等式寫成矩陣形式 $H$：
$$
H=\begin{pmatrix}
1 & 1 & 1 & 0 & 1 & 0 & 0 \\
1 & 1 & 0 & 1 & 0 & 1 & 0 \\
1 & 0 & 1 & 1 & 0 & 0 & 1\\
\end{pmatrix}
$$
因為 $H\vec{x} = 0$，我們要求的矩陣 $C=(x1, x2,..., x7)$ 就為矩陣 $H$ 的 null space。而 $rank(H)=3,\ rank(C) = 7 - rank(H) = 4$，所以 $C$ 可以選擇 4 個變數，而剩餘 3 個則可以透過這些變數來決定。

我們只需透過矩陣乘法就能找出錯誤的 bit，假設傳出的資料 $\vec{y}=\vec{x}+\vec{e}$，其中 $\vec{x}$ 為正確的資料，$\vec{e}$ 為錯誤向量，所以可以寫成 $H\vec{y}=H\vec{x}+H\vec{e}=H\vec{e}$，只要相乘的結果不為 $\vec{0}$，則表示有錯誤。

在習慣上，我們可以將 parity bit 放置在 2 的次方的位置上，這樣矩陣相乘的結果向量就是錯誤 bit 的位置，影片參考 [Hamming codes part 2](https://youtu.be/b3NxrZOu_CE?t=392)。

除了利用矩陣外，hamming code 其實可以用 XOR 來達成上面的運算，因為 parity bit 用來判斷奇偶數個 1，可以用對應的 bit 做 XOR 來取得，如果將 parity bit 放置在 2 的指數次方位置，parity bit 檢測的範圍為特定的 bit 位置，如 (7,4) hamming code 中：
* 001 位置的 parity bit 用來檢查 011, 101, 111 位置 1 的個數
* 010 位置的 parity bit 用來檢查 011, 110, 111 位置 1 的個數
* 100 位置的 parity bit 用來檢查 101, 110, 111 位置 1 的個數

其他資料的 bit 按順序放入剩餘的 bit，這樣檢查時只需對每個為 1 的 bit 的二進位位置做 XOR，結果就為錯誤 bit 的位置。這是因為錯誤的位置會剛好影響到對應的 parity bit，這樣可以有效的減少實作的程式碼。

    void decoder(int receive)
    {
        int error = 0;
        for (int i = 0; i < 32; ++i) {
            if (receive & (1 << i)) {
                error ^= i;
            }
        }
        if (!error) {
            receive ^= (1 << error);
        }
        return receive;
    }

使用 32 bit integer 來儲存 32 bit 的資料，可以看到程式碼相當簡潔，只需將每個為 1 的 bit 的位置做 XOR，其結果就是 error 的位置。

* 研讀 [Reed–Solomon error correction](https://en.wikipedia.org/wiki/Reed%E2%80%93Solomon_error_correction)，再閱讀 Linux 核心原始程式碼的 [lib/reed_solomon 目錄](https://github.com/torvalds/linux/tree/master/lib/reed_solomon)，抽離後者程式碼為單獨可執行的程式，作為 ECC 的示範

### 測驗 2

#include <stdlib.h>

typedef struct {
    AAA;
    int max_level;
} TreeAncestor;

TreeAncestor *treeAncestorCreate(int n, int *parent, int parentSize)
{
    TreeAncestor *obj = malloc(sizeof(TreeAncestor));
    obj->parent = malloc(n * sizeof(int *));

    int max_level = 32 - __builtin_clz(n) + 1;
    for (int i = 0; i < n; i++) {
        obj->parent[i] = malloc(max_level * sizeof(int));
        for (int j = 0; j < max_level; j++)
            obj->parent[i][j] = -1;
    }
    for (int i = 0; i < parentSize; i++)
        obj->parent[i][0] = parent[i];

    for (int j = 1;; j++) {
        int quit = 1;
        for (int i = 0; i < parentSize; i++) {
            obj->parent[i][j] = obj->parent[i][j + BBB] == -1
                                    ? -1
                                    : obj->parent[obj->parent[i][j - 1]][j - 1];
            if (obj->parent[i][j] != -1) quit = 0;
        }
        if (quit) break;
    }
    obj->max_level = max_level - 1;
    return obj;
}

int treeAncestorGetKthAncestor(TreeAncestor *obj, int node, int k)
{
    int max_level = obj->max_level;
    for (int i = 0; i < max_level && node != -1; ++i)
        if (k & (CCC))
            node = obj->parent[node][i];
    return node;
}

void treeAncestorFree(TreeAncestor *obj)
{
    for (int i = 0; i < obj->n; i++)
        free(obj->parent[i]);
    free(obj->parent);
    free(obj);
}

#### `AAA`

- [x] `(b)` `int **parent`

#### `BBB`

- [x] `(b)` `(-1)`

#### `CCC`

- [x] `(f)` `1 << i`

`AAA` 是用來記憶每個 node i 的第 k 個 parent，可看到下方程式碼使用到二為陣列，所以這邊應為 `**parent`。

使用 dynamic programming 的技巧、Top down 的方式來解題。從樹的 root 節點開始向下，因每個 i 節點的第 j 個 parent 就是其 parent 的第 j - 1 個 parent，所以更新公式可以寫成 `obj->parent[i][j] = obj->parent[obj->parent[i][j - 1]][j - 1];`，但還要包括邊界條件，當第 j 個 parent 的 parent 沒有時，設定為 -1。

所以 `BBB` 是用來檢查是否到達邊界，故為 -1。

`CCC` 用來判斷 k 的 bit，運用每個 bit 來取得對應數字的 parent，再透過迴圈 loop k 的 bit，最後就能取得到第 k 個 parent node，所以 `CCC` 為 `1 << i`。

#### 延伸問題

* 指出可改進的地方，並實作對應的程式碼;
* `max_level` 在實際上多分配一個 node 的空間，因為下方 `for (int j = 1;; j++)` 是用此位置來判斷是否結束，實在是多此一舉，將 `max_level` 改為真正的高度並設定迴圈條件 `for (int j = 1; j < max_level; j++)`。
* 在 `treeAncestorFree()` 中有用到 `obj->n`，但結構中沒有這項成員，將之加入結構中並新增 `obj->n = n` 在 `treeAncestorCreate()` 當中。
* `Create` 中第三個迴圈是以 row 為主，可能會有 locality 的問題，將陣列紀錄順序對調。
* 在 treeAncestorCreate 函式內部，若干記憶體被浪費，請撰寫完整測試程式碼，並透過工具分析;
* 在 LeetCode [1483. Kth Ancestor of a Tree Node](https://leetcode.com/problems/kth-ancestor-of-a-tree-node/) 提交你的實作，你應該要在 C 語言項目中，執行時間領先 75% 的提交者;

將上述內容實作後並提交：

    TreeAncestor *treeAncestorCreate(int n, int *parent, int parentSize)
    {
        TreeAncestor *obj = malloc(sizeof(TreeAncestor));
        int max_level = 32 - __builtin_clz(n);
        obj->parent = malloc(max_level * sizeof(int *));

        for (int i = 0; i < max_level; i++) {
            obj->parent[i] = malloc(n * sizeof(int));
            for (int j = 0; j < n; j++)
                obj->parent[i][j] = -1;
        }
        for (int i = 0; i < parentSize; i++)
            obj->parent[0][i] = parent[i];

        for (int i = 1; i < max_level; i++) {
            for (int j = 0; j < parentSize; j++) {
                obj->parent[i][j] = obj->parent[i - 1][j] == -1
                                        ? -1
                                        : obj->parent[i - 1][obj->parent[i - 1][j]];
            }
        }
        obj->max_level = max_level;
        obj->n = n;
        return obj;
    }

Runtime: **244 ms**, faster than **98.45%** of C online submissions for Kth Ancestor of a Tree Node.
Memory Usage: **62.9 MB**, less than **6.50%** of C online submissions for Kth Ancestor of a Tree Node.

### 測驗 3

#define MSG_LEN 8

static inline bool is_divisible(uint32_t n, uint64_t M) {
    return n * M <= M - 1;
}

static uint64_t M3 = UINT64_C(0xFFFFFFFFFFFFFFFF) / 3 + 1;
static uint64_t M5 = UINT64_C(0xFFFFFFFFFFFFFFFF) / 5 + 1;

int main(int argc, char *argv[]) {
    for (size_t i = 1; i <= 100; i++) {
        uint8_t div3 = is_divisible(i, M3);
        uint8_t div5 = is_divisible(i, M5);
        unsigned int length = (2 << div3) << div5;

        char fmt[MSG_LEN + 1];
        strncpy(fmt, &"FizzBuzz%u"[(MSG_LEN >> KK1) >> (KK2 << KK3)], length);
        fmt[length] = '\0';

        printf(fmt, i);
        printf("\n");
    }
    return 0;
}

:::info
`is_divisible` 可以參考 quiz2 的 [第三題](https://hackmd.io/85fHzF42SOasoMatafgZkQ?view#%E6%B8%AC%E9%A9%97-3) 與 [第四題](https://hackmd.io/85fHzF42SOasoMatafgZkQ?view#%E6%B8%AC%E9%A9%97-4)。
:::

根據題意：我們若能精準從給定輸入 i 的規律去控制 start 及 length，即可符合 FizzBuzz 題意，可以把 div3, div5 與對應的值寫下：

| div3 | div5 | (MSG_LEN >> KK1) >> (KK2 << KK3) |
|:----:|:----:|:--------------------------------:|
| 0 | 0 | 8 |
| 1 | 0 | 0 |
| 0 | 1 | 4 |
| 1 | 1 | 0 |

分析：

可以看到對 `MSG_LEN` 的運算元只有 `>>`，所以在 `div3` 與 `div5` 皆為 0 的情形時，可得到 `(KK1 == 0) && ((KK2 << KK3) == 0)` 的關係。

在 `div3` 為 0 且 `div5` 為 1 時，`MSG_LEN` 應要右移一個 bit，此有兩種可能：
1. KK1 = 1, KK2 = 0, KK3 = ?
2. KK1 = 0, KK2 = 1, KK3 = 0

在 `div3` 為 1 且 `div5` 為 0 時，`MSG_LEN` 應要右移四個 bit 以上，此有七種可能：
1. KK1 = 3, KK2 = 0, KK3 = ?
2. KK1 = 2, KK2 = 1, KK3 = 1
3. KK1 = 2, KK2 = 2, KK3 = 0
4. KK1 = 1, KK2 = 3, KK3 = 0
5. KK1 = 0, KK2 = 1, KK3 = 2
6. KK1 = 0, KK2 = 2, KK3 = 1
7. KK1 = 0, KK2 = 4, KK3 = 0

在 `div3` 與 `div5` 皆為 1 時，與上方結果一樣。

帶數字就能得到結果。

#### `KK1`

- [x] `(a)` `div5`

#### `KK2`

- [x] `(d)` `div3`

#### `KK3`

- [x] `(c)` `2`

#### 延伸問題
* 評估 naive.c 和 bitwise.c 效能落差
* 避免 stream I/O 帶來的影響，可將 printf 更換為 sprintf

    #include <stdio.h>
    #include <time.h>
    #include <string.h>
    #include <stdint.h>
    #include <stdbool.h>
    
    #define MSG_LEN 8

    static inline bool is_divisible(uint32_t n, uint64_t M) {
        return n * M <= M - 1;
    }

    static uint64_t M3 = UINT64_C(0xFFFFFFFFFFFFFFFF) / 3 + 1;
    static uint64_t M5 = UINT64_C(0xFFFFFFFFFFFFFFFF) / 5 + 1;

    void naive() {
        for (unsigned int i = 1; i < 1000000; i++) {
            char s[MSG_LEN + 1];
            if (i % 3 == 0) 1sprintf(s, "Fizz");
            if (i % 5 == 0) sprintf(s, "Buzz");
            if (i % 15 == 0) sprintf(s, "FizzBuzz");
            if ((i % 3) && (i % 5)) sprintf(s, "%u", i);
            sprintf(s, "\n");
        }
    }

    void bitwise() {
        for (size_t i = 1; i <= 1000000; i++) {
            uint8_t div3 = is_divisible(i, M3);
            uint8_t div5 = is_divisible(i, M5);
            unsigned int length = (2 << div3) << div5;

            char fmt[MSG_LEN + 1];
            strncpy(fmt, &"FizzBuzz%u"[(MSG_LEN >> div5) >> (div3 << 2)], length);
            fmt[length] = '\0';
        }
    }

    int main(int argc, char *argv[]) {
        double time_taken;
        struct timespec start, end;
        clock_gettime(CLOCK_MONOTONIC, &start);
        bitwise();
        clock_gettime(CLOCK_MONOTONIC, &end);
        time_taken = (end.tv_sec - start.tv_sec) * 1e9;
        time_taken = (time_taken + (end.tv_nsec - start.tv_nsec)) * 1e-9;
        printf("%f,", time_taken);
        clock_gettime(CLOCK_MONOTONIC, &start);
        naive();
        clock_gettime(CLOCK_MONOTONIC, &end);
        time_taken = (end.tv_sec - start.tv_sec) * 1e9;
        time_taken = (time_taken + (end.tv_nsec - start.tv_nsec)) * 1e-9;
        printf("%f\n", time_taken);

        return 0;
    }


因 bitwise 少了分支與餘數運算，效能比 naive 好非常多。

* 分析 [Faster remainders when the divisor is a constant: beating compilers and libdivide](https://lemire.me/blog/2019/02/08/faster-remainders-when-the-divisor-is-a-constant-beating-compilers-and-libdivide/) 一文的想法 (可參照同一篇網誌下方的評論)，並設計另一種 bitmask，如「可被 3 整除則末位設為 1」「可被 5 整除則倒數第二位設定為 1」，然後再改寫 bitwise.c 程式碼，試圖運用更少的指令來實作出 branchless;
* 參照 [fastmod: A header file for fast 32-bit division remainders on 64-bit hardware](https://github.com/lemire/fastmod)

### 測驗 4

#include <stdint.h>
#define ffs32(x)  ((__builtin_ffs(x)) - 1)
size_t blockidx;
size_t divident = PAGE_QUEUES;
blockidx = offset >> ffs32(divident);
divident >>= ffs32(divident);
switch (divident) {
    CASES
}

#### `CASES` 至少該包含哪些數字: (複選)

- [x] `(b)` `3`
- [x] `(d)` `5`
- [x] `(f)` `7`

* `__builtin_ffs`
>Returns one plus the index of the least significant 1-bit of x, or if x is zero, returns zero.

回傳從右邊數來第一個 1 的位置

* `#define ffs32(x) ((__builtin_ffs(x)) - 1)`

回傳從右邊數來第一個 1 以前有幾個 0。
i.e. `ffs32() == __builtin_ctz()`

考慮到整數 PAGE_QUEUES 可能有以下數值:

* (1 or 2 or 3 or 4)
* (5 or 6 or 7 or 8) × ($2^n$), n from 1 to 14

可以寫成以下：

* (1 or 2 or 3 or 4)
* (5 or 7) × ($2^n$), n from 1 to 14
* (3 or 4) × ($2^n$), n from 2 to 15

可以把相同倍數的合併：

* (1 or 3 or 5 or 7) × ($2^n$), n from 0 to 15

這邊雖然某些數的 n 範圍不同，但我們可以擴充成 n 皆為相同，只要有涵蓋到原本的範圍就好。

blockidx = offset >> ffs32(divident);
divident >>= ffs32(divident);

程式碼中這兩步是把上面公式中的 $2^n$ 部份移除，剩下部份就要依靠 `CASE` 來處理：

switch (divident) {
    case 3: blockidx /= 3; 
              break;
    case 5: blockidx /= 5; 
              break;
    case 7: blockidx /= 7; 
              break;
}

因除以 1 沒有意義，所以不必寫在 `CASE` 中。

#### 延伸問題
* 解釋程式運算原理，可搭配 Microsoft Research 的專案 [snmalloc](https://github.com/microsoft/snmalloc) 來解說;
* 對應的原始程式碼 [src/mem/sizeclass.h](https://github.com/microsoft/snmalloc/blob/master/src/mem/sizeclass.h#L54-L145)

`snmalloc` 裡的 `round_by_sizeclass()` 用來計算 $\lfloor\frac{x}{y}\rfloor \times y$ 的值，而其中除法的部分使用 [fastdiv](https://hackmd.io/85fHzF42SOasoMatafgZkQ?view#%E6%B8%AC%E9%A9%97-3) 來計算，並與本題一樣利用 1, 3, 5, 7 當作 special case，便能依照對應的 case 加速運算。
* 練習 LeetCode [51. N-Queens](https://leetcode.com/problems/n-queens/)，應選擇 C 語言並善用上述的 `__builtin_ffs`，大幅降低記憶體開銷;

參考 [使用位元計算加速求解 N-Queen](https://medium.com/fcamels-notes/%E4%BD%BF%E7%94%A8%E4%BD%8D%E5%85%83%E8%A8%88%E7%AE%97%E5%8A%A0%E9%80%9F%E6%B1%82%E8%A7%A3-n-queen-d20442f34110) 與 [RinHizakura](https://hackmd.io/@RinHizakura/BkoGM5V8v#%E6%B8%AC%E9%A9%97%E9%A1%8C-3) 同學的程式碼：

    int len;
    char cpy[] = ".................";
    int sizes[] = { 1,      1,       0,        0,        2,        10,    4,
                   40,     92,      352,      724,      2680,     14200, 73712,
                   365596, 2279184, 14772512, 95815104, 666090624 };

    bool recursive(char*** ans, int n, int row, uint16_t l, uint16_t m, uint16_t r) {
        if (row == n) {
            ++len;
            if (len == sizes[n])
                return true;
            for (int i = 0; i < n; ++i) {
                memcpy(ans[len][i], ans[len - 1][i], n);
            }
            return false;
        }
        uint16_t full = 0xffffffff;
        uint16_t accept_cols = ~(l | m | r);
        for (int col = __builtin_ffs(accept_cols) - 1; col < n; col = __builtin_ffs(accept_cols) - 1) {
            ans[len][row][col] = 'Q';
            if (recursive(ans, n, row + 1, (l | (1 << col)) << 1, (m | (1 << col)), (r | (1 << col)) >> 1))
                return true;
            ans[len][row][col] = '.';
            accept_cols ^= (1 << col);
        }
        return false;
    }

    char*** solveNQueens(int n, int* returnSize, int** returnColumnSizes) {
        len = 0;
        char*** ans;
        ans = malloc(sizeof(char**) * sizes[n]);
        *returnColumnSizes = malloc(sizes[n] * sizeof(int));
        *returnSize = sizes[n];
        for (int i = 0; i < sizes[n]; ++i) {
            (*returnColumnSizes)[i] = n;
            ans[i] = malloc(sizeof(char*) * n);
            for (int j = 0; j < n; ++j) {
                ans[i][j] = malloc(sizeof(char) * (n + 1));
                memcpy(ans[i][j], cpy, n);
                ans[i][j][n] = '\0';
            }
        }
        if (!sizes[n])
            return ans;
        recursive(ans, n, 0, 0, 0, 0);
        return ans;        
    }

將可放置皇后的位置 bit 設為 1，使用三個變數 `l`, `m`, `r` 記錄目前不能放置的位置，因為是從 `row 0` 迴圈到 `row n`，所以不用管 `row` 上重複的情形，而 `col` 則是將前一個 `row` 皇后的位置的左右一個 bit，加上中間的 bit 加到 bitmask 中防止選到不能選的位置。

因 bit 為 1 表示可選擇之位置，0 表示不能選擇之位置，故使用 `__builtin_ffs` 能快速找到可選擇的 `col` 位置。

Leetcode 結果：
Runtime: **4 ms**, faster than **96.97%** of C online submissions for N-Queens.
Memory Usage: **7.2 MB**, less than **27.27%** of C online submissions for N-Queens.

###### tags: `linux2020`