---
tags: linux2022
---
# 2022q1 Homework2 (quiz2)
contributed by ???
> [第 2 周測驗題](https://hackmd.io/@sysprog/linux2022-quiz2)
### 測驗`1` 兩個無號整數取平均值
#### 說明及補完程式碼
* 考慮以下對二個無號整數取平均值的程式碼:
```c
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
return (a + b) / 2;
}
```
這個直覺的解法會有 overflow 的問題，若我們已知 `a`, `b` 數值的大小，可用下方程式避免 overflow:
```c
#include <stdint.h>
uint32_t average(uint32_t low, uint32_t high)
{
return low + (high - low) / 2;
}
```
接著我們可改寫為以下等價的實作:
```c
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
return (a >> 1) + (b >> 1) + (EXP1);
}
```
其中 `EXP1` 為 `a`, `b`, `1` (數字) 進行某種特別 bitwise 操作，限定用 OR, AND, XOR 這三個運算子。
兩個整數 `a`,`b` 做 `/` 運算，在 C 語言當中為整數除法，e.g. 5 / 2 = 2.
而對兩數分別作 `>>` 運算，形同對兩數作 `/ 2` 運算。因此若是遇到奇數，作除法後小數會被捨去。本來 `（a+b）/ 2` 只捨去一次，現在分別對`a`,`b` 作除法，最多可能被捨去兩次，因此 `EXP1` 是為了補足被捨去超過一次的情況。
在 `a/2 + b/2` 之下，探討 `a`,`b` 的奇偶：
|a 的最右邊bit|b 的最右邊bit|小數被捨去次數 |期望補上的值|
|---|---|---|---|
|0|0|0|0|
|1|0|1|0|
|0|1|1|0|
|1|1|2|1|
期望補上的值，可以看出是 `a 的最右邊bit` `&` `b 的最右邊bit` ，取得 `a 的最右邊bit` 為 `a & 1` ，因此 ==EXP1== = `(a & 1) & (b & 1)` = `a & b & 1`
>參考 [laneser](https://hackmd.io/@laneser/linux2022-quiz2)
* 我們再次改寫為以下等價的實作:
```c
uint32_t average(uint32_t a, uint32_t b)
{
return (EXP2) + ((EXP3) >> 1);
}
```
其中 `EXP2` 和 `EXP3` 為 `a` 和 `b` 進行某種特別 bitwise 操作，限定用 OR, AND, XOR 這三個運算子。
`EXP2` 和 `EXP3` 必定包含 `a` 和 `b` 字元
參考你所不知道的 [C 語言：數值系統篇](https://hackmd.io/@sysprog/c-numerics#%E7%AE%97%E8%A1%93%E5%AE%8C%E5%85%A8%E5%8F%AF%E7%94%A8%E6%95%B8%E4%BD%8D%E9%82%8F%E8%BC%AF%E5%AF%A6%E5%81%9A)其中的「算術完全可用數位邏輯」節，`a + b` 可拆成兩個運算，`a & b` 算的是進位（carry），`a ^ b` 算的是相加但不進位，可用以下真值表說明：
|a|b|a&b|a^b
|-|-|-|-|
|0|1|0|1|
|1|0|0|1|
|0|0|0|0|
|1|1|1|0|
因此 `a + b` = `((a & b) << 1) + (a ^ b)`，而根據題目，要求的是 `(a + b)/2` ，所以將上述式子` / 2`：`(((a & b) << 1) + (a ^ b)) / 2` = `(((a & b) << 1) + (a ^ b)) >> 1` = `(a & b) + ((a ^ b) >> 1)`，所以 ==EXP2== = `a & b` ， ==EXP3== = `a ^ b`
:::warning
TODO:
* 比較上述實作在編譯器最佳化開啟的狀況，對應的組合語言輸出，並嘗試解讀 (可研讀 CS:APP 第 3 章)
* 研讀 Linux 核心原始程式碼 include/linux/average.h，探討其 Exponentially weighted moving average (EWMA) 實作，以及在 Linux 核心的應用
> 移動平均（Moving average)，又稱滾動平均值、滑動平均，在統計學中是種藉由建立整個資料集合中不同子集的一系列平均數，來分析資料點的計算方法。
移動平均通常與時間序列資料一起使用，以消除短期波動，突出長期趨勢或周期。短期和長期之間的閾值取決於應用，移動平均的參數將相應地設置。例如，它通常用於對財務數據進行技術分析，如股票價格、收益率或交易量。它也用於經濟學中研究國內生產總值、就業或其他宏觀經濟時間序列。
:::
---
### 進入測驗 `2` 前，先來看取絕對值與 min / MAX 問題
參考[解讀計算機編碼](https://hackmd.io/@sysprog/binary-representation#%E5%B8%B8%E6%95%B8%E6%99%82%E9%96%93%E5%AF%A6%E4%BD%9C)其中的「常數時間實作」篇，提到**絕對值**的 bit-wise 操作：
$$
abs(x) =
\begin{cases}
x & \text{if } x \geq 0 \newline
(\sim x) + 1 =\ \sim (x-1)& \text{if } x < 0
\end{cases}
$$
對應的實作可以為 (注意此實作 `x` 不能為 `INT32_MIN` )：
```c
#include <stdint.h>
int32_t abs(int32_t x) {
int32_t mask = (x >> 31);
return (x + mask) ^ mask; //if x < 0, return ~(x-1), x 對 -1 作 ^, 形同 x 反轉位元
}
```
另外也提到**給定兩個有號整數，找出其中較小的數值**（注意這實作僅在 `INT32_MIN <= (a - b) <= INT32_MAX` 才成立）：
```c
#include <stdint.h>
int32_t min(int32_t a, int32_t b) {
int32_t diff = a - b;
return b + (diff & (diff >> 31));
}
```
若 `a - b = diff < 0`，則 `a` 較小，`diff >> 31` 為 `0xffffffff`，與 `diff` `&` 運算後仍為 `diff`，回傳 `b + diff = a`
若 `a - b = diff > 0`，則 `b` 較小，`diff >> 31` 後為`0x00000000`，與 `diff` `&` 運算後為 `0x00000000`，回傳 `b + 0x00000000 = b`
---
### 測驗`2` 找 MAX
#### 說明及補完程式碼
* 改寫[〈解讀計算機編碼〉](https://hackmd.io/@sysprog/binary-representation#%E5%B8%B8%E6%95%B8%E6%99%82%E9%96%93%E5%AF%A6%E4%BD%9C)一文的「不需要分支的設計」一節提供的程式碼 min，我們得到以下實作 (max):
```c
#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b)
{
return a ^ ((EXP4) & -(EXP5));
}
```
其中 `EXP4` 為 `a` 和 `b` 進行某種特別 bitwise 操作，限定用 OR, AND, XOR 這三個運算子。
`EXP5` 為 `a` 和 `b` 的比較操作，亦即 logical operator，限定用 >, <, ==, >=, <=, != 這 6 個運算子。
`EXP4`, `EXP5` 必定包含 `a` 和 `b` 字元。
`EXP5` 前有個負號，直覺 `EXP5` 應為 `0` 或 `1`，而若要找較大者，必定是比較 `a`, `b` 大小，假設 ==EXP5== = `a < b`，再探討下述兩種情況：
`a > b`：應該回傳 `a = a ^ 0`（a 對 0 xor，形同無作用），所以 `0 = (EXP4) & -(EXP5) = (EXP4) & 0` --- 一式。
`a < b`：應該回傳 `b = a ^ (a ^ b
)`（b 對 a 作兩次 xor 形同沒作用），所以 `a ^ b = (EXP4) & -(EXP5) = (EXP4) & -1` --- 二式。
解二式 ==EXP4== = `a ^ b`。代入一式 `(a ^ b) & 0` 也確實等於 `0` 。
因此答案為：
```c
return a ^ ((a ^ b) & -(a < b));
```
:::warning
TODO:
* 針對 32 位元無號/有號整數，撰寫同樣 branchless 的實作
不論是有號或無號上述的方法都可以使用，只需要將型別改成有號即可：`uint32_t` -> `int32_t`
:::
---
### 測驗`3` GCD
#### 說明及補完程式碼
* 考慮以下 64 位元 GCD (greatest common divisor, 最大公因數) 求值函式:
```c
#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;//回傳非 0 的參數
/*輾轉相除法概念*/
while (v) {
uint64_t t = v;
v = u % v;
u = t;
}
return u;
}
```
* 註：GCD 演算法
1. If both $x$ and $y$ are $0$, gcd is $0$. $gcd(0,0)=0$, $gcd(0,0)=0$.
2. $gcd(x,0)=x$ and $gcd(0,y)=y$, because everything divides $0$.
3. If $x$ and $y$ are both even, $gcd(x,y)=2∗gcd(\dfrac{x}{2},\dfrac{y}{2})$, because $2$ is a common divisor. Multiplication with $2$ can be done with bitwise shift operator.
4. If $x$ is even and $y$ is odd, $gcd(x,y)=gcd(\dfrac{x}{2},y)$.
5. On similar lines, if $x$ is odd and $y$ is even, then $gcd(x,y)=gcd(x,\dfrac{y}{2})$. It is because $2$ is not a common divisor.
* 改寫為以下等價實作:
```c
#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 (COND);
return RET;
}
```
* 第一個 if 為上述註的第一及第二點。
* 執行 for loop 的條件為 `u`,`v` 兩數皆為偶數（`u`, `v` 任一數二進制的最低位元為 1，條件便無法成立），若皆為偶數兩數都 `/= 2`。此為上述註的第三點。
* while loop 執行的為上述註的第四點。（$x$, $y$ 對應 `u`,`v`）
* do while loop 當中，第一個 while loop 執行的是註的第五點。接下來的 if else 為輾轉相除法的實作，結束條件為當較小的數為 0 的時後，這裡較小的數為 `v`，因此 ==COND== = `v`，且當較小數為 0 時，最大公因數即為另一較大的數。但需要注意的是，在執行上面的 for loop 時，有計算`多少 2 的冪` = `shift` ，使得 $2^{shift}$ 為公因數， 因此 ==RET== = `u << shift`。
#### 在 x86_64 上透過 `__builtin_ctz` 改寫 GCD，分析對效能的提升
參考你所不知道的 [C 語言：數值系統篇](https://hackmd.io/@sysprog/c-numerics#Count-Leading-Zero)其中的「Count Leading Zero」節，`__builtin_ctz` 算的是一數的二進制表示法中，從最低位開始遇到第一個 `1` 之前，共有幾個 `0`。
知道一數的 ctz（count tailing zero）後，便可以知道其中的因數為 `2的多少次冪`。
以相同上述實作的邏輯改寫後：
```c=
#include <stdint.h>
uint64_t mygcd64(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
int ctz_u, ctz_v, shift;
ctz_u = __builtin_ctz(u);
ctz_v = __builtin_ctz(v);
if (ctz_u > ctz_v) {
shift = ctz_u - (ctz_u - ctz_v);
} else {
shift = ctz_v - (ctz_v - ctz_u);
}
u = u >> shift;
v = v >> shift;
if (ctz_u = __builtin_ctz(u)) u = u >> ctz_u;
do {
if (ctz_v = __builtin_ctz(v)) v = v >> ctz_v;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while(v);
return u << shift;
}
```
* 第 10 到第 16 行，做的是上份實作的 for loop 的功能，算出公因數：` 2 的多少次冪`。
* 第 18 及 21 行，執行的是上份實作兩個 while loop 的功能，若有偶數，將 2 的因數除掉。
* 參考[Linux效能分析工具：Perf](https://hackmd.io/@sysprog/gnu-linux-dev/https%3A%2F%2Fhackmd.io%2Fs%2FB11109rdg)，根據以下 `main function` ：
```c
int main()
{
uint64_t u = 64, v = 52;
printf("%ld\n",gcd64(u,v));
//printf("%ld\n",mygcd64(u,v));
return 0;
}
```
改寫**前**實作效能如下：
```shell
$ perf stat --repeat 5 -e cache-misses,cache-references,instructions,cycles ./quiz
> Performance counter stats for './quiz' (5 runs):
22,466 cache-misses # 35.419 % of all cache refs ( +- 4.84% )
63,429 cache-references ( +- 1.88% )
825,880 instructions # 0.56 insn per cycle ( +- 0.40% )
1,474,393 cycles ( +- 3.79% )
0.0007654 +- 0.0000326 seconds time elapsed ( +- 4.26% )
```
改寫**後**實作效能如下：
```shell
$ perf stat --repeat 5 -e cache-misses,cache-references,instructions,cycles ./mine
> Performance counter stats for './mine' (5 runs):
23,806 cache-misses # 36.953 % of all cache refs ( +- 2.95% )
64,423 cache-references ( +- 5.30% )
830,068 instructions # 0.69 insn per cycle ( +- 0.50% )
1,199,807 cycles ( +- 5.34% )
0.001786 +- 0.000135 seconds time elapsed ( +- 7.53% )
```
兩者數據結果差不多，接著把執行次數提升到 100，並多偵測 branch event：
```shell
$ perf stat --repeat 100 -e cache-misses,cache-references,instructions,cycles,branch-misses,branch-instructions ./quiz
> Performance counter stats for './quiz' (100 runs):
22,362 cache-misses # 40.160 % of all cache refs ( +- 1.78% )
55,682 cache-references ( +- 0.76% )
830,198 instructions # 0.62 insn per cycle ( +- 0.10% )
1,339,491 cycles ( +- 1.86% )
6,600 branch-misses # 3.91% of all branches ( +- 0.39% )
168,835 branch-instructions ( +- 0.10% )
0.0007404 +- 0.0000235 seconds time elapsed ( +- 3.17% )
```
```shell
$ perf stat --repeat 100 -e cache-misses,cache-references,instructions,cycles,branch-misses,branch-instructions ./mine
> Performance counter stats for './mine' (100 runs):
22,054 cache-misses # 40.183 % of all cache refs ( +- 1.38% )
54,885 cache-references ( +- 0.61% )
830,690 instructions # 0.65 insn per cycle ( +- 0.12% )
1,270,811 cycles ( +- 1.71% )
6,630 branch-misses # 3.92% of all branches ( +- 0.42% )
168,972 branch-instructions ( +- 0.11% )
0.0006431 +- 0.0000144 seconds time elapsed ( +- 2.23% )
```
兩份實作 cache miss 約 40.1%，branch miss 約 3.9％，數據結果差異不大，改寫後的程式尚須精進。
:::warning
TODO:
* 為何改寫後效能仍尚未改進？
* Linux 核心中也內建 GCD (而且還不只一種實作)，例如 lib/math/gcd.c，請解釋其實作手法和探討在核心內的應用場景。
:::
---
### 測驗`4` 計算 bitmap 中有多少個1, 及分別在第幾位
#### 說明及補完程式碼
在影像處理中，[bit array](https://en.wikipedia.org/wiki/Bit_array) (也稱 bitset) 廣泛使用，考慮以下程式碼：
```c
#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;
}
```
此段程式中，`pos` 代表整個 `bitmap` 中有多少個 `1`。`out` 紀錄的是每個`為 1 的 bit` 在 `bitmap` 的第幾位，這裡可以看出即使 `bitmap` 為 64 * `bitmapsize`，在儲存空間上仍是一維的。
考慮 GNU extension 的 [`__builtin_ctzll`](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html) 的行為是回傳由低位往高位遇上連續多少個 `0` 才碰到 `1`。
用以改寫的程式碼如下：
```c=
#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 = EXP6;
int r = __builtin_ctzll(bitset);
out[pos++] = k * 64 + r;
bitset ^= t;
}
}
return pos;
}
```
第 9 行的作用是僅保留最低位元的 1，並紀錄到 `t` 變數。舉例來說，若目前 bitset 為 $000101_b$，那 `t` 就會變成 $000001_b$，接著就可以透過 `XOR` 把該位的 `1` 清掉，其他保留 (此為 XOR 的特性)。
若 bitmap 越鬆散 (即 1 越少)，於是 **improved** 的效益就更高。
Hint：需考慮到 `-bitwise` 的特性
根據此目的，可以先算出 `bitset` 的 ctz，把 `1L` 左移 ctz 個，如此可以得出只保留最低位 `1` 的 `bitset`。因此 ==EXP6== = `1L << __builtin_ctzll(bitset)`。
但若要利用 `-bitwise` 的特性，需要換個作法：若把 `bitset 取負`，與 `bitset` 差異為除了最低位的 `1` 以及更低位的 `0` 相同之外，其他位皆相反。因此若要得到只保留最低位 `1` 的 `bitset`，把 `bitset & -bitset` 即可。==EXP6== = `bitset & -bitset`。
第 12 行做的便是把 `bitset` 最低位 `1` 消除，讓下一個迴圈再次去找最低位 `1` 在哪。
:::warning
TODO:
* 這樣的程式碼用在哪些真實案例中
* 設計實驗，檢驗 ctz/clz 改寫的程式碼相較原本的實作有多少改進？應考慮到不同的 [bitmap density](http://www.vki.com/2013/Support/docs/vgltools-3.html)
* 思考進一步的改進空間
* 閱讀 [Data Structures in the Linux Kernel](https://0xax.gitbooks.io/linux-insides/content/DataStructures/linux-datastructures-3.html) 並舉出 Linux 核心使用 bitmap 的案例
:::
Sign in with Wallet
Connect another wallet