# 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
```
Sign in with Wallet
Connect another wallet