# 2022q1 Homework2 (quiz2)
contributed by < [`Kevin-Shih`](https://github.com/Kevin-Shih) >
## 測驗 1 對二個無號整數取平均
`EXP1`:
> 考慮以下對二個無號整數取平均值的程式碼:
> ```c
> #include <stdint.h>
> uint32_t average(uint32_t a, uint32_t b)
> {
> return (a + b) / 2;
> }
> ```
>
> 考慮到整數溢位的因素，以下為第一種可行的實做:
> ```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 >> 1)`, `(b >> 1)` 先將 `a`、`b` 分別除以 2 再相加，避免了整數溢位的可能性，然而當 `a`、`b` 都是奇數時最終的結果會因為兩者在除 2 時各少了 0.5 (右移時均丟失最後 1 bit) 而比正確的平均值少 1 ，因此需要透過加上 `EXP1` 來補償。
有了這些瞭解就可以得出 `EXP1` 只有在 `a`, `b` 都是奇數，也就是 `a`, `b` 最低位均為 1 時為 1 ，其餘為 0 。因此 `EXP1` = `a & b & 1`，其中 `& 1` 用來取最後一個 bit 的結果。
`EXP2`, `EXP3`:
> 再次改寫為以下等價的實作:
> ```c
> uint32_t average(uint32_t a, uint32_t b)
> {
> return (EXP2) + ((EXP3) >> 1);
> }
> ```
> 其中 EXP2 和 EXP3 為 a 和 b 進行某種特別 bitwise 操作，限定用 OR, AND, XOR 這三個運算子。
看到 return 右半的 `(EXP3) >> 1` 可以把 `EXP3` 除 2，可以看出 `EXP3` 主要是做 `a + b` 的動作，又因題目限用 OR, AND, XOR 這三個運算子，其中只有 XOR 能做到半加器的功能，因此 `EXP3` = `a ^ b`。 接著因半加器缺少進位的動作，所以應該透過加上 `EXP2` 來補償，因進位只發生在 `a`, `b` 均為 1 的 bit 上，故應對 `a`, `b` 做 AND 取得需進位的 bit，且最終求的結果是 $(a + b) / 2$ 因此進位留在原本的 bit 就好不須左移 1 位 (與除 2 時的右移 1 位抵銷)，依此推論得到 `EXP2` = `a & b`。
### 比較上述實作在編譯器最佳化開啟時的差異
gcc 版本：
```
$ gcc --version
gcc (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
```
> gcc -S 所得到的 assembly 預設是 AT&T syntax
- ***(a + b) / 2***
以 `gcc -S -O3 q1_1.c` 得到第一種實作經編譯器最佳化後的組合語言輸出如下:
```assembly
// rdi: @a, rsi: @b
average:
.LFB0:
.cfi_startproc
endbr64
leal (%rdi,%rsi), %eax // eax = rdi + rsi
shrl %eax // eax >>= 1
ret
.cfi_endproc
```
第一種實作得到的結果，很直覺的將 `rdi`, `rsi` (`a`, `b`) 中的值相加放入 `eax` 後，將 `eax` 右移 1 位並返回該結果。 從 `rdi`, `rsi` 來看是以 64-bit 來做加法再取 lower 32 bits 存入 `eax`。
- ***low + (high - low) / 2***
以 `gcc -S -O3 q1_2.c` 得到第二種實作經編譯器最佳化後的組合語言輸出如下:
```assembly
// edi: @low, esi: @high
average:
.LFB0:
.cfi_startproc
endbr64
movl %esi, %eax // eax = esi
subl %edi, %eax // eax -= edi
shrl %eax // eax >>= 1
addl %edi, %eax // eax += edi
ret
.cfi_endproc
```
相較於第一種實做，全程採用 32-bit 方式存取暫存器，不像第一種以 64-bit 來做加法再取 lower 32 bits 可能導致整數溢位的問題。
- ***(a >> 1) + (b >> 1) + (a & b & 1)***
以 `gcc -S -O3 q1_3.c` 得到第三種實作經編譯器最佳化後的組合語言輸出如下:
```assembly=
// edi: @a, esi: @b
average:
.LFB0:
.cfi_startproc
endbr64
movl %edi, %eax // eax = edi
movl %esi, %edx // edx = esi
andl %esi, %edi // edi &= esi (a = a & b)
shrl %eax // eax >>= 1
shrl %edx // edx >>= 1
andl $1, %edi // edi &=1
addl %edx, %eax // eax += edx
addl %edi, %eax // eax += edi
ret
.cfi_endproc
```
基本上是以 `eax`, `edx`, `edi` 分別來做 `(a >> 1)`, `(b >> 1)`, `(a & b & 1)` ，但有個疑問是為何 #7 要將 `esi` 移到 `edx` ， #12 直接 `addl %esi, %eax` 應會得到相同的結果。
:::warning
TODO
釐清上述疑問
:::
雖改用 `Shift`, `AND` 操作但最佳化後的組合語言輸出反而比前一種多了 4 行，似乎更耗時了(多了 4 個指令)，有些出乎意料。
- ***(a & b) + ((a ^ b) >> 1)***
以 `gcc -S -O3 q1_4.c` 得到第四種實作經編譯器最佳化後的組合語言輸出如下:
```assembly
// edi: @a, esi: @b
average:
.LFB0:
.cfi_startproc
endbr64
movl %edi, %eax // eax = edi
andl %esi, %edi // edi &= esi
xorl %esi, %eax // eax ^= esi
shrl %eax // eax >>= 1
addl %edi, %eax // eax += edi
ret
.cfi_endproc
```
以 `eax`, `edi`分別來做 `a ^ b`, `(a ^ b) >> 1` ，最後在相加至 `eax` 並返回該值。
### Linux 核心原始程式碼 [include/linux/average.h](https://github.com/torvalds/linux/blob/master/include/linux/average.h) 的 EWMA 實作及 Linux 核心中的應用
EWMA 即 Exponentially weighted moving average ，計算移動平均時較舊的觀測量(數值)權值以指數形式遞減，但永遠不會到 0 ，其函數定義如下：
$S_t = \begin{cases}
Y_0, &t = 0\\
\alpha Y_t + (1-\alpha)\cdot S_{t-1}, &t > 0
\end{cases}$
> - $Y_{t}$ is the value at a time period $t$.
> - $S_{t}$ is the value of the EWMA at any time period $t$. -- Wiki: [EWMA](https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average)
接著看到 [include/linux/average.h](https://github.com/torvalds/linux/blob/master/include/linux/average.h) 中定義的 macro
`#define DECLARE_EWMA(name, _precision, _weight_rcp)`， 其三個參數的用途如下
- `name`: 定義一個對應的 struct ，名為 `ewma_##name##` ，以及對應的名稱的 `init`, `read`, `add` 函式
- `_precision`： 計算時 fraction part 要使用多少 bit ，不能大於 30
- `_weight_rcp`: 相當於公式中的 $1/\alpha$ ，為了效率考量須為 2 的冪次(可以用 shift 搞定)
DECLARE_EWMA 所定義的 struct
```c
struct ewma_##name {
unsigned long internal;
};
```
其實只有一個 unsigned long 變數 `internal` 用於儲存 EWMA 的結果。
接著看到第一個函式 `ewma_##name##_init`
```c
static inline void ewma_##name##_init(struct ewma_##name *e)
{
BUILD_BUG_ON(!__builtin_constant_p(_precision));
BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp));
/*
* Even if you want to feed it just 0/1 you should have
* some bits for the non-fractional part...
*/
BUILD_BUG_ON((_precision) > 30);
BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp);
e->internal = 0;
}
```
前 4 行的 `BUILD_BUG_ON` 均是當輸入的參數不符合要求時可以在 compile 階段就報錯，以保證執行階段時的效率與正確性。 最後就是將 `internal` 初始化為 0 ，計算 EWMA 前應先呼叫此函式以確保 `internal` 初始化為 0 。
> `__builtin_constant_p(x)` 用於檢驗是 `x` 否在 compile time 就已知為常數，若否則中斷 compile 並報錯。 -- [gcc onlinedocs](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html)
> `BUILD_BUG_ON(condition)` break compile if a condition is true. -- [include/linux/build_bug.h](https://github.com/torvalds/linux/blob/master/include/linux/build_bug.h)
第二個函式 `ewma_##name##_read`
```c
static inline unsigned long ewma_##name##_read(struct ewma_##name *e)
{
BUILD_BUG_ON(!__builtin_constant_p(_precision));
BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp));
BUILD_BUG_ON((_precision) > 30);
BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp);
return e->internal >> (_precision);
}
```
再度出現同樣 4 行的 `BUILD_BUG_ON` 有點怪，似乎寫一遍即可。
:::warning
TODO
find out why
:::
讀出 `e->internal` 並右移 `_precision` 位，已得到整數位的 EWMA (因 `internal` 保留後 `_precision` bit 作為運算時的 fraction part )。
第三個函式 `ewma_##name##_add`
```c=
static inline void ewma_##name##_add(struct ewma_##name *e,
unsigned long val)
{
unsigned long internal = READ_ONCE(e->internal);
unsigned long weight_rcp = ilog2(_weight_rcp);
unsigned long precision = _precision;
BUILD_BUG_ON(!__builtin_constant_p(_precision));
BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp));
BUILD_BUG_ON((_precision) > 30);
BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp);
WRITE_ONCE(e->internal, internal ?
(((internal << weight_rcp) - internal) +
(val << precision)) >> weight_rcp :
(val << precision));
}
```
計算 EWMA 的主體，用於計算新增一筆資料 `val` 後的新 EWMA 並存入 `e->internal` 中。
- #4, #13 使用了 `READ_ONCE`, `WRITE_ONCE` 應該是為了保護在多個 thread 的讀寫操做產生 data race 。
- #5 以 `ilog2(_weight_rcp)` 求得稍後要位移 `weight_rcp` 個 bit 。
- #6 先將 `_precision` 的值存入 `precision` 主要是避免 macro 展開後可能出現的 double evaluation 的情形。
#13 當 `internal` 為 0 時， `e->internal` = `val << precision` 符合 EWMA 定義中 $t = 0$ 的情形。 而當 `internal` 不為 0 時，其敘述相當於 $internal = (internal (2^{weight\_rcp}-1)+val)/2^{weight\_rcp}$
> $2^{weight\_rcp}$ = `_weight_rcp`, `1/_weight_rcp` = $\alpha$
符合 EWMA 定義中 $t > 0$ 的情形。
**EWMA 在 Linux 核心中的應用**
直接查找 `DECLARE_EWMA` 關鍵字在 kernel 中出現的地方，之前修「通訊網路」時在計算 RTT 時出現過 EWMA ，預期 Linux 核心中也會有類似的應用。
Linux 核心中的應用:
- [net/mac80211/ieee80211_i.h](), line 435
- [net/mac80211/sta_info.h](https://github.com/torvalds/linux/blob/master/net/mac80211/sta_info.h), line 126, 373, 374, 433
- 還有不少也是在 drivers/net/wireless 底下的應用...
在 sta_info.h 中 EWMA 被用於計算 TX bit rate(傳輸的 bit rate), failed MSDUs(link layer 的 payload) 的百分比等用途，。
## 測驗 2 不需要分支的 max
`EXP4`, `EXP5`:
> 考慮以下「不需分支的 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 個運算子。
>
題目為了取最大值，因此當 `a` 較大的時候回傳值應等價於 `a ^ 0` = `a`。
而當 `b` 較大時回傳值應等價於 `a ^ (a ^ b)` = `0 ^ b` = `b`
因此得出 `EXP4` 應為 `a ^ b`，且 `-(EXP5)` 在 a > b 時為 0，在 a < b 時為 -1 即 0xFFFFFFFF，其中 0xFFFFFFFF 是透過 `-` 對 1 取 2 的補數來的，因此 `EXP5` 在 a > b 時應為 0，在 a < b 時應為 1，故 `EXP5` = `a < b`。
### 針對 32 位元有號整數，撰寫同樣 branchless 的實作
```c
#include <stdint.h>
int32_t max(int32_t a, int32_t b)
{
return a ^ ((a ^ b) & -(a < b));
}
```
將輸入的 parameters 改為 `int32_t` 型態就能以相同的方式的到正確的 max 結果了。
### 在 Linux 核心原始程式碼中，找出更多類似的實作手法
[arch/x86/kvm/mmu.h](https://github.com/torvalds/linux/blob/master/arch/x86/kvm/mmu.h) 中：
```c
/*
* If CPL < 3, SMAP prevention are disabled if EFLAGS.AC = 1.
*
* If CPL = 3, SMAP applies to all supervisor-mode data accesses
* (these are implicit supervisor accesses) regardless of the value
* of EFLAGS.AC.
*
* This computes (cpl < 3) && (rflags & X86_EFLAGS_AC), leaving
* the result in X86_EFLAGS_AC. We then insert it in place of
* the PFERR_RSVD_MASK bit; this bit will always be zero in pfec,
* but it will be one in index if SMAP checks are being overridden.
* It is important to keep this branchless.
*/
unsigned long smap = (cpl - 3) & (rflags & X86_EFLAGS_AC);
int index = (pfec >> 1) +
(smap >> (X86_EFLAGS_AC_BIT - PFERR_RSVD_BIT + 1));
```
git log 中還能找到一些其他的實做，但有些似乎已經消失了，在目前 github 上的 repo 找不到了
- [Replace int2float() with an optimized version.](https://github.com/torvalds/linux/commit/747f49ba67b8895a5831ab539de551b916f3738c)
- [Optimize pte permission checks](https://github.com/torvalds/linux/commit/97d64b788114be1c4dc4bfe7a8ba2bf9643fe6af)
## 測驗 3 64 位元 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;
> while (v) {
> uint64_t t = v;
> v = u % v;
> u = t;
> }
> return u;
> }
> ```
> 可以改寫為以下等價實做:
> ```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;
> }
> ```
根據題目中提到的 GCD 演算法，當 `u`, `v` 任一者為 0 時， gcd 即為非零者的值，當均為 0 時，gcd 為 0 。 這三種情形透過 #4 的 if 敘述來處理。
接著 #6 的 for 迴圈會在 `u`, `v` 均為偶數時，將其同除以 2 ，並將除的次數紀錄在 `shift` 以便最後乘回 gcd 中(左移 `shift` 次)。 #9 的 while 迴圈則是當 `u` 還是偶數時，除 2 直到 `u` 變為奇數。
接著 #11 開始的 do while 就是輾轉相除法的本體了，當 `u` 小於 `v` 就將 `v` 減去 `u` ，大於的話就將 `u - v` 存到 `v` ， `v` 存到 `u` ，直到其中一者為 0 為止，然而由於這個實做中，相減後的結果都會存到 `v` ，因此中止條件 `COND` 只須判斷 `v` 是否為 0 即可，故 `COND` = `v` 。
最終的回傳值 `RET` 則是將先前輾轉相除法最後留下的 `u` 值左移 `shift` 位(先前同除以 2 的次數)得到答案，因此 `RET` = `u << shift` 。
### 在 x86_64 上透過 __builtin_ctz 改寫 GCD，分析對效能的提升
**透過 __builtin_ctz 改寫的 gcd64_ctz**
```c
uint64_t gcd64_ctz(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
int shift = __builtin_ctz(u | v);
u >>= __builtin_ctz(u);
do {
v >>= __builtin_ctz(v);
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u << shift;
}
```
將原先同除以 2 的部份改為紀錄 `u`, `v` 中最低位的 1 作為 `shift` 。
```c
while (!(u & 1))
u /= 2;
```
取代為 `u` 直接右移到無法被 2 整除。 `v` 的部份則留到迴圈內處理。
其餘部份與原先的 `gcd64` 相同。
**效能測試方式**
先以 random 方式準備測試資料
```c
void prepare_data(uint64_t *data1, uint64_t *data2, size_t size){
srand(time(NULL));
for (size_t i = 0; i < size; i++) {
data1[i] = rand();
data2[i] = rand();
}
}
```
兩者分別測試
```c
prepare_data(data1, data2, size);
clock_t start = clock();
for (size_t i = 0; i < size; i++)
gcd64(data1[i], data2[i]);
clock_t finish = clock();
clock_t start2 = clock();
for (size_t i = 0; i < size; i++)
gcd64_ctz(data1[i], data2[i]);
clock_t finish2 = clock();
double time1 = (double)(finish - start) / CLOCKS_PER_SEC;
double time2 = (double)(finish2 - start2) / CLOCKS_PER_SEC;
```
在本機連續執行 10,000,000 次、 3 round 的結果如下
(5 round 取中間 3 round)
||gcd64|gcd64_ctz|
|-|-|-|
| 1 | 2.125378 secs | 0.936948 secs |
| 2 | 2.126757 secs | 0.936914 secs |
| 3 | 2.121405 secs | 0.935645 secs |
| avg.| 2.124513 secs | 0.936502 secs |
採用 ctz 的實做 3 round 的平均快了約 2.268562 倍，顯示在測試環境中 `ctz` 顯然比除法運算快上不少。
### 請解釋 Linux 核心中內建 GCD (不只一種實作)，的實作手法並探討在核心內的應用場景。
**[lib/math/gcd.c](https://github.com/torvalds/linux/blob/master/lib/math/gcd.c) 中的 GCD 實作手法**
```c
unsigned long gcd(unsigned long a, unsigned long b)
{
unsigned long r = a | b;
if (!a || !b)
return r;
b >>= __ffs(b);
if (b == 1)
return r & -r;
for (;;) {
a >>= __ffs(a);
if (a == 1)
return r & -r;
if (a == b)
return a << __ffs(r);
if (a < b)
swap(a, b);
a -= b;
}
}
```
當 ffs 可用時會採用以上實做，當 `a`, `b` 中一者為 0 或均為 0 時，gcd 為 `a | b` (三種情形分別對應 `b`, `a`, `0`)。
若 `b` 是 2 的冪次，則回傳 `r & -r`。 `r & -r` 的作法是透過 2 補數特性找出 `r` 中最低位的 1 並使其他 bit 為 0 (就像測驗 4 中 `bitset & -bitset`)，而由於 `r` = `a | b` 因此相當於找出 `min(a & -a, b & -b)` ，也就是兩者的最大公因數(因 `b` 是 2 的冪次，兩者 gcd 是小於等於 `b` 的 2 的冪次)。
若 `b` 不是 2 的冪次， `b >>= __ffs(b)` 也會如 gcd 算法中所述，將 `b` 除到剩下奇數的狀態。
for 迴圈內有 4 步：
- 將 `a` 除到剩下奇數的狀態。
- 若 `a` 是 2 的冪次，若是則 gcd 為 `r & -r` (同 `b` 的情形，代表沒有 2 的冪次以外的公因數)。
- 若剩下的 `a == b` 代表剩下的 `a` 也是兩者的公因數，將 `a` 左移先前同除的 2 的次數(`ffs(r)` = `min(ffs(a), ffs(b))`)即為 gcd 。 (再減一次就是輾轉相除法，其中一個減到變 0 的情形)
- 最後當上述情形都不成立時，執行輾轉相除法中相減的步驟，若 a < b 就先交換 `a`, `b`，確保 `a` 是被減的那個。
**核心內的應用場景**
- [block/blk-settings.c](https://github.com/torvalds/linux/blob/master/block/blk-settings.c)
- [drivers/media/tuners/mt2063.c](https://github.com/torvalds/linux/blob/master/drivers/media/tuners/mt2063.c)
- [sound/core/pcm_timer.c](https://github.com/torvalds/linux/blob/master/sound/core/pcm_timer.c)
- ...
其中在 `mt2063.c` 中一個用於檢測 [LO](https://en.wikipedia.org/wiki/Local_oscillator) spur 是否發生的函數 `IsSpurInBand()` 中出現，用於計算某種 scale factor 。
```c
/*
** For each edge (d, c & f), calculate a scale, based on the gcd
** of f_LO1, f_LO2 and the edge value. Use the larger of this
** gcd-based scale factor or f_Scale.
*/
lo_gcd = gcd(f_LO1, f_LO2);
gd_Scale = max((u32) gcd(lo_gcd, d), f_Scale);
hgds = gd_Scale / 2;
gc_Scale = max((u32) gcd(lo_gcd, c), f_Scale);
hgcs = gc_Scale / 2;
gf_Scale = max((u32) gcd(lo_gcd, f), f_Scale);
hgfs = gf_Scale / 2;
```
另外在 `git log --grep` 則發現不少 commit 是捨棄 local 的 gcd implementation 而改用 [lib/math/gcd.c](https://github.com/torvalds/linux/blob/master/lib/math/gcd.c) 的實作，例如剛剛提到的 sound/core/pcm_timer.c 在 [commit a96053](https://github.com/torvalds/linux/commit/a9605391cfab2bf9a73e51faac5147622f73c6d5) 時捨棄。
## 測驗 4 bitset
> 在影像處理中，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;
> }
> 以 __builtin_ctzll 改寫的程式碼如下:
> ```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]; // a row of bitmap per iteration
> while (bitset != 0) { // till this row has no more `set bit`
> 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 的特性)。
題目提到 `EXP6` 應保留 `bitset` 中最低位的 1 而其餘均為 0 ，又提示要 `-bitset` 的特性，以 bitset 為 $010100_b$ 為例， `-bitset` 為 $101100_b$ ，兩者只有原先最低位的 1 的 bit 均為 1 ，因此善用 2 的補數特性， `EXP6` = `bitset & -bitset`。
`-bitset` 特性:
| | binary |
|---------|------------|
| bitset | $010100_b$ |
| ~bitset | $101011_b$ |
| -bitset | $101100_b$ |
`-bitset` = `~bitset + 1` 因此只有最低位的 1 與 bitset 均為 1 。
### 設計實驗，檢驗 ctz/clz 改寫的程式碼相較原本的實作有多少改進？應考慮到不同的 bitmap density
**準備測試資料的函式**
```c
void prepare_data(uint64_t *bitmap, size_t size, int density){
srand(time(NULL));
// control bitmap density with 0 ~ 4 level 0%, 25%, 50%, 75%, 100%
int count = 16 * size * density;
// 100% density
if (density == 4){
for (size_t i = 0; i < size; i++)
bitmap[i] = 0xFFFFFFFFFFFFFFFFul;
return;
}
// init bitmap
for (size_t i = 0; i < size; i++)
bitmap[i] = 0;
// 0% density
if (density == 0)
return;
// select `count` (pseudo) random `0` bits set them to `1`
while (count) {
int i = rand64() % (size * 64);
if (bitmap[i / 64] & 1ul << (i % 64))
continue;
else{
bitmap[i / 64] |= 1ul << (i % 64);
count--;
}
}
}
```
**測試方式**
```c
size_t bitmapsize = strtoul(argv[1], NULL, 10);
uint64_t *bitmap = (uint64_t *)malloc(bitmapsize * sizeof(uint64_t));
uint32_t *out_naive = (uint32_t *)malloc(bitmapsize * 64 * sizeof(uint32_t));
uint32_t *out_improved = (uint32_t *)malloc(bitmapsize * 64 * sizeof(uint32_t));
for(size_t i = 0; i < 5; i++){ // 5 density level 0%, 25%, 50%, 75%, 100%
prepare_data(bitmap, bitmapsize, i);
for (size_t j = 0; j < bitmapsize * 64; j++)
out_naive[j] = out_improved[j] = 0;
clock_t start = clock();
naive(bitmap, bitmapsize, out_naive);
clock_t finish = clock();
clock_t start2 = clock();
improved(bitmap, bitmapsize, out_improved);
clock_t finish2 = clock();
double time1 = (double)(finish - start) / CLOCKS_PER_SEC;
double time2 = (double)(finish2 - start2) / CLOCKS_PER_SEC;
printf("density: %3ld%%, naive: %f secs, improved: %f secs\n", i * 25, time1, time2);
}
free(bitmap);
free(out_naive);
free(out_improved);
```
測試 5 種不同的 bitmap density 的情形下 `naive` 與 `ctz` 改寫的版本的執行時間。
在我的環境中執行 `bitmapsize` = 1,000,000 的結果:
| density | naive (secs) | improved (secs) |
|-|-|-|
| 0% | 0.129309 | 0.002487 |
| 25% | 0.276713 | 0.069022 |
| 50% | 0.390655 | 0.121114 |
| 75% | 0.273711 | 0.174631 |
| 100% | 0.144110 | 0.229330 |
在絕大多數時間 `ctz` 改寫的 `improved` 版本確實比 `naive` 快上不少，但有趣的是在 100% density (all set) 的情形中 `naive` 反而比 `improved` 快上一些。
兩個版本的主要差異：
**naive**
```c
for (int i = 0; i < 64; i++) {
if ((bitset >> i) & 0x1)
out[pos++] = p + i;
}
```
**improved**
```c
while (bitset != 0) { // till this row has no more `set bit`
uint64_t t = EXP6;
int r = __builtin_ctzll(bitset);
out[pos++] = k * 64 + r;
bitset ^= t;
}
```
在 100% density 兩者的內層迴圈都須 iterate 64 次且 `naive` 的 if 條件恆為真，因此 prediction miss 幾乎不會出現，因此迴圈內單純對 `out[pos++]` 賦值的 `naive` 反而比較快。
另外， `naive` 版本在 50% density 時的執行時間最久，應是因為 prediction miss 較高， `improved` 版本則是隨著 density 上升執行時間隨之遞增 (`ctz` 所能省下內層迴圈 iteration 次數越來越少)。
### 思考進一步的改進空間
```c=
size_t improved_v2(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) {
uint64_t t = bitset & -bitset;
out[pos++] = k << 6 | __builtin_ctzll(bitset);
bitset ^= t;
}
}
return pos;
}
```
#9 中避免使用 `*`, `+` 運算改用 `<<`, `|` ，並將原本 `__builtin_ctzll(bitset)` 不再存入 `r` 而直接移至 #9 使用讓程式更簡潔。
:::spoiler **測試程式**
```c
int main(int argc,char **argv){
if (argc != 2)
return 1;
size_t bitmapsize = strtoul(argv[1], NULL, 10);
uint64_t *bitmap = (uint64_t *)malloc(bitmapsize * sizeof(uint64_t));
uint32_t *out_naive = (uint32_t *)malloc(bitmapsize * 64 * sizeof(uint32_t));
uint32_t *out_improved = (uint32_t *)malloc(bitmapsize * 64 * sizeof(uint32_t));
uint32_t *out_improved_v2 = (uint32_t *)malloc(bitmapsize * 64 * sizeof(uint32_t));
for(size_t i = 0; i < 5; i++){ // 5 density level 0%, 25%, 50%, 75%, 100%
prepare_data(bitmap, bitmapsize, i);
for (size_t j = 0; j < bitmapsize * 64; j++)
out_naive[j] = out_improved[j] = out_improved_v2[j] = 0;
clock_t start = clock();
naive(bitmap, bitmapsize, out_naive);
clock_t finish = clock();
clock_t start2 = clock();
improved(bitmap, bitmapsize, out_improved);
clock_t finish2 = clock();
clock_t start3 = clock();
improved_v2(bitmap, bitmapsize, out_improved_v2);
clock_t finish3 = clock();
// test output correctness
int equal = 1;
for (int i = 0; i < bitmapsize * 64; i++){
if (out_naive[i] != out_improved[i]) {
printf("improved failed\n");
equal = 0;
break;
}
if (out_naive[i] != out_improved_v2[i]) {
printf("improved_v2 failed\n");
equal = 0;
break;
}
}
double time1 = (double)(finish - start) / CLOCKS_PER_SEC;
double time2 = (double)(finish2 - start2) / CLOCKS_PER_SEC;
double time3 = (double)(finish3 - start3) / CLOCKS_PER_SEC;
printf("density: %3ld%%, naive: %f secs, improved: %f secs, improved_v2: %f secs\n", i * 25, time1, time2, time3);
}
free(bitmap);
free(out_naive);
free(out_improved);
free(out_improved_v2);
return 0;
}
```
:::
(與前面測試 `ctz` 版本的方式相同，但加入檢查 output 正確性的程式，避免 `improved_v2` 改出 bug )
在我的環境中執行 `bitmapsize` = 1,000,000 的結果如下
| density | naive (secs) | improved (secs) | improved_v2 (secs) |
|-|-|-|-|
| 0% | 0.135698 | 0.002625 | 0.002510 |
| 25% | 0.291539 | 0.071203 | 0.067933 |
| 50% | 0.395848 | 0.121448 | 0.117784 |
| 75% | 0.274069 | 0.173835 | 0.169411 |
| 100% | 0.143882 | 0.213872 | 0.206336 |
**聚焦在 `improved`, `improved_v2` 的折線圖**
從折線圖可以看到避免使用 `*`, `+` 而改用 bitwise operator `<<`, `|` 可以節省一小部份的執行時間，而且在 density 較高時比較明顯。
### 閱讀 [Data Structures in the Linux Kernel](https://0xax.gitbooks.io/linux-insides/content/DataStructures/linux-datastructures-3.html) 並舉出 Linux 核心使用 bitmap 的案例
- drivers/net/ipa/ipa_endpoint.h : [commit c1aaa0](https://github.com/torvalds/linux/commit/c1aaa01dbf4cef95af3e04a5a43986c290e06ea3)
- net/mptcp/pm_netlink.c : [commit 86e39e](https://github.com/torvalds/linux/commit/86e39e04482b0aadf3ee3ed5fcf2d63816559d36)
- [sound/usb/clock.c #260](https://github.com/torvalds/linux/blob/16f73eb02d7e1765ccab3d2018e0bd98eb93d973/sound/usb/clock.c#L260)
- ...
大多數是用來儲存一種狀態如 flags, visted 等等，相較使用 bool 佔用 1 byte 儲存 1 個狀態 (true / false)， bitmap 1 byte 最多儲存 8 個狀態，若是需要節省空間， bitmap 是很好的選擇。
> [include/linux/types.h](https://github.com/torvalds/linux/blob/master/include/linux/types.h) 中 bitmap 是以 unsigned long array 方式宣告的，因此若要存 65 個狀態就需要 128 bits 似乎也是會浪費不少空間 (但同時也保留了不少未來擴充更多狀態的空間)
Sign in with Wallet
Connect another wallet