# 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)