# 2020q3 Homework (quiz4)
contributed by < `hankluo6` >
### 測驗 1
```c
int hammingDistance(int x, int y) {
return __builtin_popcount(x OP y);
}
```
#### `OP`
- [x] `(c)` `^`
使用 XOR 會將兩個有差異的 bit 設為 1,相同的設為 0,所以回傳的值就是 hamming distance。
#### 延伸問題
* 舉出不使用 gcc extension 的 C99 實作程式碼
```c
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 相同的情形,進而減少程式碼。
```c
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 外圍的部分。
![](https://i.imgur.com/VuedbNJ.png)
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,這樣可以有效的減少實作的程式碼。
```c
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
```c
#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% 的提交者;
將上述內容實作後並提交:
```c
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
```c
#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
```c
#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;
}
```
![](https://i.imgur.com/NRfjKP9.png)
因 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
```c
#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 皆為相同,只要有涵蓋到原本的範圍就好。
```c
blockidx = offset >> ffs32(divident);
divident >>= ffs32(divident);
```
程式碼中這兩步是把上面公式中的 $2^n$ 部份移除,剩下部份就要依靠 `CASE` 來處理:
```c
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) 同學的程式碼:
```c
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`