# 2020q3 Homework3 (quiz3)
contributed by < `Holycung` >
[2020q3 Homework3 (quiz3) 題目](https://hackmd.io/@sysprog/2020-quiz3)
## 目錄
[TOC]
## 測驗 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，包含 Cppcheck 在內的靜態分析器會嘗試找出 C 程式中這項 unspecified behavior。考慮到程式碼的相容性，我們想實作一套可跨越平台的 ASR，假設在二補數系統，無論標的環境的 shift 究竟輸出為何，我們的 ASR 總可做到:
$−5 ≫^{arith} 1 = {−3}$
上述 $−3$ 正是 $\frac{−5}{2}$ 的輸出，並向 $−∞$ 進位 (rounding)。
針對有號整數的跨平台 ASR 實作如下:
```cpp=
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));
}
```
請補完程式碼。
### 作答
回答這題前要先了解到算術右移 (arithmetic right shift) 有所了解，$−5 ≫^{arith} 1 = {−3}$，在 binary 當中的表示如下，以下皆以 8 個位元方便表示。
| decimal | binary |
| -------- | -------- |
| -3 | 11111101 |
| -5 | 11111011 |
觀察程式碼一開始要先檢查平台的實作 shift 方式才能達到支援跨平台，第 3 行就透過將 $-1$ 往右位移，判斷是否大於零，如果是 ars 了話就會是 $-1$，如果是邏輯位移的話就會變成大於零的數如下，藉此就可以用 `logical` 判斷出平台是否支援算術位移。
|| decimal | binary |
|--------| -------- | -------- |
|原本的值| -1 | 11111111 |
|算數右移一位| -1 | 11111111 |
|邏輯右移一位| 127 | 01111111 |
第 4 行看到 `fixu`，這邊是判斷是否要做 fix signed bit 的動作，除了剛剛的 `logical` 判斷，還要是今天 input 的 $m < 0$ 才需要 fix，所以這邊 OP1 = m < 0。如果兩個都成立的話，就會得到 $-1$，不成立的話是 $0$，但是這邊變數 `fixu` 是 `unsigned int`，所以要想一下這邊的型態轉換，可以看到 [c99 規格書](http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1124.pdf)當中說明。
**6.3.1.3 Signed and unsigned integers**
> 1 When a value with integer type is converted to another integer type other than _Bool, if
the value can be represented by the new type, it is unchanged.
>
> 2 Otherwise, if the new type is unsigned, the value is converted by repeatedly adding or subtracting one more than the maximum value that can be represented in the new type
until the value is in the range of the new type.
在 C 語言當中做型別轉換，假設要把一個 M 位的類型（值為X）轉換成一個 N 位的類型，會先判斷兩個變數的大小做 interger promotion，然後在根據 M N 大小，會有不同的做法，細節可以參考 [C 類型轉換](http://shihyu.github.io/books/ch15s03.html)。那我們今天只要探討 signed integer to unsigned integer ，對應的做法就會是 $X\ \%\ 2^N$，就如同規格書第 2 點說的一樣。
所以如果是得到 $0$ 了話沒差，得到 $-1$ 的話 `fixu` 就會轉換成每個 bit 都是 1，因為原本的 `-1` 在 binary 當中是 `1111 1111`。
第五行 `int fix = *(int *) &fixu;` 把 fixu 轉回 signed int 這邊不太懂，看了老師在 [RinHizakura](https://hackmd.io/@RinHizakura/SkjPS8vBP) 的註解是**可避免編譯器進行過度 (aggressive) 最佳化**。
最後 `return (m >> n) | (fix ^ (fix >> n));` ，對 m shift 完後，or 上 fix 的值就會得到 ars 的結果了。
## 測驗 2
### 題目
在 C 語言：數值系統篇 中，我們介紹過 power of 2 (漢字可寫作「二的冪」)。
我們可用 __builtin_ctz 來實作 LeetCode 342. Power of Four，假設 int 為 32 位元寬度。實作程式碼如下:
```cpp=
bool isPowerOfFour(int num)
{
return num > 0 && (num & (num - 1))==0 &&
!(__builtin_ctz(num) OPQ);
}
```
### 作答
分成三個部分討論 `num > 0`、`(num & (num - 1))==0`、`!(__builtin_ctz(num) OPQ)`
- `num > 0` 要確保 num 是大於 0， 0 不是 power of four
- `(num & (num - 1))==0` 這個是判斷 num 是否為 power of two，因為如果是 power of four 就必定是 power of two，而這個是在 [你所不知道的 C 語言：數值系統篇](https://hackmd.io/@sysprog/c-numerics#%E9%81%8B%E7%94%A8-bit-wise-operator) 有講到的，如果 num 是 power of two，其 binary 的最高位是 1 其餘為 0，利用這個特性把 `num - 1`，就是借位的概念把後面的 bit 全部反轉，因此 `num & (num - 1)` 的值必定為 0。
- `!(__builtin_ctz(num) OPQ)` 從 [GCC extension](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html) 當中可以看到 `builtin_ctz` 會回傳 trailing zero 的數量，那我們可以透過觀察下表發現 power of four 的 ctz 數量必定是偶數，偶數的最小位元必定是 0，所以透過 bitwise and `0x1` 會得到 0 奇數則會得到 1，所以最後再加上 not 就完成判斷是否為 power of four。
> **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.
| decimal | binary |
| -------- | -------- |
| $4^1$ | 00000100 |
| $4^2$ | 00010000 |
| $4^3$ | 01000000 |
### 延伸問題
改寫上述程式碼，提供等價功能的不同實作，儘量降低 branch 的數量
#### 作答
參考 [RinHizakura](https://hackmd.io/@RinHizakura/SkjPS8vBP) 使用 gcc extension 的 `__builtin_popcount` 和 `__builtin_ffs`，前者 popcount 是計算有多少個位元為 1，如果是 power of four 了話就要只會是 1。後者 `ffs` 的意思是 `find first set`，會回傳從右起第一個 1 的位置，如果是 power of four 了話就必定是奇數。
```cpp=
bool isPowerOfFour_new(int num)
{
return !(__builtin_popcount(num) ^ 1) && __builtin_ffs(num) & 1;
}
```
:::info
TODO
- 練習 LeetCode 1009. Complement of Base 10 Integer 和 41. First Missing Positive，應善用 clz;
- 研讀 2017 年修課學生報告，理解 clz 的實作方式，並舉出 Exponential Golomb coding 的案例說明，需要有對應的 C 原始程式碼和測試報告;
:::
## 測驗 3
### 題目
針對 [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 位元，程式碼列表如下
```cpp=
int numberOfSteps (int num)
{
return num ? AAA + 31 - BBB : 0;
}
```
### 作答
這一題從位元的角度來看，只要看最低位元是 0 或 1 來判斷。
- 最低位是 0，代表直接向右 shift，就等於除二，需要一個 step。
- 最低位是 1，要先做減一再向右 shift，需要兩個 steps。除非是最後只剩一個時候，只要做減一一個 step 就好。
所以我們最多只會 shift `32 - __builtin_clz(num)` 次，然後如果有遇到 1 了話要多一個 step，所以要加上 `__builtin_popcount(num)` 的次數，但是最後一定會遇到 1，所以要再減掉一個 step，最後就是 `32 - __builtin_clz(num) + __builtin_popcount(num) - 1`，化減後就是 `__builtin_popcount(num) + 31 - __builtin_clz(num)`。
### 延伸問題
:::info
TODO
- 看完 popcnt_branchless 的實作
- 改寫上述程式碼，提供等價功能的不同實作並解說;
- 提示: Bit Twiddling Hacks 中提及許多 bitwise operation 的使用，如 bit inverse、abs 等等，是極好的參考材料。
- 避免用到編譯器的擴充功能，只用 C99 特徵及 bitwise 運算改寫 LeetCode 1342. Number of Steps to Reduce a Number to Zero，實作出 branchless 解法，儘量縮減程式碼行數和分支數量;
:::
## 測驗 4
### 題目
考慮以下 64-bit GCD (greatest common divisor, 最大公因數) 求值函式:
```cpp=
#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;
}
```
改寫為以下等價實作:
```cpp=
#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;
}
```
### 作答
作答之前要先了解 gcd 的一些特性，參考 [Greatest Common Divisor 特性和實作考量](https://hackmd.io/W8BtJFQdSZekR_O53ICTyQ?view) 還有 [Euclidean_algorithm](https://en.wikipedia.org/wiki/Euclidean_algorithm)。
第 5 行 for 迴圈做的事情是對 $u$、$v$ 同除 $2$，`!((u | v) & 1)` 的作用在於檢查是不是兩個都是偶數，只要有一個不是偶數就停止，只要兩個數都是偶數就代表可以同時被 $2$ 整除，並且用 `shift` 記下有幾次。
第 8 行跟第 11 行的 while 迴圈則都是再判斷 $u$、$v$ 是否為偶數，如果是了話就除二直到變成奇數，因為 $u$、$v$ 已經不會有 2 的公因數了，所以 $\text{gcd}(x, y)=\text{gcd}(x, \frac{y}{2})$。
第 10 行開始的 while 迴圈其實就是再做輾轉相除法(Euclidean algorithm)，這個特性如下，假設 $b_0 > a_0$
\begin{aligned}
\begin{split}
\gcd(a_0, b_0) &= \gcd(b_0 \% a_0 = \underline{b_1}, a_0) \\
\gcd(\underline{b_1}, a_0) &= \gcd(a_0 \% b_1 = \underline{a_1}, b_1) \\
\gcd(\underline{a_1}, b_1) &= \gcd(b_1 \% a_1 = b_2, a_1) \\
&\vdots \\
\gcd(a_i, b_i) &= \gcd(b_i \% a_i = a_{i+1}, b_i) \\
&\vdots \\
\gcd(\mathbf{0}, b_n) &= |b_n|
\end{split}
\end{aligned}
所以這邊可以寫成，當 $u > v$， $\text{gcd}(u, v) = \text{gcd}(u-v, v)$，每次都把大的數減掉小的，不過要注意的是 $u、v$ 都是奇數相減一定會得到偶數，所以要確保下次進入迴圈的時候 $u$ 必定是奇數，所以每次迴圈只要把 $v$ 除以 2 再繼續，直到 u 跟 v 兩數相同，最後相減完把 0 放在 v，所以可以知道終止條件就是 v = 0，也就是 XXX = `v`。
最後 return 的時候就要把輾轉相除法的結果 u 乘上最開始除 2 的次數，也就是 YYY = `u << shift`。
### 延伸問題
在 x86_64 上透過 `__builtin_ctz` 改寫 GCD，分析對效能的提升;
#### 作答
把前面除二的動作換成 `__builtin_ctz` 的作法，並且把除二都換成了 shift 的動作加速。
```cpp=
#define min(a,b) \
({ __typeof__ (a) _a = (a); \
__typeof__ (b) _b = (b); \
_a > _b ? _b : _a; })
uint64_t faster_gcd64(uint64_t u, uint64_t v) {
if (!u || !v) return u | v;
int shift = min(__builtin_ctz(u), __builtin_ctz(v));
u >>= shift;
v >>= shift;
while (!(u & 1))
u >>= __builtin_ctz(u);
do {
while (!(v & 1))
v >>= __builtin_ctz(v);
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u << shift;
}
```
對原本的 `gcd` 跟這個加速過的版本 `faster_gcd64` 進行比較，範圍從 1~50000，對所有的組合找最大公因數並計算其所花費的時間(單位 ns)，詳細的程式可以參考 [plot.py](https://github.com/Holychung/2020q3_System_Software/blob/main/quiz3/plot.py)，結果如下。
可以看出 `faster_gcd64` 在數字越來越大的時候成長的斜率明顯低於 `gcd`，明顯有比較快。
## 測驗 5
### 題目
在影像處理中，[bit array](https://en.wikipedia.org/wiki/Bit_array) (也稱 bitset) 廣泛使用，考慮以下程式碼:
```cpp=
#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 改寫為效率更高且等價的程式碼:
```cpp=
#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;
}
```
請補完程式碼。
### 作答
先觀察 `naive` 這段程式碼在做的事情是把 `bitmap` 當中是 1 的 bit 位置記錄到 `out` 然後 return 總共有幾個，也就是 `pos` 這個變數。
這段程式可以藉由 ctz 的數量來更快速的找到下一個是 1 的 bit 位置，也就是最右邊為 1 的 bit，找到之後記錄下來把那個 bit 設為 0，再繼續找下一個，直到整個 `bitset` 等於 0。
所以 KKK = `bitset & -bitset`，這邊 `-bitset` 就是 `~bitset+1`，然後 `bitset & -bitset` 就會只剩下最右邊為 1 的 bit 是 1 其他 bit 都是 0，最後 `bitset ^= t` 的時候就會把最右邊是 1 的 bit 變為 0，如下方例子。
| variable | binary |
| -------- | -------- |
| bitset | 0010 0100 |
| -bitset | 0001 1100 |
| bitset & -bitset (t) | 0000 0100 |
| bitset ^ t | 0010 0000 |
:::info
TODO
- 設計實驗，檢驗 clz 改寫的程式碼相較原本的實作有多少改進？應考慮到不同的 bitmap density;
- 思考進一步的改進空間;
:::
## 參考資料
- [c99 規格書](http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1124.pdf)
- [C 類型轉換](http://shihyu.github.io/books/ch15s03.html)
- [你所不知道的 C 語言：數值系統篇](https://hackmd.io/@sysprog/c-numerics#%E9%81%8B%E7%94%A8-bit-wise-operator)
- [GCC extension](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html)
- [Greatest Common Divisor 特性和實作考量](https://hackmd.io/W8BtJFQdSZekR_O53ICTyQ?view)
- [Euclidean_algorithm](https://en.wikipedia.org/wiki/Euclidean_algorithm)
- [bit array](https://en.wikipedia.org/wiki/Bit_array)
Sign in with Wallet
Connect another wallet