# Linux 核心專題: 重作第 3, 4, 7 週測驗題
> 執行人: dockyu
> [專題解說影片](https://youtu.be/wHDFin4Sr7w)
### Reviewed by `Daniel-0224`
假設目標是求 $\sqrt{x}=N, x=N^2$
$N$ 可以視為
$$
N=a_n+ a_{n-1}+ a_{n-2}+ \dots+ a_2+ a_1+ a_0\text{ , }a_i=2^i \text{ or } a_i=0
$$
$N=\sum_{i=0}^n a_{i}$
多打了 `$$`。
### Reviewed by `yu-hsiennn`
\begin{split}
\left[
\begin{array}{c}
a_na_n \\
\end{array}
\right] = a_{n}\left[ 2\left(\right)+a_{n}\right]
\end{split}
這個 2 是? 乘上 1 還是 0? 如果不看後面的推導會不知道這個等式怎麼成立的。
> 這裡的 2 要乘上 0,我已經把 0 加上去了,我的原意是想讓讀者了解每一個 L 型都可以寫成這個形式
> 顯然 $P_m = P_{m+1} + a_m$
>
因為前方寫的是 $P_m = a_n+a_{n-1}+...+a_m$ ,以直覺來說,會有點卡住為什麼會成立(通常數字會遞增呈現),或許可以加一些敘述,如 $P_m$ 指的是由 $m$ 到 $n$ 的總和。
> 現在這樣就夠清楚了
第 3 週測驗 5 的程式碼,若將 1 帶入,結果會為 1,但是 `ceil_ilog(1)` 應為 0。
> 我只依照 extra 的說明處理了輸入為 0 的情況,經過你的提醒之後才發現題目給的程式碼並不能處理輸入 1 的情況,可以將目前遇到的問題簡化為: 在 $x=1$ 時結果要減 1,換句話說就是最後一行的 `+1` 不要做就可以了,所以如果將 `+1` 替換為 `+(x != 1)` 就只會在 $x=1$ 時 `+0` (比原本的版本少 1),可以在輸入 1 時正確回傳 0,完整程式碼已經更新在 [測驗五](#%E6%B8%AC%E9%A9%97-5)
```diff
int ceil_ilog2(int x)
{
- return ((r | shift | u > 1) + 1 | z;
+ return ((r | shift | u > 1) + (x != 1)) | z;
}
```
### Reviewed by `marsh-fish`
建議在探討 SWAR 的時間可以把相關的函式獨立計算時間,比較好看出時間差異。
### Reviewed by `ChenFuhuangKye`
> $x/10$ 等同於 $x \times \frac{1}{10}$ 約等於 $x \times \frac{13}{128}$ ,之所以使用 13 跟 128 是因為 除以 128 可以透過位移做到,乘以 13 也可以透過 bitwise 做到,且這個數字在 $1\le x \le 19$ 時的誤差在小數點後一位之下,當然也有其他數字更接近 $\frac{1}{10}$ 不過題目已經確定被除數最大為 19,所以 $\frac{13}{128}$ 已經足夠
被除數**最大**為 19 有錯字要更正,可以列出 0 ~ 19 的誤差嗎?
>錯字已更正
>數字 $x$ 的誤差會是 $x*(\cfrac{13}{128}-\cfrac{1}{10})=x*\cfrac{1}{640}$ ,當 $x=19$ 時誤差為 $19*\cfrac{1}{640}=0.0296875$ ,若 $x<19$ 則誤差會比 $x=19$ 來的小
### Reviewed by `david965154`
>```c
>while (!IS_ALIGNED((long)src, sizeof(long))) {
> if (!length--)
> return NULL;
> if (*src == d)
> return (void *)src;
> src++;
> }
>```
此處用途為何,為何只在前方做判斷,能保證記憶體位置的後方無需做此種處理嗎?
## 任務簡介
重作[第 3 週測驗題](https://hackmd.io/@sysprog/linux2024-quiz3)、[第 4 週測驗題](https://hackmd.io/@sysprog/linux2024-quiz4), [第 7 週測驗題](https://hackmd.io/@sysprog/linux2024-quiz7)和延伸問題,探究 bitops 在 Linux 核心的應用。
## TODO: [第 3 週測驗](https://hackmd.io/@sysprog/linux2024-quiz3)
> 包含延伸問題
### 測驗 1
[第 3 週測驗題-測驗 1](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-1)
#### 版本一
版本一對 $N$ 取 $log_2$ , 參考以下的表格可以發現其實 $\log_2 N$ 就是在取 $N$ 的 most significant bit (MSB)
| $N$ | $N$ 十六進位 | $\log_2 N$=msb |`1 << msb` |
|-|-|-|-|
|1|`0001`|0|`0001`|
|2|`0010`|1|`0010`|
|3|`0011`|1|`0010`|
|4|`0100`|2|`0100`|
|6|`0110`|2|`0100`|
|8|`1000`|3|`1000`|
從上面的表格可以看出再開平方根時完全不需要考慮 `1 << msb` 以上的位元,所以只要<s>遍歷</s>逐一走訪比 `1 << msb` 小的所有可能,版本一的方法是每次都將 `(1 << msb) >> 1` 加上去,直到變為0
:::danger
注意用語,詳見: https://hackmd.io/@sysprog/it-vocabulary
:::
###### 版本一 程式碼
```c
int msb = (int) log2(N);
int a = 1 << msb;
int result = 0;
while (a != 0) {
if ((result + a) * (result + a) <= N)
result += a;
a >>= 1;
}
return result;
```
#### 版本二
因為 $\log_2 N$ 的目的是得到 msb,版本一使用 math library 裡的 log function,但是也可以透過 bitwise 做到
###### 版本二 程式碼
```c
int msb = 0;
int n = N;
while (n > 1) {
n >>= 1;
msb++;
}
```
> [你所不知道的 C 語言:數值系統](https://hackmd.io/@sysprog/c-numerics#Count-Leading-Zero)中也有提到,當我們計算 (以 2 為底的對數) 時, 其實只要算高位有幾個 0's bits. 再用 31 減掉即可
```c
int BITS = 31;
while (BITS) {
if (N & 0x80000000) break;
N <<= 1;
BITS--;
}
```
所以$\log_2N$也可以改寫成以上形式
---
假設目標是求 $\sqrt{x}=N$,$x=N^2$,$N$ 可以視為
$$N=a_n+ a_{n-1}+ a_{n-2}+ \dots+ a_2+ a_1+ a_0\text{ , }a_i=2^i \text{ or } a_i=0$$
$$N=\sum_{i=0}^n a_{i}$$
$$
\begin{split}
N^2
&= N*N \\
&= \left(\sum_{i=0}^n a_{i}\right)*\left(\sum_{i=0}^n a_{i}\right) \\
&= \left(a_n+ a_{n-1}+ a_{n-2}+ \dots+ a_2+ a_1+ a_0 \right)*\left(a_n+ a_{n-1}+ a_{n-2}+ \dots+ a_2+ a_1+ a_0 \right) \\
&= \left[
\begin{array}{ccccc}
a_na_n & a_na_{n-1} & a_na_{n-2} & ... & a_na_{1} & a_na_0 \\
a_{n-1}a_n & a_{n-1}a_{n-1} & a_{n-1}a_{n-2} & ... & a_{n-1}a_1 & a_{n-1}a_0 \\
a_{n-2}a_n & a_{n-2}a_{n-1} & a_{n-2}a_{n-2} & ... & a_{n-2}a_1 & a_{n-2}a_0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
a_1a_n & a_1a_{n-1} & a_1a_{n-2} & ... & a_1a_1 & a_1a_0 \\
a_0a_n & a_0a_{n-1} & a_0a_{n-2} & ... & a_0a_1 & a_0a_0 \\
\end{array}
\right]
\end{split}
$N^2$ 是上面的矩陣的每一項相加如果將矩陣拆如下可以找到規律
\begin{split}
\left[
\begin{array}{c}
a_na_n \\
\end{array}
\right] = a_{n}\left[ 2\left( 0 \right)+a_{n}\right]
\end{split}
\begin{split}
\left[
\begin{array}{cc}
& a_na_{n-1} \\
a_{n-1}a_n & a_{n-1}a_{n-1}
\end{array}
\right] = a_{n-1}\left[ 2\left(a_n\right)+a_{n-1}\right]
\end{split}
\begin{split}
\left[
\begin{array}{ccc}
& & a_na_{n-2} \\
& & a_{n-1}a_{n-2} \\
a_{n-2}a_n & a_{n-2}a_{n-1} & a_{n-2}a_{n-2}
\end{array}
\right] = a_{n-2}\left[ 2\left(a_n+a_{n-1}\right)+a_{n-2}\right]
\end{split}
\begin{split} \\
N^2 =& a_n^2+[2a_n+a_{n-1}]a_{n-1}+[2(a_n+a_{n-1})+a_{n-2}]a_{n-2}+...+[2(\displaystyle\sum_{i=1}^{n}a_i)+a_0]a_0 \\
=& a_n^2+[2P_n+a_{n-1}]a_{n-1}+[2P_{n-1}+a_{n-2}]a_{n-2}+...+[2P_1+a_0]a_0 \\
P_m =& a_n+a_{n-1}+...+a_m \\
\end{split}
顯然 $P_m = P_{m+1} + a_m$
目標是求 $P_0$, 但是 $P_m$ 會由前一個得到,所以要先知道 $P_n$ -> $P_{n-1}$ -> $P_{n-2}$ -> $\dots$ -> $P_0$,而 $P_n$ 也是由 $P_{n+1}$ 得到的,所以可以初始化 $P_{n+1}=0$ 因為從 $P$ 的定義裡 $P_{n+1}$ 不會是任何數值的相加
每一次都判斷 $P_m^2 \leq N^2$ 是否成立,也就是版本一&版本二的每一次迴圈做的事
\begin{cases}
P_m = P_{m+1} + 2^m, & \text{if $P_m^2 \leq N^2$} \text{(第m個 bit 是 1)} \\
P_m = P_{m+1}, & \text{otherwise} \text{(第m個 bit 是 0)}
\end{cases}
但是每一次都計算 $P_m^2 \leq N^2$ 太麻煩,我們只需要判斷 $X_m = N^2 - P_m^2$ 有沒有大於 0 就可以了
+ $X_m = N^2 - P_m^2 = X_{m+1} - Y_m$
+ $Y_m = P_m^2 - P_{m+1}^2 = 2P_{m+1}a_m+a_m^2$
每一輪都只需要 計算 $Y_m$ & $X_m$ 並判斷 $X_m$ 就可以了,需要紀錄 $P_m$ 做為下一輪的 $P_{m+1}$
![圖片](https://hackmd.io/_uploads/HJ9QhiTVA.png)
![圖片](https://hackmd.io/_uploads/ryD9t3pER.png)
+ $Y_m = P_m^2 - P_{m+1}^2 = 2P_{m+1}a_m+a_m^2$
+ $c_m=2P_{m+1}a_m=P_{m+1}(a_m*2)=P_{m+1}(2^m*2)=P_{m+1}2^{m+1}$
+ $d_m=(a_m)^2=(2^m)^2$
\begin{gather*}
Y_m = \begin{cases}
c_m + d_m \ \ ,if \ \ a_m=2^m\\
0 \ \ ,if \ \ a_m=0\\
\end{cases}
\end{gather*}
可以用 $c_m$, $d_m$ 用 bitwise 快速求出 $c_{m+1}$, $d_{m-1}$
+ $c_{m-1}=P_m2^m=(P_{m+1}+a_m)2^m=P_{m+1}2^m+a_m2^m=\begin{cases}
c_m/2+d_m & \text{if }a_m=2^m \\
c_m/2 & \text{if }a_m=0
\end{cases}$
+ $d_{m-1}=\dfrac{d_m}{4}$
初始化 `n`
+ $X_{n+1}$: $N^2$
+ $c_n$: $0$
+ $d_n$: $(2^n)^2$ ,這裡的求法是讓 $d_n<N^2$ 也就是 $4^n<N^2$,所以 1~3 是 0,4~15 是 1,16~63 是 2,發現在二進位最高位每移動兩位 n 就加 1
$c_{-1}=P_0$ 即為所求
```c
int i_sqrt(int x)
{
if (x <= 1) /* Assume x is always positive */
return x;
int z = 0;
for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
int b = z + m;
z >>= 1;
if (x >= b)
x -= b, z += m;
}
return z;
}
```
這裡的 m 就是 $d_m$,每次迭代都會乘上 $\frac{1}{4}$ 也就是右移 2,z 則是 $c_m$ 每次迭代都會乘上 $\frac{1}{2}$ 也就是右移 1
commit: [84ac45f](https://github.com/dockyu/End-of-term/blob/84ac45fca06759b85c63e3ffff1e0b84d3341012/3_1/main.c)
:::success
延伸問題:
1. 解釋上述程式碼運作原理並重新整理數學式 (原題目故意略去某些細節),並嘗試用第 2 週測驗題提到的 `ffs` / `fls` 取代 `__builtin_clz`,使程式不依賴 GNU extension,且提供分支和無分支 (branchless) 的實作。
2. 在 Linux 核心原始程式碼找出對整數進行平方根運算的程式碼,並解說其應用案例,至少該包含 block 目錄的程式碼。
3. 閱讀〈Analysis of the lookup table size for square-rooting〉和〈Timing Square Root〉,嘗試撰寫更快的整數開平分根運算程式。
:::
#### 使用 fls 取代 `__builtin_clz`
clz 等於 `31 - fls`,所以 `31 - __builtin_clz(x)` 等於 `31 - (31 - fls)` 等於 `fls`
[fls 程式碼](https://github.com/torvalds/linux/blob/master/include/asm-generic/bitops/fls.h)
commit: [40b7d29](https://github.com/dockyu/End-of-term/blob/40b7d294ab26541bb3323002d17ce4810ead3608/3_1/main.c)
#### fls 使用 branchless
我發現 fls 每次 x 要位移以及 r 要減掉的數值都一樣,所以可以用同一個變數來表示,而這個數值都是 2 的冪,所以都可以用1位移得到
```c
if (!(x & 0xffff0000u)) {
x <<= 16;
r -= 16;
}
```
:::danger
注意書寫規範:
* 無論標題和內文中,漢字和英語/數值字元之間要有空白字元 (對排版和文字搜尋有利)
:::
原本的版本為 0 才需要操作,所以如果將如果為 0 時可以得到 1 非 0 時得到 0 ,就可以用這個結果位移使的變數的值只有 0 或是目標值(要位移、要減掉的值)
```c
shift = ((x & 0xffff0000u) == 0) << 4; /* 16 */
x <<= shift;
r -= shift;
```
commit: [12f6334](https://github.com/dockyu/End-of-term/commit/12f63340c7ec4fff0125ba628bfeaeb04ac3d866)
---
### 測驗 2
[第 3 週測驗題-測驗 2](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-2)
:::danger
提供數學證明。
參見 [Linux 核心原始程式碼的整數除法](https://hackmd.io/@sysprog/linux-intdiv) 和 [從 √2 的存在談開平方根的快速運算](https://hackmd.io/@sysprog/sqrt)
:::
$x/10$ 等同於 $x \times \frac{1}{10}$ 約等於 $x \times \frac{13}{128}$ ,之所以使用 13 跟 128 是因為 除以 128 可以透過位移做到,乘以 13 也可以透過 bitwise 做到,且這個數字在 $1\le x \le 19$ 時的誤差在小數點後一位之下,當然也有其他數字更接近 $\frac{1}{10}$ 不過題目已經確定被除數對大為 19,所以 $\frac{13}{128}$ 已經足夠
:::danger
注意書寫規範:
* 無論標題和內文中,漢字和英語/數值字元之間要有空白字元 (對排版和文字搜尋有利)
:::
### 測驗 3
[第 3 週測驗題-測驗 3](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-3)
log2 就是在算二進位表示下最左邊的1代表的數值是2幾次方,例如 $7(0111)=2^2+2^1+2^0$,最大的是$2^2$所以$\log_{2}7=2$
版本一
```c
int ilog2(int i)
{
int log = -1;
while (i) {
i >>= 1;
log++;
}
return log;
}
```
這個版本不斷右移直到 i 為0,此時位移的次數是最左邊的1的位置,假設最左邊的 1 是第 x 位,代表的數值是 $2^{x-1}$,所以只要紀錄位移次數再減1就會是答案
:::danger
注意書寫規範:
* 無論標題和內文中,漢字和英語/數值字元之間要有空白字元 (對排版和文字搜尋有利)
:::
版本二
```c
static size_t ilog2(size_t i)
{
size_t result = 0;
while (i >= 2^16) {
result += 16;
i >>= 16;
}
while (i >= 2^8) {
result += 8;
i >>= 8;
}
while (i >= 2^4) {
result += 4;
i >>= 4;
}
while (i >= 2) {
result += 1;
i >>= 1;
}
return result;
}
```
這個版本用二分搜尋法每次只判斷一半的長度的最大值,如果是就代表最左邊的1在右邊那一半的左邊,所以直接加上一半的長度
因為每次一半都只需要執行一遍,所以不需要用到迴圈,可以改寫為
```c
static size_t ilog2(size_t i)
{
size_t result = 0;
if (i >= 2^16) {
result += 16;
i >>= 16;
}
if (i >= 2^8) {
result += 8;
i >>= 8;
}
if (i >= 2^4) {
result += 4;
i >>= 4;
}
while (i >= 2) {
result += 1;
i >>= 1;
}
return result;
}
```
最下面並不是一半,所以還是必須一個位元一個位元檢查,但是前三次已經節省了非常多次檢查,這裡最多只需要檢查8次
版本三
```c
int ilog32(uint32_t v)
{
return (31 - __builtin_clz(v | 1));
}
```
計算最高的 set bit 其實就是算左邊有幾個 0,可以使用 clz 函數來得到,`v|1` 是為了必面傳入 0 讓 clz 回傳未定義的結果
:::success
延伸問題:
2. 在 Linux 核心原始程式碼找出 log2 的相關程式碼 (或類似的形式),並解說應用案例
:::
```c
/**
* ilog2 - log base 2 of 32-bit or a 64-bit unsigned value
* @n: parameter
*
* constant-capable log of base 2 calculation
* - this can be used to initialise global variables from constant data, hence
* the massive ternary operator construction
*
* selects the appropriately-sized optimised version depending on sizeof(n)
*/
#define ilog2(n) \
( \
__builtin_constant_p(n) ? \
((n) < 2 ? 0 : \
63 - __builtin_clzll(n)) : \
(sizeof(n) <= 4) ? \
__ilog2_u32(n) : \
__ilog2_u64(n) \
)
```
[drivers/net/ethernet/amd/pds_core/core.c](https://elixir.bootlin.com/linux/latest/source/drivers/net/ethernet/amd/pds_core/core.c#L355)
```c
cidi.adminq_ring_size = ilog2(pdsc->adminqcq.q.num_descs);
cidi.notifyq_ring_size = ilog2(pdsc->notifyqcq.q.num_descs);
```
看這個檔名應該是和網路相關的功能會用到 ilog ,但是我看不懂他的用途
### 測驗 4
[第 3 週測驗題-測驗 4](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-4)
$S_t = \left\{
\begin{array}{l}
Y_0& t = 0 \\
\alpha Y_t + (1 - \alpha)\cdot S_{t-1}& t > 0 \\
\end{array}
\right.$
```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;
}
```
一開始我覺得程式碼和公式對不上,參考 [vax-r 的筆記後](https://hackmd.io/@vax-r/linux2024-homework4#%E6%B8%AC%E9%A9%97%E5%9B%9B) 才知道這是定點數運算,證據是 ewma_read 時應該要回傳目前的$S_t$,但是實際回傳的卻是 `avg->internal >> avg->factor` 所以 `avg-factor` 應該就是定點數的 scaling factor 了
如果把 $2^{avg->weight}$ 暫時替換成 $w$,不為 0時可以看作
$\cfrac{S_{t-1}*w-S_{t-1}+Y_t}{w}=\cfrac{(w-1)S_{t-1}+Y_t}{w}=(1-\cfrac{1}{w})S_{t-1}+\cfrac{Y_t}{w}$ ,可已看出 $\cfrac{1}{w}$ 是 $\alpha$
:::success
延伸問題:
2. 在 Linux 核心原始程式碼找出 EWMA 的相關程式碼 (或類似的形式),並解說應用案例,至少該涵蓋無線網路裝置驅動程式。
:::
[drivers/net/wireless/ralink/rt2x00/rt2x00link.c](https://elixir.bootlin.com/linux/latest/source/drivers/net/wireless/ralink/rt2x00/rt2x00link.c#L25)
`rt2x00link_get_avg_rssi` 函式傳入一個 `struct ewma_rssi *ewma` , 他的功能是取得無線網路的 RSSI,並且函式中有一行 `avg = ewma_rssi_read(ewma);` 看起來在取得平均值,但是奇怪的是我在 linux 的程式碼裡找不到這個函數的定義
### 測驗 5
[第 3 週測驗題-測驗 5](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-5)
因為要對 $\log_{2}x$ 取 ceil,可以先觀察 $x$ 持續增加時什麼時候 $\lceil \log_{2}x \rceil$ 會增加,當 x 剛好是 2 的冪時$\lceil \log_{2}x \rceil$會是這個結果的最大的 x ,例如當 $x=8=2^{3}$ 時 $\lceil \log_{2}x \rceil=3$,當 $x=9=2^{3}+1$ 時 $\lceil \log_{2}x \rceil=4$ 觀察 8 (1000) 和 9 (1001) 的二進位,在 `x--;` 時8的對高位 set bit 會變成 3 而 9 不會變(還是 4),在 `shift = (x > 0x3) << 1;` 後此時 r 和 shift 都是被位移過後的數值所以第 1 個位元一定是0,`x >>= shift;` 因為 shift 一定是 2 如果 $x$ 只有 3 位元那位移後只有 1 位元最大就是 1,判斷 x 是否大於1就可以用來當作要不要+1的標準
:::success
延伸問題:
2. 改進程式碼,使其得以處理 x = 0 的狀況,並仍是 branchless
:::
$\log_{2}0$ 應該是無解的所以我希望 可以輸出 -1
:::danger
為何選 `-1`?在工程上有無額外風險?
> 因為 0 以及正數是 $\log$ 的可能輸出值,所以無解的問題我認為應該不能是 0 或正數,那就只能輸出負數了
:::
[Branchless code that maps zero, negative, and positive to 0, 1, 2](https://stackoverflow.com/questions/1610836/branchless-code-that-maps-zero-negative-and-positive-to-0-1-2) 有一則回應
```c
int c = (n > 0) - (n < 0);
```
可以把正數對應到 1,0 對應到 0,-1 對應到 -1
在這一題中只會有正數以及 0,分別討論這兩種情況
||正數|0|
|-|-|-|
|c|0x00000001|0x00000000|
|c-1|0x00000000|0xFFFFFFFF|
如果把最後的結果和 c-1 做 or 運算,正數不會變,而0則會變成-1
```c
int ceil_ilog2(int x)
{
uint32_t r, shift;
uint32_t u = x;
uint32_t z = (x > 0) - (x < 0);
z--;
u--;
r = (u > 0xFFFF) << 4;
u >>= r;
shift = (u > 0xFF) << 3;
u >>= shift;
r |= shift;
shift = (u > 0xF) << 2;
u >>= shift;
r |= shift;
shift = (u > 0x3) << 1;
u >>= shift;
return ((r | shift | u > 1) + 1) | z;
}
```
上面的程式碼在 $x=1$ 時會回傳 1,但是正確答案是 0,所以 $x=1$ 時要減 1,根據程式碼也可以理解為 $x=1$ 時不要做 +1 這件事,如果將 `+1` 改為 `+(x!=1)` 就只會在 $x=1$ 時變成 `+0`
```c
int ceil_ilog2(int x)
{
uint32_t r, shift;
uint32_t u = x;
uint32_t z = (x > 0) - (x < 0);
z--;
u--;
r = (u > 0xFFFF) << 4;
u >>= r;
shift = (u > 0xFF) << 3;
u >>= shift;
r |= shift;
shift = (u > 0xF) << 2;
u >>= shift;
r |= shift;
shift = (u > 0x3) << 1;
u >>= shift;
return ((r | shift | u > 1) + (x != 1)) | z;
}
```
## TODO: [第 4 週測驗](https://hackmd.io/@sysprog/linux2024-quiz4)
> 包含延伸問題
### 測驗 1
[第 4 週測驗題-測驗 1](https://hackmd.io/@sysprog/linux2024-quiz4#%E6%B8%AC%E9%A9%97-1)
$popcount(x) = x - \left \lfloor{{\dfrac{x}{2}}}\right \rfloor - \left \lfloor{{\dfrac{x}{4}}}\right \rfloor - ... - \left \lfloor{{\dfrac{x}{2^{{31}}}}}\right \rfloor$
這個方法可以只看其中一個 set bit 來解釋,如果x的第5個 bit $b_5$是1,那在$\left \lfloor{{\dfrac{x}{2}}}\right \rfloor$ (x>>1) 的 $b_4$ 會是1,在$\left \lfloor{{\dfrac{x}{4}}}\right \rfloor$ (x>>2) 的 $b_3$ 會是1,在$\left \lfloor{{\dfrac{x}{8}}}\right \rfloor$ (x>>3) 的 $b_2$ 會是1,在$\left \lfloor{{\dfrac{x}{16}}}\right \rfloor$ (x>>4) 的 $b_1$ 會是1
因為 $b_4+b_3+b_2+b_1=b_5-1$ 所以 $b_5-b_4-b_3-b_2-b_1=1$由此可知這個方法每個 set bit 會帶給結果 1,所以有幾個 set bit 結果就會是這個值
如果對每個 nibble 這麼做最後每個 nibble 中的 set bit 數量就會存在這個 nibble 裡,下面的程式碼就是在將每一個 nibble 相加
```c
v = (v + (v >> 4)) & 0x0F0F0F0F;
v *= 0x01010101;
return v >> 24;
```
```c
int totalHammingDistance(int* nums, int numsSize)
{
int total = 0;;
for (int i = 0;i < numsSize;i++)
for (int j = 0; j < numsSize;j++)
total += __builtin_popcount(nums[i] ^ nums[j]);
return total >> 1;
}
```
因為兩個迴圈會讓一個 pair 被算到兩次,所以結果要除以2,等同右移 1 bit,這個方法在 [LeetCode 477. Total Hamming Distance](https://leetcode.com/problems/total-hamming-distance/submissions/1295714268/) Time Limit Exceeded
可以將第二個迴圈的起始設為 i+1 就不會有重複計算的問題,所以也不用除以2,可以通過 LeetCode 的測試
```c
int totalHammingDistance(int* nums, int numsSize) {
int total = 0;;
for (int i = 0;i < numsSize;i++)
for (int j = i+1; j < numsSize;j++)
total += __builtin_popcount(nums[i] ^ nums[j]);
return total;
}
```
:::success
2. 不依賴 popcount,嘗試以上述歸納的規則,針對 LeetCode 477. Total Hamming Distance 撰寫出更高效的程式碼。
> 所有 Number 中從第 0 個 bit 開始計算全部該 bit 的 Hamming Distance ,在全部累加起來就是 Total Hamming Distance
:::
```c
int totalHammingDistance(int* nums, int numsSize) {
int total = 0;;
int bits;
for (int i = 0; i < 32; i++) {
bits = 0;
for (int j = 0; j < numsSize; j++) {
bits += nums[j] >> i & 0x1;
}
total += bits * (numsSize-bits);
}
return total;
}
```
[LeetCode submit](https://leetcode.com/problems/total-hamming-distance/submissions/1296382756/)
我本來想直接寫32個判斷式,後來同學提醒我可以位移後和1做 and 運算,就可以判斷出特定位元是不是1
### 測驗 2
[第 4 週測驗題-測驗 2](https://hackmd.io/@sysprog/linux2024-quiz4#%E6%B8%AC%E9%A9%97-2)
##### 版本一
經知道 $x \mod 3$ 其實就是奇數與偶數 bit 的1的數量的差
並且可以推導出下面的公式
\begin{align}
popcount(x \And \overline{m}) - popcount(x \And m) = popcount(x \bigoplus m) - popcount(m)
\end{align}
```c=
int mod3(unsigned n)
{
n = popcount(n ^ 0xAAAAAAAA) + 23;
n = popcount(n ^ 0x2A) - 3;
return n + ((n >> 31) & 3);
}
```
第3行其實就等於`奇偶的差+39`,因為只有 32 位元,所以奇偶的差的範圍是 `-16 ~ 16` ,第一行的 n 的範圍是 `23 ~ 55`,所以要再做一次但是因為最大才55(只會用到最前面6個位元)所以可以用 0x2A 當作 m ,做完第二行後 n 的範圍是 `-3 ~ 3`(其實不會有3這個結果,因為只有在 $n=21$ 時才會是 3,但是這不在第一行的結果範圍裡),如果是負的後面 `((n >> 31) & 3)` 會是3,把 n 補回正的(如下表)
|$n$|$n_{16}$|`((n >> 31) & 3)`|`n + ((n >> 31) & 3)`|
|-|-|-|-|
|-3|0xFFFFFFFD|3|0|
|-2|0xFFFFFFFE|3|1|
|-1|0xFFFFFFFF|3|2|
|0|0x00000000|0|0|
|1|0x00000001|0|1|
|2|0x00000002|0|2|
##### 版本二
```c=
int mod3(unsigned n)
{
static char table[33] = {2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1 };
n = popcount(n ^ 0xAAAAAAAA);
return table[n];
}
```
第3行明顯跟上面推導的公式不一樣少了 `-16` ,所以 n 的範圍是 `0 ~ 32` ,而 $n=0$ 實際上 $\mod 3$ 的結果是-16,對照下表,因為範圍是 `0 ~ 32` 而且我們已經知道每一個數的答案,那麼直接用一個 lookup table 比較快
|$n$|$\mod 3$|正確答案|
|:-:|:-:|:-:|
|0|-16|2|
|1|-15|0|
|2|-14|1|
|3|-13|2|
|4|-12|0|
|$\vdots$|$\vdots$|$\vdots$|
|31|15|0|
|32|16|1|
:::success
延伸問題:
1. 參照《[Hacker's Delight](https://en.wikipedia.org/wiki/Hacker%27s_Delight)》和 CS:APP 第二章,解釋上述程式碼運作原理。
:::
[Hackers Delight](https://doc.lagout.org/security/Hackers%20Delight.pdf) 給出兩個版本
##### 版本一
```c
int remu7(unsigned n) {
static char table [75] = {0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4};
n = (n >> 15) + (n & 0x7FFF);
n = (n >> 9) + (n & 0x001FF);
n = (n >> 6) + (n & 0x0003F);
return table[n];
}
```
\begin{split}
n=d*k+r
\end{split}
\begin{split}
n &\equiv d*k+r \mod 7 \\
&\equiv d*k \mod 7 + r \mod 7
\end{split}
\begin{split}
n &\equiv d*8^{k}+r \mod 7 \\
&\equiv (d*8^k \mod 7) + (r \mod 7) \\
&= (d \mod 7) + (r \mod 7)
\end{split}
$(n>>15) = \cfrac{n}{2^{15}} = \cfrac{n}{8^5} = d$ 所以經過三次操作之後的 $n' \equiv n \mod 7$ 只是 $n'$ 的範圍被縮小到 `0x4A` 以下,已經可以直接用 table 列出所有解
##### 版本二
```c
int remu7(unsigned n) {
n = (0x24924924*n + (n >> 1) + (n >> 4)) >> 29;
return n & ((int)(n - 7) >> 31);
}
```
\begin{split}
n &\equiv x (\mod 7) \\
n &= d*7 + x \\
\cfrac{8n}{7} &= \cfrac{8*(d*7+x)}{7} = \cfrac{56d+8x}{7} = 8d+\frac{8}{7}x \\
\cfrac{8n}{7} (\mod 8) &= 8d+\frac{8}{7}x (\mod 8) = \frac{8}{7}x (\mod 8)
\end{split}
$x$ 只會是 0~6
|$x$|$\frac{8}{7}x$|$\lfloor\frac{8}{7}x\rfloor$|
|-|-|-|
|0|0|0|
|1|1.1|1|
|2|2.2|2|
|3|3.4|3|
|4|4.5|4|
|5|5.7|5|
|6|6.8|6|
由上表可知在所有情況下 $\lfloor \cfrac{8}{7}n \rfloor \mod 8 = n \mod 7$,而$\lfloor \cfrac{8}{7}n \rfloor$ 可以透過 $n*(\cfrac{2^{32}}{7})/2^{29}$ 得到
\begin{split}
n*(\cfrac{2^{32}}{7})/2^{29} &= n*(613566756.571_{10})/2^{29} \\
&= n*(613566756_{10} + 0.5_{10} + 0.0625_{10} + 0.0085)/2^{29} \\
&= n*(\text{0x}24924924 + \cfrac{1}{2} + \cfrac{1}{2^4} + 0.0085)/2^{29} \\
&= (\text{0x}24924924*n + \cfrac{n}{2} + \cfrac{n}{2^4} + 0.0085*n)/2^{29} \\
&= (\text{0x}24924924*n + (n>>1) + (n>>4) + 0.0085*n) >> 29 \\
&\approx (\text{0x}24924924*n + (n>>1) + (n>>4)) >> 29
\end{split}
所以第一行 `n = (0x24924924*n + (n >> 1) + (n >> 4)) >> 29` 就是在做 $\lfloor \cfrac{8}{7}n \rfloor \mod 8$ ,因為一開始 `0x24924924*n` 等同 $n*\cfrac{2^{32}}{7}$ 如果將 $n$ 拆成 $b_{31}*2^{31}+b_{30}*2^{30}+\dots+b_3*2^3+b_2*2^2+b_1*2^1+b_0*2^0$,$b_2*2^2$以上的部分會溢出,相當於被捨棄掉了所以 $n$ 的範圍只會在 $b_2*2^2+b_1*2^1+b_0*2^0$ 可以表示的範圍裡也就是 0~7
|$n$|$(n - 7)$|`(int)(n - 7) >> 31`||
|-|-|-|-|
|0~6|負數|`0xFFFFFFFF`|因為有號數負數右移會補1|
|7|0|`0x00000000`||
所以當 $n<7$ 時和 `0xFFFFFFFF` 做 `&` 運算結果不變,當$n=7$ 時和 `0x0` 做 `&` 運算結果為 0,因為 $7\equiv0\mod 7$
:::info
理解這個程式碼之後我突然冒出一個問題,這個函數適用於所有數字嗎?因為他是靠捨棄 $0.0085n$ 當作取 floor ,但是如果 $n$ 很大,導致 $0.0085n>1$ 那不就會多捨棄了1嗎?這樣做出來的 remainder 應該會少1?
> 本手法有輸入範圍的限制。
:::
---
## TODO: [第 7 週測驗](https://hackmd.io/@sysprog/linux2024-quiz7)
> 包含延伸問題
### 測驗 1
[第 7 週測驗題-測驗 1](https://hackmd.io/@sysprog/linux2024-quiz7#%E6%B8%AC%E9%A9%97-1)
SIMD within a register (SWAR) 是軟體最佳化技巧之一,讓程式可以在特定的資料格式上平行處理
一次判斷一個 unsigned long,所以 mask 裡要放滿要找的字 d ,因為32位元和64位元的系統的 unsigned long 長度不同,所以需要在 mask 不夠長時將它複製成兩倍長
```c
for (unsigned int i = 32; i < LBLOCKSIZE * 8; i <<= 1)
mask |= mask << i;
```
每次都判斷一個 LBLOCKSIZE 有沒有 d
```c
while (len >= LBLOCKSIZE) {
if (DETECT_CHAR(*asrc, mask))
break;
```
判斷方式是,如果有 d ,那這個 unsigned char 的空間和 ~d 應該要完全顛倒,and 運算出來為0
```c
#define DETECT_NULL(X) (((X) - (0x01010101)) & ~(X) & (0x80808080))
```
接著再在這個 LBLOCKSIZE 裡找到 d 即可
DETECT_NULL 的判斷方式
```c
#define DETECT_NULL(X) \
(((X) - (0x0101010101010101)) & ~(X) & (0x8080808080808080))
```
拆開成 8 bit 來看,$y = 8 bit$ ,如果 $y$ 非0 那麼有兩種可能
1. $y<0x80$ 最左邊的 bit 為0,此時 $y-1$ 最左邊的 bit 為0
2. $y\ge 0x80$ 最左邊的 bit 為1,此時 $~y$ 最左邊的 bit 為0
這兩種可能和 0x80808080 做 and 一定是 0,所以非0的結果一定是0
:::success
延伸問題:
2. 比較 Linux 核心原本 (與平台無關) 的實作和 memchr_opt,設計實驗來觀察隨著字串長度和特定 pattern 變化的效能影響
:::
這是 [linux memchr](https://github.com/torvalds/linux/blob/master/lib/string.c#L787) 的程式碼,一次只判斷一個 char
當可以找到 d 時,執行時間只和 d 的位置有關,當找不到 d 時執行時間應該和字串長度有關,所以我的實驗會分成兩個情況討論
###### 可以找到 d 的情況
固定 string 的長度為 100005 ,紀錄第一次出現的 d 的位置增長會花費的時間的變化
![with_d](https://hackmd.io/_uploads/S1mKTkv8C.png)
###### 找不到 d 的情況
觀察不同長度的字串在不包含 d 的情況下的執行時間
![without_d](https://hackmd.io/_uploads/HJGoT1D8C.png)
:::success
延伸問題:
研讀 [2022 年學生報告](https://hackmd.io/@arthur-chang/linux2022-quiz8),在 Linux 核心原始程式碼找出 x86_64 或 arm64 對應的最佳化實作,跟上述程式碼比較,並嘗試舉出其中的策略和分析
> 下一步:貢獻程式碼到 Linux 核心
> [Phoronix 報導](https://www.phoronix.com/news/Linux-Kernel-Faster-memchr)
:::
[Phoronix 報導](https://www.phoronix.com/news/Linux-Kernel-Faster-memchr)
> We use word-wide comparing so that 8 characters can be compared at the same time on CPU.
可以同時比對8個 char ,速對會快差不多4倍,不過我做的實驗沒有快到這麼多
[[PATCH 0/2] x86: Optimize memchr() for x86-64](https://lore.kernel.org/lkml/20220528081236.3020-1-arthurchang09@gmail.com/T/#Z2e.:..:20220528081236.3020-2-arthurchang09::40gmail.com:1arch:x86:include:asm:string_64.h) 裡有對於 memchr 的最佳化的程式碼,可以直接跳到 [arch/x86/lib/string_64.c](https://lore.kernel.org/lkml/20220528081236.3020-1-arthurchang09@gmail.com/T/#Z2e.:..:20220528081236.3020-2-arthurchang09::40gmail.com:1arch:x86:lib:string_64.c) ,(但是我在現在的 linux 裡已經找不到這個檔案了只有 [arch/x86/lib/string_32.c](https://github.com/torvalds/linux/blob/f2661062f16b2de5d7b6a5c42a9a5c96326b8454/arch/x86/lib/string_32.c))
```diff
+// SPDX-License-Identifier: GPL-2.0
+#include <linux/string.h>
+#include <linux/export.h>
+#include <linux/align.h>
+
+/* How many bytes are loaded each iteration of the word copy loop */
+#define LBLOCKSIZE (sizeof(long))
+
+#ifdef __HAVE_ARCH_MEMCHR
+
+void *memchr(const void *cs, int c, size_t length)
+{
+ const unsigned char *src = (const unsigned char *)cs, d = c;
+
+ while (!IS_ALIGNED((long)src, sizeof(long))) {
+ if (!length--)
+ return NULL;
+ if (*src == d)
+ return (void *)src;
+ src++;
+ }
+ if (length >= LBLOCKSIZE) {
+ unsigned long mask = d << 8 | d;
+ unsigned int i = 32;
+ long xor, data;
+ const long consta = 0xFEFEFEFEFEFEFEFF,
+ constb = 0x8080808080808080;
+
+ /*
+ * Create a 8-bytes mask for word-wise comparing.
+ * For example, a mask for 'a' is 0x6161616161616161.
+ */
+
+ mask |= mask << 16;
+ for (i = 32; i < LBLOCKSIZE * 8; i <<= 1)
+ mask |= mask << i;
+ /*
+ * We perform word-wise comparing with following operation:
+ * 1. Perform xor on the long word @src and @mask
+ * and put into @xor.
+ * 2. Add @xor with @consta.
+ * 3. ~@xor & @constb.
+ * 4. Perform & with the result of step 2 and 3.
+ *
+ * Step 1 creates a byte which is 0 in the long word if
+ * there is at least one target byte in it.
+ *
+ * Step 2 to Step 4 find if there is a byte with 0 in
+ * the long word.
+ */
+ asm volatile("1:\n\t"
+ "movq (%0),%1\n\t"
+ "xorq %6,%1\n\t"
+ "lea (%1,%4), %2\n\t"
+ "notq %1\n\t"
+ "andq %5,%1\n\t"
+ "testq %1,%2\n\t"
+ "jne 2f\n\t"
+ "add $8,%0\n\t"
+ "sub $8,%3\n\t"
+ "cmp $7,%3\n\t"
+ "ja 1b\n\t"
+ "2:\n\t"
+ : "=D"(src), "=r"(xor), "=r"(data), "=r"(length)
+ : "r"(consta), "r"(constb), "r"(mask), "0"(src),
+ "1"(xor), "2"(data), "3"(length)
+ : "memory", "cc");
+ }
+
+ while (length--) {
+ if (*src == d)
+ return (void *)src;
+ src++;
+ }
+ return NULL;
+}
+EXPORT_SYMBOL(memchr);
+#endif
```
1. Perform xor on the long word @src and @mask and put into @xor.
2. Add @xor with @consta.
3. ~@xor & @constb.
4. Perform & with the result of step 2 and 3.
第一步跟原本的差不多,如果有 NULL 就會有一個 8 byte 是 0x00,接著的步驟是找出有沒有 0x00
### 針對其他學員 (含授課教師和社會人士) 在開發紀錄頁面提出的問題或建議
[Linux 核心專題: 重做 lab0](https://hackmd.io/@sysprog/BkltCg5IA)
[Linux 核心專題: 井字遊戲改進](https://hackmd.io/@sysprog/rkF2M6FIA)
[Linux 核心專題: 浮點數運算案例探討](https://hackmd.io/@sysprog/BJDY-YoL0)
[Linux 核心專題: 高效網頁伺服器](https://hackmd.io/@sysprog/Byq9C2q80)
[Linux 核心專題: 重作第四次作業](https://hackmd.io/@sysprog/r1WImnwIC)