2022q1 Homework3 (quiz3)

# 2022q1 Homework3 (quiz3) contributed by < [`tinyynoob`](https://github.com/tinyynoob) > > [作業要求](https://hackmd.io/@sysprog/BJJMuNRlq) ## 測驗 3 ### 解釋程式碼 ```c= #include <stdint.h> uint8_t rev8(uint8_t x) { x = (x >> 4) | (x << 4); x = ((x & 0xCC) >> 2) | (EXP2); x = ((x & 0xAA) >> 1) | (EXP3); return x; } ``` 為了對 8 bit 數字進行翻轉，這裡的作法以類似 top-down 的方式往下處理，每次交換左右半部，直到每個半部的 size 都為 1 即完成所有翻轉，圖解如下： ```graphviz graph { node[shape=record]; parti8[label="7|6|5|4|3|2|1|0"]; } ``` ```graphviz graph { node[shape=record]; parti4a[label="3|2|1|0"]; parti4b[label="7|6|5|4"]; } ``` ```graphviz graph { node[shape=record]; parti2a[label="1|0"]; parti2b[label="3|2"]; parti2c[label="5|4"]; parti2d[label="7|6"]; } ``` ```graphviz graph { node[shape=record]; 0; 1; 2; 3; 4; 5; 6; 7; } ``` 因為 `0xcc` = `11001100` ，所以第 5 行需要補上的就是另一半的操作，因此 `EXP2` 為 `(x & 0x33) << 2` 。同理 `EXP1` 為 `(x & 0x55) << 1` 。 ## 測驗 5 ### 解釋程式碼 ```c= #include <limits.h> int divide(int dividend, int divisor) { int signal = 1; unsigned int dvd = dividend; if (dividend < 0) { signal *= -1; dvd = ~dvd + 1; } unsigned int dvs = divisor; if (divisor < 0) { signal *= -1; dvs = ~dvs + 1; } int shift = 0; while (dvd > (EXP6)) shift++; unsigned int res = 0; while (dvd >= dvs) { while (dvd < (dvs << shift)) shift--; res |= (unsigned int) 1 << shift; EXP7; } if (signal == 1 && res >= INT_MAX) return INT_MAX; return res * signal; } ``` 可以看出第 4 至 15 行在把 `dvd` 跟 `dvs` 都轉成非負的值，並把 sign 值存放到 `signal` 。另外也可以看出缺乏對 divisor 為 0 的檢查。這裡由於是 int 轉成 unsigned int ，所以取二補數並不會因 INT_MIN 導致 undefined behaviour ，以 4 bit 舉例說明： ``` 1000 # INT_MIN of 4-bit signed integer, which is -8 1000 # two's complement of 1000, it is unsigned 8 ``` 因此第 8 及 14 行會得出正確的結果。接著後面就是要想辦法用 shift 和 subtract 來做出除法。因為除法是從最高位開始，所以這裡的第 17~19 行即是要推出最高位並以 `shift` 紀錄，大概會想到兩個可能的答案 1. `dvs << shift` 2. `dvs << shift + 1` 有例子會比較容易推測，所以我們先假定 `dvd` = 16 和 `dvs` = 3 ，則相除結果應為 0101 。除了第 26 行未知以外，唯一在更新答案的就是第 25 行，其功能為將特定位元設為 1 ，因此可以推斷第一次執行到第 25 行時 `shift` 應該要為 2 。而因為第 23 ~ 24 行結束後會保證 `dvd >= (dvs << shift)` ，所以第 17~19 行把 `shift` 起始值稍微設超過也沒關係，因此 `EXP6` 可以選成較為簡潔的 `dvs << shift` 。 ``` 1 -------- 11 / 10000 1100 # (11) << 2 ------ ``` 在繼續除之前，我們需要更新 `dvd` ，所以 `EXP7` 就是 `dvd -= (dvs << shift)` ，得： ``` 1 -------- 11 / 10000 1100 ------ 100 ``` 接著繼續想，會發現 23~24 行幫助我們找到哪一位可以繼續除，全部做完就會像這樣： ``` 101 -------- 11 / 10000 1100 ------ 100 11 ---- 1 ``` 不過因為我們是用 unsigned int 在除，而商的回傳型態是 int ，所以前面所提及的 INT_MIN 還是可能造成問題。一樣以 4-bit 為例， 1000 除以 1 得 1000 ，但因為 1000 原本代表的是負數所以 `signal` 會被標記為 -1 ，那麼到第 31 行就會發生 1000 * (-1) 的狀況 :::warning `abs(INT_MIN)` 是未定義行為，那 `(INT_MIN) * (-1)` 呢？測試似乎會得到原值 ::: :::spoiler 補齊後程式碼 ```c= int divide(int dividend, int divisor) { int signal = 1; unsigned int dvd = dividend; if (dividend < 0) { signal *= -1; dvd = ~dvd + 1; } unsigned int dvs = divisor; if (divisor < 0) { signal *= -1; dvs = ~dvs + 1; } int shift = 0; while (dvd > (dvs << shift)) shift++; unsigned int res = 0; while (dvd >= dvs) { while (dvd < (dvs << shift)) shift--; res |= (unsigned int) 1 << shift; dvd -= (dvs << shift); } if (signal == 1 && res >= INT_MAX) return INT_MAX; return res * signal; } ``` ::: :::warning 不理解 28, 29 行的功能 ::: ### 改進並實作 #### 錯誤處理我發現它沒有對 `divisor` 為 0 的檢查，結果會使程式陷入無窮迴圈。使用一些 corner cases 測試可以發現： ```shell 2147483647 / 1 = 2147483647 2147483647 / -1 = -2147483647 -2147483648 / 1 = -2147483648 -2147483648 / -1 = 2147483647 ``` 在 `INT_MIN / -1` 的情境中，也因為結果出現 32 位有號數無法表示的 $2^{32}$ 導致答案不如預期。為了了解在慣例上這兩種情況該如何處置，我去測試了 standard library 中的 `div()` 。呼叫 `div(7, 0)` 以及 `div(INT_MIN, -1)` 皆得到相同結果： ```shell 浮點數例外 (核心已傾印) ``` 看來它們都是直接用 exception 的方式處理，查詢 `man signal` 找到這段敘述： > Integer division by zero has undefined result. On some architectures it will generate a SIGFPE signal. (Also dividing the most negative integer by -1 may generate SIGFPE.) Ignoring this signal might lead to an endless loop. 我想仿造使用相同的方式，因此在 `divide()` 函式的最前方加入： ```c if (divisor == 0) raise(SIGFPE); else if (dividend == INT_MIN && divisor == -1) raise(SIGFPE); ``` 測試後 signal 可被正確地觸發。 #### 效能改善我認為這裡 while 迴圈判斷 `shift` 的方式，可以使用 `clz()` 來取而代之，程式碼更動如下： ```diff + int shift = __builtin_clz(dvs) - __builtin_clz(dvd); - int shift = 0; - while (dvd > (dvs << shift)) - shift++; unsigned int res = 0; while (dvd >= dvs) { + if (dvd < (dvs << shift)) + shift--; - while (dvd < (dvs << shift)) - shift--; res |= (unsigned int) 1 << shift; dvd -= (dvs << shift); } ``` 嘗試測量 10000 組隨機數進行除法所耗費的時間，得： ```shell number of cases: 10000 origial consumption: 966675 v1 consumption: 148825867 ``` 結果顯示改進版程式碼的耗費時間竟為原版的 ${10}^2$ 倍，與我的預期截然不同。 :::info working on... ::: ## 測驗 7 ### 解釋程式碼給定 32 位的 unsigned `v` ，這題想計算出需要幾個 bit 才能表示 `v` ，顯然 `ret`$\in [0, 32]$ ，且 `ret` 是由 `v` 的 most significant set bit 所決定。根據二進制的原理，我們知道任意 `v` 必定能表示為 $$ \sum_{i=0}^{31}b_i\cdot 2^i \quad \text{where } b_i\in\{0,1\} $$ $i$ 值最大且為 1 的 $b_i$ 就是我們要的 `ret` ，這裡想透過類似二分搜尋的方式來進行查找。首先把 32 個 bit 分為兩半，數字代表 most significant set bit 在由下往上的第幾位。 ```graphviz digraph { whole [label="32 ~ 1"]; left [label="32 ~ 17"]; right [label="16 ~ 1"]; whole -> right; whole -> left; } ``` `v` 至少要是 $$ \text{binary } 1 \underbrace{00\cdots00}_{16個0} $$ 才能落在左邊，因此當 `v` $\geq2^{16}$ 時 `v` 屬於左半，否則 `v` 屬於右半。在屬於左半的情況若是我們先把數字 16 存到某個變數，那麼左右兩個 branch 又變成一樣的問題，那就是求「最大的 1 在 16 位中的何處？」，於是繼續： ```graphviz digraph { whole [label="16 ~ 1"]; left [label="16 ~ 9"]; right [label="8 ~ 1"]; whole -> right; whole -> left; } ``` 完成後一樣把結果加到 `m` 並繼續 4 位的搜尋。 ```c= int ilog32(uint32_t v) { int ret = v > 0; int m = (v > 0xFFFFU) << 4; v >>= m; ret |= m; m = (v > 0xFFU) << 3; v >>= m; ret |= m; m = (v > 0xFU) << 2; v >>= m; ret |= m; m = EXP10; v >>= m; ret |= m; EXP11; return ret; } ``` 在表示上，我們將前面的 `v` $\geq2^{16}$ 更改為等價的 `v` $>2^{16}-1 = \mathtt{0xFFFFU}$ 以利程式看起來的對稱性和二分關係。在每一輪中我們藉由比較和左移把結果暫存到 `m` ，並用 `|` 來立起 `ret` 中對應的 bit ，接著右移 `v` 來讓兩種結果回歸同一個新問題。繼續二分的結果， `EXP10` 顯然為 `(v > 0x3U) << 1` 。最後第 16 行的 `EXP11` 就比較特別了，在程式走到這裡的時候僅剩下 2 個 bit ，所以要判斷是否 $>1$ ，我們可以列出 `v` 為 `00`, `01`, `10`, `11` 的情況來思考，結果應該要分別對應到 0, 1, 2, 2。然而我們目前不能產生兩位數 (已於 13 至 15 行佔用第 2 位)，此外仔細觀察可以發現第 1 位也早在第 3 行就已經被佔用，因此 `0` 和 `1` 的 case 其實已經判斷完畢。那麼 `EXP11` 就應該以加法更新，並只在最後是 `10` 或 `11` 時更新，即為 `ret += (v > 1)` 。 ### 探索 kernel 中的應用 ### 研讀《[Using de Bruijn Sequences to Index a 1 in a Computer Word](http://supertech.csail.mit.edu/papers/debruijn.pdf)》雖然目前 architecture 上大都有提供 fls 或相似的 instruction，但為了 cross-platform 和 cross-compiler 的情境，我們有時仍需要 software 的 fls 處理。上方的 `ilog32()` 程式碼仍可能有 branch 產生，因此我們想找到更好的 branch-free 版本。考慮 $n$-bit 無號整數，此篇文章以 $n=8$、求 ffs 為例。首先，我們需要一組長度為 $n$ 的 [de Bruijn sequence](https://en.wikipedia.org/wiki/De_Bruijn_sequence)，這類特殊 binary sequence 所有 $\log n$ 長度的 cyclic subsequence 完全不重複，也就是每個 $\log n$ 長度的整數都會剛好出現一次。例如這裡給定為 `debru` = `00011101`，它所有長 $\log n$ 的 cyclic subsequence 為：`000`, `001`, `011`, `111`, `110`, `101`, `010`, `100` (由左而右)，我們可以建立 look-up table 以快速查詢位置： |subsequence|index| |:-:|:-:| |`000`|0| |`001`|1| |`010`|6| |`011`|2| |`100`|7| |`101`|5| |`110`|4| |`111`|3| 接著要如何實現 ffs 呢？我們先針對輸入值 weight 為 1 的 case，做法如下： 1. 輸入值 `x` 乘以 `debru` 2. `>> 5` 3. `& 7` 4. 查表 #### 原理以 `x` = `00010000` 為例，首先因為 `x` 的 weight 為 1，將其乘以 `debru` 等同於對 `debru` 左移 ffs(x) ``` 00010000 * 00011101 ----------------- 0000000111010000 @@@ ``` 步驟 2 和 3 固定取出 `@@@` 位置的 bits 可以發現這些操作就等同取出 `debru` 的第 ffs(x) 個 cyclic subsequence，於是反查就可以輕易找到 ffs(x)。 --- 以上僅適用於 weight 為 1 的 `x`，那要如何擴充到任意 `x` 呢？只要先對 `x` 使用 isolate first set 操作就行了，一個經典的方法是 `x = x & (-x)`。 :::success 需要注意到我們這裡的 ffs 並未考量 `x` = 0 ::: 解釋至此， fls() 也可以用相似的技巧配合 isolate last set 進行快速求解。我們取 [stackoverflow 討論串](https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers) 裡的 `debru`，並嘗試自己撰寫出 fls() 。以 32-bit 為目標，首先撰寫輔助程式建出 look-up table ```c static const int debruijn_fls32_table[32] = { 32, 31, 22, 30, 21, 18, 10, 29, 2, 20, 17, 15, 13, 9, 6, 28, 1, 23, 19, 11, 3, 16, 14, 7, 24, 12, 4, 8, 25, 5, 26, 27}; ``` 於是可以寫出： ```c /* most significant = 1 * least significant = 32 * 0 = 0 */ int fls32(uint32_t x) { return (0 - !!x) & debruijn_fls32_table[(ils(x) * 0x07C4ACDDu) >> 27 & 31]; } ``` :::spoiler `ils()` function ```c // isolate last set uint32_t ils(uint32_t n) { n |= n >> 1; n |= n >> 2; n |= n >> 4; n |= n >> 8; n |= n >> 16; return n ^ (n >> 1); } ``` ::: \ 依據之前 branch 版本的慣例，我將 MSB 視為 1、LSB 為 32，並補上 fls(0) = 0。簡單用亂數進行測試，將 `fls32()` 的結果與 `__builtin_clz() + 1` 比對，完全沒出現錯誤回報，唯兩者對 0 會產生不同結果。 :::spoiler test code ```c #define TESTNUM 10000000 int main() { srand(time(NULL)); puts("start test:"); for (int i = 0; i < TESTNUM; i++) { uint32_t x = rand() << 16 | rand(); if (__builtin_clz(x) + 1 != fls32(x)) printf("different answer when x = %d\n", x); } for (int i = 0; i < TESTNUM; i++) { uint32_t x = rand(); if (__builtin_clz(x) + 1 != fls32(x)) printf("different answer when x = %d\n", x); } puts("end test."); return 0; } ``` ::: ### 編譯時期的 `ilog2()` 實作 ## 測驗 10 ### 解釋程式碼 ```c DIVIDE_ROUND_CLOSEST(x, divisor) ``` 此題想在整數除法中選擇較接近的整數作為答案，而非原本的無條件捨去。先考慮正數的例子，從數學上來思考，假設除法結果落在 `@` 會得到 1，而落在 `#` 會得 2，從數線看起來如： ``` 0 1 2 3 |------|------|------| @@@@@@@####### traditional computer result @@@@@@@####### our goal ``` 假設除出來位在 ``` 0 1 2 3 |------|------|------| x ``` 如果我們能夠把它右移 0.5 ，那麼 `x` 就會落入傳統 `@` 的範圍並被電腦映射到 1。然而在整數中根本沒有 0.5 的概念，而追究 0.5 的來源即是整數除法中無法整除的部份，因此我們應想辦法去調整被除數。假設正整數 $a, b$ 在電腦上做除法得 $a/b=c$ ，我們希望在算出 $c$ 之後可以將它右移 $0.5$ ，那麼可以將原算式改造成 $(a+b*0.5)/b$ 達成目的， $*0.5$ 在電腦上可以使用 `>> 1` 。 --- 接著我們需要將負整數除法也納入討論，在本機使用 gcc 編譯測試顯示 `/` 運算的結果會對稱於原點 0 ： ``` -6 / 3 = -2 -5 / 3 = -1 -4 / 3 = -1 -3 / 3 = -1 -2 / 3 = 0 -1 / 3 = 0 0 / 3 = 0 1 / 3 = 0 2 / 3 = 0 3 / 3 = 1 4 / 3 = 1 5 / 3 = 1 6 / 3 = 2 ``` 注意負數部份，例如 `-4/3` 得 -1 而非 -2 。 :::warning 嘗試搜尋 C99 規格書，看到除法運算寫在 6.5.5 章節，但沒有查到對於商或餘數選擇的敘述 ::: 由於這種對稱性，因此前述對於正數右移 0.5 的工作，到了負數情況便改成要左移 0.5 。回頭來看題目： ```c= #define DIVIDE_ROUND_CLOSEST(x, divisor) \ ({ \ typeof(x) __x = x; \ typeof(divisor) __d = divisor; \ (((typeof(x)) -1) > 0 || ((typeof(divisor)) -1) > 0 || \ (((__x) > 0) == ((__d) > 0))) \ ? ((RRR) / (__d)) \ : ((SSS) / (__d)); \ }) ``` 因為 `-1 = UINT_MAX > 0` ，所以表達式 `((typeof(x)) -1) > 0` 可用於判斷該變數是否為 unsigned ；而 `(((__x) > 0) == ((__d) > 0)))` 檢驗兩者是否同號。因為 C 語言會將 unsigned 和 signed 混合運算的結果轉為 unsigned ，因此只要被除數或除數中有 unsigned number ，那結果必定非負。若兩者同號那相除結果也會非負。於是乎我們可以推測 `RRR` 處理 $\geq 0$ 的結果， `SSS` 處理 $<0$ 的結果。根據之前的討論，填入： - `RRR` : `__x + (__d >> 1)` - `SSS` : `__x - (__d >> 1)` ### 探索 kernel 中的應用 :::info 待補 ::: ## 測驗 11 ### 解釋程式碼這裡我跳過 `fls()` 的部份，直接針對開平方根的部份進行解釋，其功能等同 $\lfloor\sqrt{x}\rfloor$ 。 ```c= unsigned long i_sqrt(unsigned long x) { unsigned long b, m, y = 0; if (x <= 1) return x; m = 1UL << (fls(x) & ~1UL); while (m) { b = y + m; XXX; if (x >= b) { YYY; y += m; } ZZZ; } return y; } ``` 對於 $x$ 滿足 $\mathrm{fls}(x) = n$ ，我可以猜到 $\sqrt{x}\geq \lfloor n/2\rfloor$ ，但是我不知道要怎麼算出更細節的結果。我在 [wiki](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 找到相關的說明，它的概念是利用 $(a+b)^2 = a^2+2ab+b^2$ 公式，想辦法在已知 $a$ 的情況下去湊出 $b$ 。這題的變數和要填空的格子有點多，很難依可見的資訊判斷作答，我研究了好一陣子都很茫然，所以我直接去看答案。 ```c= unsigned long i_sqrt(unsigned long x) { unsigned long b, m, y = 0; if (x <= 1) return x; m = 1UL << (fls(x) & ~1UL); while (m) { b = y + m; y >>= 1; if (x >= b) { x -= b; y += m; } m >>= 2; } return y; } ``` 以 $x = (120)_{10} = (1111000)_2$ 為例來說明：首先`fls(x)` 會算出 6 ，接著 `& ~1` 的作用是把最低的位元改成 0 ，也就是取出小於等於 `fls(x)` 的最大偶數，以我們的例子會輸出 6 。最後把 `m` 初始化為 $2^6$ ，稍候就會發現 `m` 變數是用來指示目前運算到第幾位。因為 return `y` 所以 `y` 應是存開根號結果的地方。 ``` y = 1 ----------------- \/ 1 11 10 00 = x b = 1 ----------- 0 11 ``` 每次它會算出 `b` 之後檢查 `x >= b` ，如果成立那 `x` 便可以減 `b` 並得該位的 `y` 是 1 ，否則該位 `y` 是 0 。繼續 loop 之前執行第 16 行，目的是把 `m` 對齊下一個要運算的位。 ``` y = 1 ----------------- \/ 1 11 10 00 1 ----------- 0 11 = x b = 1 01 ```

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.