# 2024q1 Homework4 (quiz3+4)
contributed by < [st10740](https://github.com/st10740) >
## [第三週測驗題](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-4)
### 測驗三
這段程式碼的目的是計算以 2 為底的對數,其做法是找到該數字的二進位表示中,最高位元為 1 的索引值,而該索引值即為以 2 為底的對數值。
然而,若是透過線性搜尋 (Linear search) 的方式會找太慢,因此以下方法採取二分搜尋法 (Binary search) 的方式進行搜索。首先會先利用 `while (i >= ...)` 確認左半邊的位元是否存在 1,若存在則搜尋左半邊,並將 `result` 加上一半的位元數,若左半邊不存在為 1 的位元,則搜尋右半邊,並不斷縮小範圍直到找到為止。
```c
static size_t ilog2(size_t i)
{
size_t result = 0;
while (i >= 65536) {
result += 16;
i >>= 16;
}
while (i >= 256) {
result += 8;
i >>= 8;
}
while (i >= 16) {
result += 4;
i >>= 4;
}
while (i >= 2) {
result += 1;
i >>= 1;
}
return result;
}
```
而另外一個作法的概念也相同,只是利用了 GCC 提供的內置函式 `__builtin_clz` ,其作用是計算無號整數 `v` 的二進位表示式中,從最高有效位開始的連續零位的数量。因此透過將 31 減去該值就可以得到最高位元為 1 的索引值。
```c
int ilog32(uint32_t v)
{
return (31 - __builtin_clz(v));
}
```
#### 在 Linux 核心原始程式碼中有使用 log2 的相關程式碼
在 [tools/include/linux/log2.h](https://elixir.bootlin.com/linux/latest/source/include/linux/log2.h) 中定義了巨集 `ilog2`。
這個方法會先利用 ` __builtin_constant_p` 判斷 `n` 在編譯時期是否為常數,若是則會在編譯時期完成如下面 `(n) & (1ULL << 63) ? 63` 的運算來求得最高位為 1 的索引值,因此 `ilog2` 的運算結果可以在編譯時期就確定;但若不是常數則需要等到執行時期再進行運算,並且會根據 `n` 的大小決定呼叫 `__ilog2_u32` 或是 `__ilog2_u64`。因此這個作法可以節省輸入為常數時的執行時間。
```c
#define ilog2(n) \
( \
__builtin_constant_p(n) ? ( \
(n) < 2 ? 0 : \
(n) & (1ULL << 63) ? 63 : \
(n) & (1ULL << 62) ? 62 : \
(n) & (1ULL << 61) ? 61 : \
(n) & (1ULL << 60) ? 60 : \
(n) & (1ULL << 59) ? 59 : \
(n) & (1ULL << 58) ? 58 : \
(n) & (1ULL << 57) ? 57 : \
(n) & (1ULL << 56) ? 56 : \
(n) & (1ULL << 55) ? 55 : \
// ...
(n) & (1ULL << 10) ? 10 : \
(n) & (1ULL << 9) ? 9 : \
(n) & (1ULL << 8) ? 8 : \
(n) & (1ULL << 7) ? 7 : \
(n) & (1ULL << 6) ? 6 : \
(n) & (1ULL << 5) ? 5 : \
(n) & (1ULL << 4) ? 4 : \
(n) & (1ULL << 3) ? 3 : \
(n) & (1ULL << 2) ? 2 : \
1 ) : \
(sizeof(n) <= 4) ? \
__ilog2_u32(n) : \
__ilog2_u64(n) \
)
```
我進行了簡單的測試確認以上巨集處理常數的部份是否真的會在編譯時期取得運算結果,因此我拿取了巨集中處理常數的程式碼建立了一個新的巨集,並在 `main` 中使用該巨集,將常數 `1022` 作為參數輸入,並將結果印出來。
```c
int main() {
printf("log_2_%d: %d\n", 1022, ilog2(1022));
return 0;
}
```
接著我將程式碼進行編譯,取得其編譯後的組合語言,可以觀察到第 10 行的 `9` 即為 `ilog2(1022)` 的結果,在這邊被當作常數來使用。
```c=
main:
.LFB0:
.cfi_startproc
endbr64
pushq %rbp
.cfi_def_cfa_offset 16
.cfi_offset 6, -16
movq %rsp, %rbp
.cfi_def_cfa_register 6
movl $9, %edx
movl $1022, %esi
leaq .LC0(%rip), %rax
movq %rax, %rdi
movl $0, %eax
call printf@PLT
movl $0, %eax
popq %rbp
.cfi_def_cfa 7, 8
ret
.cfi_endproc
```
### 測驗四
Exponentially Weighted Moving Average (EWMA; 指數加權移動平均) 是一種取平均的作法,在訊號處理中可以作為濾波器使用,其公式如下:
$$
S_t =
\begin{cases}
Y_0, & t = 0 \\
\alpha Y_t + (1 - \alpha) S_{t - 1}, & t > 0
\end{cases}
$$
可以透過 $\alpha$ 決定新資料與歷史資料的權重,介於 0 到 1 之間,$\alpha$ 值越大代表歷史資料的權重減少地越快。
在實作方面,首先宣告一個結構體 `ewma` ,`internal` 用來儲存當前算出來的平均值;`weight` 用來存新資料與歷史資料的的權重,$\frac{1}{weight}$ 等同於上面公式中的 $\alpha$ ;`factor` 的目的是保留平均數的小數部份,其值存的是欲保留的小數點位數,因為 `internal` 的型別為 unsigned long,進行平均值運算時可能會丟失小數部份,因此透過 `factor` 來避免誤差。
```c
struct ewma {
unsigned long internal;
unsigned long factor;
unsigned long weight;
};
```
`ewma` 的初始化確保傳入的 `weight` 和 `factor` 都是 2 的冪,原因是為了提升效能,會使用位元運算子進行乘除運算,接著分別取得 `weight` 和 `factor` 欲左移或右移的位元數,並將 `interval` 初始化為 0。
```c
void ewma_init(struct ewma *avg, unsigned long factor, unsigned long weight)
{
if (!is_power_of_2(weight) || !is_power_of_2(factor))
assert(0 && "weight and factor have to be a power of two!");
avg->weight = ilog2(weight);
avg->factor = ilog2(factor);
avg->internal = 0;
}
```
`ewma_add` 用以計算加入新資料後的整體平均值,當資料為第一筆時則直接把值當作平均值,即 $S_0 = Y_0$ ,對應到程式碼中的第 6 行,注意到這邊有對 `val` 進行左移 `avg->factor` 個位元的操作,目的是預先保留 `avg->factor` 個位元當作小數部分,避免之後進行平均值運算時丟失小數資訊。
而當資料不為第一筆時,則需計算新資料與歷史資料的權重,即 $S_t = \alpha Y_t + (1-\alpha)S_{t-1}$,對應到程式碼的第 4, 5 行,這邊針對新資料 `val` 也有先左移 `avg->factor` 個位元,另外觀察這行表示式 `(((avg->internal << avg->weight) - avg->internal) + (val << avg->factor)) >> avg->weight` 可以拆成:`((avg->internal << avg->weight) - avg->internal) >> avg->weight` 和 `(val << avg->factor) >> avg->weight` 相加,分別對應到 $(1-\alpha)S_{t-1}$ 和 $\alpha Y_t$ 的運算。
```c=
struct ewma *ewma_add(struct ewma *avg, unsigned long val)
{
avg->internal = avg->internal
? (((avg->internal << avg->weight) - avg->internal) +
(val << avg->factor)) >> avg->weight
: (val << avg->factor);
return avg;
}
```
因為 EWMA 的計算過程中為了保留小數部分有將新資料都先左移 `avg->factor` 個位元再進行運算,因此若要讀取計算完成的 EWMA,需要將 `avg->interval` 右移 `avg->factor` 個位元校正回來。
```c
unsigned long ewma_read(const struct ewma *avg)
{
return avg->internal >> avg->factor;
}
```
#### 在 Linux 核心原始程式碼中有使用 EWMA 的相關程式碼
在 [/include/linux/average.h](https://elixir.bootlin.com/linux/v6.8.2/source/include/linux/average.h) 中定義了巨集 `DECLARE_EWMA` 用以進行 `EWMA` 的運算。裡面定義的函式 `ewma_##name##_init`、`ewma_##name##_read` 和 `ewma_##name##_add` 都有進行以下的確認才進行後續的操作:
```c!
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);
```
首先透過 `BUILD_BUG_ON(!__builtin_constant_p(...))` 可以確保 `_precision` 和 `_weight_rcp` 在編譯時期都是常數。
`BUILD_BUG_ON((_precision) > 30)` 用以確保 `_precision` 的值不會大於 30,根據註解說明其目的是為了避免過多位元都被當作小數部分,不禁好奇為甚麼是使用 30 作為分界呢? 但是翻過歷史 commit 紀錄,都沒有說明使用 30 的原因。
另外我發現,在 Linux 核心原始程式碼中的 `_precision` 功能與測驗題的 `factor` 相同,但是比起 `factor`,`_precision` 的命名方式更可以讓我明白它的用途,而 [commit eb1e011](https://github.com/torvalds/linux/commit/eb1e011a14748a1d9df9a7d7df9a5711721a1bdb) 做的變更就是為了這個目的把 `_factor` 改成 `_precision`。
### 測驗五
測驗三的 `ilog2` 計算的是 $\lfloor log_{2}{(x)} \rfloor$ ,而測驗五的 `ceil_ilog2` 計算的是 $\lceil log_{2}{(x)} \rceil$。在作法上,`ilog2` 需要求該數字的二進位表示中,最高位元為 1 的索引值,其索引從 0 開始;而 `ceil_ilog2` 的作法則需要分成兩種情況探討,一種是其值為 2 的冪的情況,在這個情況下作法同 `ilog2`,另外一種則是值不為 2 的冪,其作法同 `ilog2` 但是最後需要再加上 1 。
一開始 `x--` 的操作用以處理 `x` 為 2 的冪的情況,使其不再是 2 的冪,因此便可以採取處理不是 2 的冪的方式進行計算,且最終也會得到正確的結果,因此接下來的操作適用於任何大於 0 的整數。
接下來的操作與 `ilog2` 找最高位元為 1 的方式大同小異,都是採取二分搜尋法,`r` 用來存最高位的索引值,不過與 `ilog2` 不同的地方在於,這邊都是使用位元運算子來進行操作,可以達到 branchless ,提高效能。
```c
int ceil_ilog2(uint32_t x)
{
uint32_t r, shift;
x--;
r = (x > 0xFFFF) << 4;
x >>= r;
shift = (x > 0xFF) << 3;
x >>= shift;
r |= shift;
shift = (x > 0xF) << 2;
x >>= shift;
r |= shift;
shift = (x > 0x3) << 1;
x >>= shift;
return (r | shift | x > 0x1) + 1;
}
```