---
tags: Linux, linux2022
---
# 2022q1 Homework2 (quiz2)
contributed by < `NOVBobLee` >
> [作業要求](https://hackmd.io/@sysprog/BybKYf5xc)
## 測驗 1
對兩個無號 32 位元整數取平均值。
```c
/* method 1 */
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
return (a + b) / 2;
}
```
方法一會有 overflow 的問題，我們將 `a`, `b` 都代入 $4294967295$ （ `0xFFFF'FFFF` ），結果會得到 $2147483647$ ，是代入值的一半，顯然跟我們期望的平均值不同（期望輸出為 $4294967295$ ），以下用 `uint8_t` 為例。
| `uint8_t` | 二進位 | 備註 |
|:---------:| -------------:| :--: |
| `a` | `1111'1111` | |
| `b` | `1111'1111` | |
| `a + b` | `(1)'1111'1111` | 溢位 |
| `(a + b) / 2` | `0111'1111` | |
無號整數的溢位會發生在相加超過上限值或者相減變為負數的時候，若我們知道誰的值比較大，那可以使用方法二，避免以上兩者的發生。
$\dfrac{l+h}{2} = \dfrac{2l-l+h}{2} = l+\dfrac{h-l}{2} \\
\text{1. }h-l \geq 0 \\
\text{2. Suppose that }\exists\ h, l \in \text{uint32_t s.t. }\alpha = l+\dfrac{h-l}{2} > \max(\text{uint32_t}) \\
\text{Then, }l+\dfrac{h-l}{2} = \dfrac{l+h}{2} > \max(\text{uint32_t}) \\
\text{But, } \dfrac{l+h}{2} \text{ has upper bound } \max(\text{uint32_t}) \text{ , happens when } h=l=\max(\text{uint32_t}) \rightarrow\!\leftarrow$
```c
/* method 2 */
#include <stdint.h>
uint32_t average(uint32_t low, uint32_t high)
{
return low + (high - low) / 2;
}
```
但需要事先知道誰大誰小，而方法三就避免掉減法，也就不需要比大小了。在整數上做除法，得到的是商數，所以還需要對兩者的餘數做處裡，從真值表可以看出該處理 EXP1 應該為 `a & b & 1` （用 `&1` 限定 $2^0$ 位置的位元）。
$\dfrac{a+b}{2} = \dfrac{a}{2} + \dfrac{b}{2}$
| `a ＆ 1` | `b & 1` | `(a & 1 + b & 1) / 2` |
| :-: | :-: | :-: |
| `0` | `0` | `0` |
| `0` | `1` | `0` |
| `1` | `0` | `0` |
| `1` | `1` | `1` |
```c
/* method 3 */
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
return (a >> 1) + (b >> 1) + (EXP1);
}
```
那有沒有可能把方法三再縮短，可以觀察 `a + b` 的位元計算，將 `a` 自己拆開成與 `b` 相同的部份跟不同的部份，同樣對 `b` 做拆開。
```c
a = a & b + (a ^ b) & a
b = a & b + (a ^ b) & b
^^^^^ ^^^^^^^^^^^
```
然後相加 `a` 和 `b` ，可以看到相同的部份就是兩倍，而不同的部份則剛好是互補。
```c
a + b = (a & b) << 1 + (a ^ b) & a + (a ^ b) & b
= (a & b) << 1 + (a ^ b)
```
相加可以看到左移一位有溢位的可能，不過我們還會再做一個除以 $2$ 的右移，就避開了溢位的危險，而方法四即使用此運算方式。
```c
(a + b) / 2 = ((a & b) << 1 + (a ^ b)) >> 1 // 示意，左移和右移要先相消
= (a & b) + (a ^ b) >> 1
```
```c
/* method 4 */
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
return (EXP2) + ((EXP3) >> 1);
}
```
:::info
延伸問題:
1. 解釋上述程式碼運作的原理
2. 比較上述實作在編譯器最佳化開啟的狀況，對應的組合語言輸出，並嘗試解讀 (可研讀 [CS:APP 第 3 章](https://hackmd.io/@sysprog/CSAPP-ch3))
3. 研讀 Linux 核心原始程式碼 [include/linux/average.h](https://github.com/torvalds/linux/blob/master/include/linux/average.h)，探討其 [Exponentially weighted moving average](https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average) (EWMA) 實作，以及在 Linux 核心的應用
> 移動平均（Moving average)，又稱滾動平均值、滑動平均，在統計學中是種藉由建立整個資料集合中不同子集的一系列平均數，來分析資料點的計算方法。
> 移動平均通常與時間序列資料一起使用，以消除短期波動，突出長期趨勢或周期。短期和長期之間的[閾值](https://terms.naer.edu.tw/detail/1085244/)取決於應用，移動平均的參數將相應地設置。例如，它通常用於對財務數據進行技術分析，如股票價格、收益率或交易量。它也用於經濟學中研究國內生產總值、就業或其他宏觀經濟時間序列。
:::
以下為比較有開最佳化和沒開的編譯組合語言結果，使用 AT&T 語法。
### 方法一
沒有開啟最佳化結果：
```c
; gcc -S averages.c
; uint32_t average1(uint32_t a, uint32_t b) { return (a + b) / 2; }
average1:
.LFB0:
.cfi_startproc
endbr64
pushq %rbp ; push frame pointer into stack
.cfi_def_cfa_offset 16
.cfi_offset 6, -16
movq %rsp, %rbp ; frame pointer = stack pointer
.cfi_def_cfa_register 6
movl %edi, -4(%rbp) ; write 1st argument into stack
movl %esi, -8(%rbp)
movl -4(%rbp), %edx ; read from stack, assign the value to register
movl -8(%rbp), %eax
addl %edx, %eax ; %eax += %edx , i.e. (a + b)
shrl %eax ; %eax >>= 1 , right shift
popq %rbp ; pop frame pointer from stack
.cfi_def_cfa 7, 8
ret ; return (%eax is result register)
.cfi_endproc
```
:::info
1. `endbr64` 是什麼？
CET (Intel's Control-Flow Enforcement Technology) ， GCC 8 開始加入支援，用來防止 [ROP](https://en.wikipedia.org/wiki/Return-oriented_programming), [JOP/COP](https://security.stackexchange.com/questions/201196/concept-of-jump-oriented-programming-jop) 攻擊。
來源 [What does the endbr64 instruction actually do?](https://stackoverflow.com/a/69226244/16257547)
3. `%rbp`, `%rsp` 是什麼？
> [Intel® 64 and IA-32 Architectures Software Developer’s Manual](https://www.intel.com/content/dam/www/public/us/en/documents/manuals/64-ia-32-architectures-software-developer-instruction-set-reference-manual-325383.pdf)
> If the nesting level is 0, the processor pushes the frame pointer from the BP/EBP/RBP register onto the stack, copies the current stack pointer from the SP/ESP/RSP register into the BP/EBP/RBP register, and loads the SP/ESP/RSP register with the current stack-pointer value minus the value in the size operand. For nesting levels of 1 or greater, the processor pushes additional frame pointers on the stack before adjusting the stack pointer. These additional frame pointers provide the called procedure with access points to other nested frames on the stack.
> Releases the stack frame set up by an earlier ENTER instruction. The LEAVE instruction copies the frame pointer (in the EBP register) into the stack pointer register (ESP), which releases the stack space allocated to the stack frame. The old frame pointer (the frame pointer for the calling procedure that was saved by the ENTER instruction) is then popped from the stack into the EBP register, restoring the calling procedure’s stack frame.
4. `.cfi_*` 是什麼？
CFI (call frame information) directives 是屬於 GNU AS extension ， assembler directives 之一，可以管理 call frames ，在有些架構裡需要使用這些 CFI directives 來處理錯誤。
來源一 [What are CFI directives in Gnu Assembler (GAS) used for?](https://stackoverflow.com/a/2529237/16257547)
來源二 [CFI directives - Using as](https://sourceware.org/binutils/docs-2.17/as/CFI-directives.html#CFI-directives)
:::
開啟最佳化結果（之後都將 `.cfi_*`, `endbr64` 去除），其結果 `-O0`, `-O1`, `-O2` 皆相同：
```c
; gcc -S averages.c -O0
; uint32_t average1(uint32_t a, uint32_t b) { return (a + b) / 2; }
average1:
leal (%rdi,%rsi), %eax ; sum of 1st and 2nd argument, i.e. %eax = a + b
shrl %eax ; right shift
ret
```
從兩者可以看出指令數的差別，還有就是開了最佳化後，連 stack 都不見了XD 沒開最佳化的時候，他會先把引數存進 stack 裡，再從 stack 讀回到 register 裡面，不知道是不是因為 argument register 不能拿來使用才這樣做？（可以使用，見方法三使用最佳化編譯結果）
這邊有點疑問，在 [GCC Optimize Options](https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html) 是說預設就是 `-O0` ，但 `-O0` 結果跟沒加 flag 的結果不同。
### 方法二
沒開最佳化，可以看到第 7 到 10 行，剛開始使用 `%eax` 計算，最後又換到 `%edx` 上，其實可以直接使用 `%edx` 即可。
```c=
; uint32_t average2(uint32_t low, uint32_t high) { return low + (high - low) / 2; }
average2:
pushq %rbp
movq %rsp, %rbp
movl %edi, -4(%rbp) ; write 1st argument into stack
movl %esi, -8(%rbp)
movl -8(%rbp), %eax ; read from stack, assign the value to register
subl -4(%rbp), %eax ; %eax -= -4(%rbp), i.e. %eax = high - low
shrl %eax ; %eax >>= 1
movl %eax, %edx
movl -4(%rbp), %eax
addl %edx, %eax ; %eax += %edx, i.e. %eax += low
popq %rbp
ret
```
開最佳化之後（ `-O0`, `-O1`, `-O2` 三者結果相同）， stack 也不見了。
```c
; uint32_t average2(uint32_t low, uint32_t high) { return low + (high - low) / 2; }
average2:
movl %esi, %eax ; copy 2nd argument to %eax
subl %edi, %eax ; %eax -= %edi, i.e. %eax = high - low
shrl %eax ; %eax >>= 1
addl %edi, %eax ; %eax += %edi, i.e. %eax += low
ret
```
### 方法三
沒開最佳化一樣會使用 stack ，跟開最佳化的結果比較後， stack 還是多餘的動作。
```c=
; uint32_t average3(uint32_t a, uint32_t b) { return (a >> 1) + (b >> 1) + (a & b & 1); }
average3:
pushq %rbp
movq %rsp, %rbp
movl %edi, -4(%rbp) ; write 1st argument into stack
movl %esi, -8(%rbp)
movl -4(%rbp), %eax ; %eax = a
shrl %eax ; %eax >>= 1
movl %eax, %edx
movl -8(%rbp), %eax ; %eax = b
shrl %eax ; %eax >>= 1
addl %eax, %edx ; %edx = (a >> 1) + (b >> 1)
movl -4(%rbp), %eax ; %eax = a
andl -8(%rbp), %eax ; %eax &= b
andl $1, %eax ; %eax &= 1
addl %edx, %eax ; %eax += %edx , i.e. (a >> 1) + (b >> 1) + (a & b & 1)
popq %rbp
ret
```
開最佳化後（ `-O0`, `-O1`, `-O2` 結果一樣），沒有使用 stack ，而 register 的使用上比沒開最佳化的更靈活，像是開最佳化後的第 3 到 7 行對比沒開的第 7 至 12 行，還有開最佳化的第 8 到 10 行對比沒開的第 13 到 16 行。
```c=
; uint32_t average3(uint32_t a, uint32_t b) { return (a >> 1) + (b >> 1) + (a & b & 1); }
average3:
movl %edi, %eax ; %eax = a
shrl %eax ; %eax >>= 1
movl %esi, %edx ; %edx = b
shrl %edx ; %edx >>= 1
addl %edx, %eax ; %eax += %edx , i.e. (a >> 1) + (b >> 1)
andl %esi, %edi ; b &= a ,use argument register
andl $1, %edi ; b &= 1
addl %edi, %eax ; add together
ret
```
### 方法四
沒開最佳化有使用 stack ，在 register 的選擇上會優先使用 `%eax` ，所以出現第 9 行的交換。
```c=
; uint32_t average4(uint32_t a, uint32_t b) { return (a & b) + ((a ^ b) >> 1); }
average4:
pushq %rbp
movq %rsp, %rbp
movl %edi, -4(%rbp) ; write argument register into stack
movl %esi, -8(%rbp)
movl -4(%rbp), %eax ; read from stack, %eax = a
andl -8(%rbp), %eax ; %eax &= b , i.e. %eax = a & b
movl %eax, %edx
movl -4(%rbp), %eax ; %eax = a
xorl -8(%rbp), %eax ; %eax ^= b , i.e. %eax = a ^ b
shrl %eax ; %eax >>= 1
addl %edx, %eax ; %eax = (a & b) + ((a ^ b) >> 1)
popq %rbp
ret
```
這次開啟最佳化 `-O2` 有些變化，與 `-O0`, `-O1` 不一樣，兩個都沒有使用 stack 。
```c
; -O0, -O1
; uint32_t average4(uint32_t a, uint32_t b) { return (a & b) + ((a ^ b) >> 1); }
average4:
movl %edi, %eax ; %eax = a
xorl %esi, %eax ; %eax ^= b
shrl %eax ; %eax >>= 1
andl %esi, %edi ; a &= b
addl %edi, %eax ; add together
ret
```
雖然做的事跟 `-O0`, `-O1` 一樣，不過 `-O2` 在 register 上的使用上更靈活，第 4 行將 argument register `%edi` 的值複製進 `%eax` 之後就拿來存 `a & b` 的結果了，換句話說，運用 register 的手段比之前的多。
```c=
; -O2
; uint32_t average4(uint32_t a, uint32_t b) { return (a & b) + ((a ^ b) >> 1); }
average4:
movl %edi, %eax ; %eax = a
andl %esi, %edi ; a &= b
xorl %esi, %eax ; %eax ^= b
shrl %eax ; %eax >>= 1
addl %edi, %eax ; add together
ret
```
看完之後做比較，方法一二的使用指令數比較少，不過有 overflow 的可能或是需要 a, b 的大小關係。而方法三四以比較 general 的方法，且方法四的使用指令數比較少。
（TODO: 測試）
### Linux EWMA 實作與應用
EWMA ([Exponentially Weighted Moving Average](https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average)) 是一種特殊的取平均方式，在一資料集中不同的子集合各取其平均，從名稱中可以看到移動，意指不同的子集，而其平均方式會使用指數權重。給一資料集合序列 $Y$ ，而其 EWMA $S$ 計算方式如下：
$S_t = \begin{cases}
Y_0 & \text{if } t = 0 \\
\alpha Y_t + (1 - \alpha)\cdot S_{t-1} & \text{if } t > 0
\end{cases}$
其計算非常簡單，最初的 EWMA $S_0$ 為資料本身 $Y_0$ ，因為在之前還沒有計算 EWMA 。而之後每增加一筆新的資料 $Y_t$ 會乘上一個權重 $\alpha$ ，再加上之前計算的 EWMA $S_{t-1}$ 乘上 $1-\alpha$ ，亦即以前的資料也會影響到現在的 EWMA ，可參考下面的圖例（ $\alpha=0.2$ ）。
```vega
{
"data": {
"values": [
{"xl": 0, "yl": 200, "data": "1. Y(0)=200"},
{"xl": 1, "yl": 0, "data": "1. Y(0)=200"},
{"xl": 2, "yl": 0, "data": "1. Y(0)=200"},
{"xl": 3, "yl": 0, "data": "1. Y(0)=200"}
]
},
"mark": "bar",
"encoding": {
"x": {"field": "xl", "type": "ordinal", "title": "t"},
"y": {"aggregate": "sum", "field": "yl", "type": "quantitative", "title": "EWMA S(t)"},
"color": {"field": "data", "type": "nominal", "title": "Data Y(t)"}
}
}
```
```vega
{
"data": {
"values": [
{"xl": 0, "yl": 200, "data": "1. Y(0)=200"},
{"xl": 1, "yl": 160, "data": "1. Y(0)=200"},
{"xl": 1, "yl": 30, "data": "2. Y(1)=150"},
{"xl": 2, "yl": 0, "data": "1. Y(0)=200"},
{"xl": 3, "yl": 0, "data": "2. Y(1)=150"}
]
},
"mark": "bar",
"encoding": {
"x": {"field": "xl", "type": "ordinal", "title": "t"},
"y": {"aggregate": "sum", "field": "yl", "type": "quantitative", "title": "EWMA S(t)"},
"color": {"field": "data", "type": "nominal", "title": "Data Y(t)"}
}
}
```
```vega
{
"data": {
"values": [
{"xl": 0, "yl": 200, "data": "1. Y(0)=200"},
{"xl": 1, "yl": 160, "data": "1. Y(0)=200"},
{"xl": 1, "yl": 30, "data": "2. Y(1)=150"},
{"xl": 2, "yl": 128, "data": "1. Y(0)=200"},
{"xl": 2, "yl": 24, "data": "2. Y(1)=150"},
{"xl": 2, "yl": 36, "data": "3. Y(2)=180"},
{"xl": 3, "yl": 0, "data": "1. Y(0)=200"}
]
},
"mark": "bar",
"encoding": {
"x": {"field": "xl", "type": "ordinal", "title": "t"},
"y": {"aggregate": "sum", "field": "yl", "type": "quantitative", "title": "EWMA S(t)"},
"color": {"field": "data", "type": "nominal", "title": "Data Y(t)"}
}
}
```
```vega
{
"data": {
"values": [
{"xl": 0, "yl": 200, "data": "1. Y(0)=200"},
{"xl": 1, "yl": 160, "data": "1. Y(0)=200"},
{"xl": 1, "yl": 30, "data": "2. Y(1)=150"},
{"xl": 2, "yl": 128, "data": "1. Y(0)=200"},
{"xl": 2, "yl": 24, "data": "2. Y(1)=150"},
{"xl": 2, "yl": 36, "data": "3. Y(2)=180"},
{"xl": 3, "yl": 102, "data": "1. Y(0)=200"},
{"xl": 3, "yl": 19, "data": "2. Y(1)=150"},
{"xl": 3, "yl": 29, "data": "3. Y(2)=180"},
{"xl": 3, "yl": 20, "data": "4. Y(3)=100"}
]
},
"mark": "bar",
"encoding": {
"x": {"field": "xl", "type": "ordinal", "title": "t"},
"y": {"aggregate": "sum", "field": "yl", "type": "quantitative", "title": "EWMA S(t)"},
"color": {"field": "data", "type": "nominal", "title": "Data Y(t)"}
}
}
```
其中 $\alpha$ 又稱 smoothing factor ，取值介於 $[0, 1]$ 之間，越低則舊的資料的影響越大，相反地，越高則是新加入的資料影響越大。通常會取比較低的值，使得新加入的資料不會讓現有的變化趨勢波動過大，有 smoothing 的作用。
以下為 Linux 的 EWMA 實作，位置於 [include/linux/average.h](https://github.com/torvalds/linux/blob/master/include/linux/average.h) 。要使用 EWMA 需要先用此 macro 將 EWMA 專用的結構和函式建立起來，不同名稱的 EWMA 會有不同的對應結構名稱與函式。第一個引數 `name` 為 EWMA 名稱，用來識別不同的 EWMA 。第二個參數 `_precision` 為 EWMA 的精密度是在小數第幾位（二進位）。第三個參數 `_weight_rcp` 是權重的倒數，而且要為 $2$ 的次方，直接看數學式 $S_t = \dfrac{1}{w}Y_t + (1-\dfrac{1}{w})S_{t-1}\text{ , }t> 0$ ，使用倒數的原因是權重 $\dfrac{1}{w}$ 範圍可以被限定在 $(0, 1]$ 之間了。
```c=
/*
* Exponentially weighted moving average (EWMA)
*
* This implements a fixed-precision EWMA algorithm, with both the
* precision and fall-off coefficient determined at compile-time
* and built into the generated helper funtions.
*
* The first argument to the macro is the name that will be used
* for the struct and helper functions.
*
* The second argument, the precision, expresses how many bits are
* used for the fractional part of the fixed-precision values.
*
* The third argument, the weight reciprocal, determines how the
* new values will be weighed vs. the old state, new values will
* get weight 1/weight_rcp and old values 1-1/weight_rcp. Note
* that this parameter must be a power of two for efficiency.
*/
#define DECLARE_EWMA(name, _precision, _weight_rcp) \
struct ewma_##name { \
unsigned long internal; \
}; \
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; \
} \
```
所有建立的函式名稱和結構名稱都會接上 `name` 這個名稱，用以區分不同的 EWMA 。最開始會先宣告 `struct ewma_name` （之後用 name 代表接上的部份），裡面只有一個成員 `internal` ，紀錄我們的 EWMA 。再來是自己宣告結構變數，然後使用 `ewma_name_init` 初始化該變數。
在函式中可以看到很多 `BUILD_BUG_ON` ，那是 Linux 設計在 compile-time 檢查錯誤的 macro 。這邊會檢查 `_precision`, `_weight_rcp` 是否為常數表示式， `_precision` 是否在 $30$ 以內， `_weight_rcp` 是否為 $2$ 的冪（且不可為零），若沒有符合，在編譯的時候就會跳出錯誤。
這邊我有個問題，為何要三個函式都重複檢查，何不在其中一個函數裡作所有檢查就好？我目前猜測是跳出錯誤訊息的時候，可以明確指出哪些函式會使用到這些含有錯誤的表示式，在除錯時可以比較全盤的掌握錯誤所在。
```c=24
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); \
} \
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 ，另一個拿來加入新資料更新 EWMA 用的。 `ewma_name_add` 會加入一筆新資料，更新我們的 EWMA ，可以看到 `_weight_rcp` 用 $2$ 的冪是為了使用位元的左移右移，會經過 `ilog2` 轉換成次方數 `weight_rcp` 。
$S_0 = Y_0\cdot 2^p$
$\begin{split}
S_t &= \big((S_{t-1}\cdot 2^w - S_{t-1}) + Y_t\cdot 2^p\big)\cdot 2^{-w} \text{ , }t>0 \\
&=(1-\dfrac{1}{2^w})\cdot S_{t-1} + \dfrac{1}{2^w}\cdot Y_t\cdot 2^p \text{ , }t>0
\end{split}$
而精密度的實作是將每一筆資料加入的時候都預先左移，這樣紀錄在 `internal` 裡的 EWMA 其實是實際值的 $2^{precision}$ 倍，這樣多出來的右側位數就可以拿來計算小數的部份，也因為如此，查看 EWMA 的時候，使用 `ewma_name_read` 裡會做一個右移，將 `internal` 還原成實際的 EWMA 。
在 Linux 裡很多網路相關的 driver 有使用 EWMA ，有算晶片溫度、資料吞吐量、訊號強度 (rssi) 等平均，而其他像是網路協議和 GPU 的驅動也有在使用。從搜尋結果來看，這巨集設計也是不錯，可以從 `name` 看出他的作用為何。
```c
drivers/net/virtio_net.c
54:DECLARE_EWMA(pkt_len, 0, 64)
drivers/net/wireless/ath/ath5k/ath5k.h
1255:DECLARE_EWMA(beacon_rssi, 10, 8)
drivers/net/wireless/realtek/rtw88/main.h
642:DECLARE_EWMA(tp, 10, 2);
738:DECLARE_EWMA(rssi, 10, 16);
1497:DECLARE_EWMA(thermal, 10, 4);
1543:DECLARE_EWMA(evm, 10, 4);
1544:DECLARE_EWMA(snr, 10, 4);
drivers/net/wireless/realtek/rtw89/core.h
860:DECLARE_EWMA(tp, 10, 2);
1832:DECLARE_EWMA(rssi, 10, 16);
2384:DECLARE_EWMA(thermal, 4, 4);
drivers/net/wireless/ralink/rt2x00/rt2x00.h
258:DECLARE_EWMA(rssi, 10, 8)
drivers/net/wireless/mediatek/mt76/mt76x02_mac.h
34:DECLARE_EWMA(pktlen, 8, 8);
drivers/net/wireless/mediatek/mt76/mt76.h
232:DECLARE_EWMA(signal, 10, 8);
drivers/net/wireless/mediatek/mt7601u/mt7601u.h
134:DECLARE_EWMA(rssi, 10, 4);
drivers/net/wireless/intel/iwlwifi/mvm/mvm.h
560:DECLARE_EWMA(rate, 16, 16)
drivers/net/wireless/mediatek/mt76/mt7921/mt7921.h
112:DECLARE_EWMA(rssi, 10, 8);
net/batman-adv/types.h
578:DECLARE_EWMA(throughput, 10, 8)
net/mac80211/sta_info.h
126:DECLARE_EWMA(avg_signal, 10, 8)
373:DECLARE_EWMA(mesh_fail_avg, 20, 8)
374:DECLARE_EWMA(mesh_tx_rate_avg, 8, 16)
433:DECLARE_EWMA(signal, 10, 8)
net/mac80211/ieee80211_i.h
435:DECLARE_EWMA(beacon_signal, 4, 4)
drivers/gpu/drm/drm_self_refresh_helper.c
56:DECLARE_EWMA(psr_time, 4, 4)
drivers/gpu/drm/i915/gt/intel_context_types.h
24:DECLARE_EWMA(runtime, 3, 8);
drivers/gpu/drm/i915/gt/intel_engine_types.h
125:DECLARE_EWMA(_engine_latency, 6, 4)
```
## 測驗 2
以無分支方式求出最大值。
```c
#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b)
{
return a ^ ((EXP4) & -(EXP5));
}
```
先來看看 XOR 運算的其中一些特性：
1. `(a ^ b) ^ a = (b ^ a) ^ a = b ^ (a ^ a) = b`
2. `(a ^ b) ^ b = a ^ (b ^ b) = a`
3. `a ^ 0 = a`
以上特性在 `max` 的引數裡兩者選其一可以使用，當要替換 `a` 變成 `b` ，那麼 `a` 後面可以接 `a ^ b` ，而不想替換的話則是接 `0` ，這樣看起來 `&` 後面接一個 bitmask 可以辦到。所以當 max 為 `a` 時需要 mask 為 `0x0` ，而當 max 為 `b` 時則需要 mask 為 `0xFFFF'FFFF` ，而題目有提示 EXP5 要用 relational operator ，那應該填 `a > b` 。當 `a > b` 為真，則 `(a ^ b) & (-1)` 變成 `(a ^ b) & 0xFFFF'FFFF` ，得到 `(a ^ b)` ，最後 `a ^ (a ^ b)` 等於 `a` 。當 `a <= b` 同理。其中有用到 integer promotions 。
先看 relational operator 結果的型態，是 `int` 。
> N1256 (6.5.8-6)
> Each of the operators < (less than), > (greater than), <= (less than or equal to), and >= (greater than or equal to) shall yield 1 if the specified relation is true and 0 if it is false. The result has type int.
再來看 integer promotions ，由於 `(a ^ b)` 結果是 `unsigned` ，而 `-1` 是 `int` ，根據 integer promotions 的規則， `-1` 會改變型態成 `unsigned` ，而其值則為 `0xFFFF'FFFF` ，沒有改變。
> N1256 (8.3.1.8-1)
> [...] Otherwise, the integer promotions are performed on both operands. Then the following rules are applied to the promoted operands:
>> If both operands have the same type, then no further conversion is needed.
>> Otherwise, if both operands have signed integer types or both have unsigned integer types, the operand with the type of lesser integer conversion rank is converted to the type of the operand with greater rank.
>> Otherwise, if the operand that has unsigned integer type has rank greater or equal to the rank of the type of the other operand, then the operand with signed integer type is converted to the type of the operand with unsigned integer type.
>> Otherwise, if the type of the operand with signed integer type can represent all of the values of the type of the operand with unsigned integer type, then the operand with unsigned integer type is converted to the type of the operand with signed integer type.
>> Otherwise, both operands are converted to the unsigned integer type corresponding to the type of the operand with signed integer type.
:::spoiler Integer Conversion Rank
> N1256 (6.3.1.1-1)
> Every integer type has an integer conversion rank defined as follows:
— No two signed integer types shall have the same rank, even if they have the same representation.
— The rank of a signed integer type shall be greater than the rank of any signed integer type with less precision.
— The rank of `long long int` shall be greater than the rank of `long int`, which shall be greater than the rank of `int`, which shall be greater than the rank of `short int`, which shall be greater than the rank of `signed char`.
— The rank of any unsigned integer type shall equal the rank of the corresponding signed integer type, if any.
— The rank of any standard integer type shall be greater than the rank of any extended integer type with the same width.
— The rank of `char` shall equal the rank of `signed char` and `unsigned char`.
— The rank of `_Bool` shall be less than the rank of all other standard integer types.
— The rank of any enumerated type shall equal the rank of the compatible integer type.
— The rank of any extended signed integer type relative to another extended signed integer type with the same precision is implementation-defined, but still subject to the other rules for determining the integer conversion rank.
— For all integer types `T1`, `T2`, and `T3`, if `T1` has greater rank than `T2` and `T2` has greater rank than `T3`, then `T1` has greater rank than `T3`.
:::
:::info
延伸問題:
1. 解釋上述程式碼運作的原理
2. 針對 32 位元無號/有號整數，撰寫同樣 [branchless](https://en.wikipedia.org/wiki/Branch_(computer_science)) 的實作
3. Linux 核心也有若干 branchless / branch-free 的實作，例如 [lib/sort.c](https://github.com/torvalds/linux/blob/master/lib/sort.c):
```c
/*
* Logically, we're doing "if (i & lsbit) i -= size;", but since the
* branch is unpredictable, it's done with a bit of clever branch-free
* code instead.
*/
__attribute_const__ __always_inline
static size_t parent(size_t i, unsigned int lsbit, size_t size)
{
i -= size;
i -= size & -(i & lsbit);
return i / 2;
}
```
請在 Linux 核心原始程式碼中，找出更多類似的實作手法。請善用 `git log` 檢索。
:::
程式碼 `uint32_t` 也可以換成 `int32_t` ，其無號有號都適用。若要算 `min` 則將 bitmask 的生成條件改變即可，即 `a < b` 。
## 測驗 3
使用 [Euclidean algorithm](https://en.wikipedia.org/wiki/Euclidean_algorithm) 或稱輾轉相除法，求出 64 位元 GCD ([greatest common divisor](https://en.wikipedia.org/wiki/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;
}
```
GCD 演算法（又稱 binary GCD alogrithm）：
1. If both $x$ and $y$ are $0$, gcd is zero $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.
以下用可位元操作的改寫。可以看出下面改寫比原算法多了 `v/=2`, `u/=2` 的動作，其對應的是 GCD 演算法的第 3 到 5 點，相較於使用比較費力的取模操作，改使用右移與減法更適當。
```c
#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v; /* GCD 演算法第 1, 2 點 */
int shift;
for (shift = 0; !((u | v) & 1); shift++) {
u /= 2, v /= 2; /* GCD 演算法第 3 點 */
}
while (!(u & 1))
u /= 2; /* GCD 演算法第 4(5) 點 */
do {
while (!(v & 1))
v /= 2; /* GCD 演算法第 4(5) 點 */
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (COND);
return RET;
}
```
在 do-while 的區塊，其本質跟原算法相同，即計算餘數，再拿餘數當除數，換句話說就是輾轉相除法，停止條件跟原算法相同，條件為整除，餘數為零的時候，所以 COND 應填為 `v` 。而回傳值應該也跟原算法回傳商數一樣，不過我們有先做 GCD 演算法第 3 點，先提出公因數 $2^{shift}$ ，所以要再乘回到我們的商數 `u` 上，即 RET 應為 `u << shift`。雖說是輾轉相除，但實作是使用反覆的減法代替除法（意義上跟原本的輾轉相除法的商數、餘數相同）。
:::info
延伸問題:
1. 解釋上述程式運作原理;
2. 在 x86_64 上透過 `__builtin_ctz` 改寫 GCD，分析對效能的提升;
3. Linux 核心中也內建 GCD (而且還不只一種實作)，例如 [lib/math/gcd.c](https://github.com/torvalds/linux/blob/master/lib/math/gcd.c)，請解釋其實作手法和探討在核心內的應用場景。
:::
使用 [GCC built-in functions](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html) 的`__builtin_ctz` 作改寫。
> Built-in Function: int __builtin_ctz (unsigned int x)
>> Returns the number of trailing 0-bits in x, starting at the least significant bit position. If x is 0, the result is undefined.
```c
#include <stdint.h>
uint64_t min(uint64_t a, uint64_t b)
{
return a ^ ((a ^ b) & -(a < b));
}
uint64_t gcd64(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v; /* GCD 演算法第 1, 2 點 */
/* GCD 演算法第 3 點 */
int shift = min(__builtin_ctz(u), __builtin_ctz(v));
u >>= shift;
v >>= shift;
/* GCD 演算法第 4(5) 點 */
u >>= __builtin_ctz(u);
do {
/* GCD 演算法第 4(5) 點 */
v >>= __builtin_ctz(v);
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u << shift;
}
```
在 Linux 內建的 GCD 實作使用大量 bitwise 操作，有分使用 `__ffs` 的版本跟沒有使用的版本，有使用 `__ffs` 的版本會比 even/odd algorithm (binary GCD algorithm) 的速度更快。
```c=
#if !defined(CONFIG_CPU_NO_EFFICIENT_FFS)
/* If __ffs is available, the even/odd algorithm benchmarks slower. */
/**
* gcd - calculate and return the greatest common divisor of 2 unsigned longs
* @a: first value
* @b: second value
*/
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` 的版本， ffs (find first set) 是在找 LSB (least significant bit) 的位數，有分[硬體支援](https://en.wikipedia.org/wiki/Find_first_set#Hardware_support)和[函式庫支援](https://en.wikipedia.org/wiki/Find_first_set#Tool_and_library_support)，這邊指的是硬體。
方法本質也是輾轉相除法，不過有些改良，第 14 行一樣是 GCD 演算法第 1, 2 點，不過這個 `r` 就有意思了， `r & -r` 會得到 LSB ，在我們的算法中，剛開始有一步是算公因數 $2^{shift}$ ，這邊使用 `r & -r` 代替掉我們的迴圈。
第 16 行是 GCD 演算法的第 4(5) 點（第 3 點因為用 `r` 來算了，所以主要動機應該是第 4(5) 點），將 `b` 裡 $2$ 的因數都提出，假如可以被 $2$ 整除（第 17 行），則回傳 `r & -r` ，這是因為 `b` 只有 $2$ 這個因數，也就不需要進入輾轉相除的階段。
再來進入迴圈準備輾轉相除，第 22 行跟第 16 行一樣，假如把 $2$ 的因數提光之後互質，則回傳 `r & -r` ，同時也是輾轉相除法停止的條件之一。下一行 `a == b` 也是停止條件之一，即商數等於餘數，回傳商數乘以 $2$ 的公因數。
```c=32
#else
/* If normalization is done by loops, the even/odd algorithm is a win. */
unsigned long gcd(unsigned long a, unsigned long b)
{
unsigned long r = a | b;
if (!a || !b)
return r;
/* Isolate lsbit of r */
r &= -r;
while (!(b & r))
b >>= 1;
if (b == r)
return r;
for (;;) {
while (!(a & r))
a >>= 1;
if (a == r)
return r;
if (a == b)
return a;
if (a < b)
swap(a, b);
a -= b;
a >>= 1;
if (a & r)
a += b;
a >>= 1;
}
}
#endif
```
在沒有使用 ffs 的版本裡，可以看到使用迴圈取代 ffs 。最初第 38, 39 行一樣是 GCD 演算法第 1, 2 點。然後將 `r` 取只剩 LSB 一個位元，意思是 `r` 是 `a`, `b` 兩者最大的 $2$ 的公因數，即 GCD 演算法第 3 點。
從註解中可以看到 normalization ，在這裡的 norm 是 [$p$-adic norm ($p$-adic absolute value)](https://en.wikipedia.org/wiki/P-adic_order#p-adic_absolute_value) ，特別指 $2$-adic norm 。簡單說 $2$-adic order 就是一整數 $n$ 所擁有 $2$ 的因數最高次方數，而 $n$ 的 $2$-adic norm 則是 $2$ 的 $2$-adic order 次方之後取倒數，然後整數的長度和兩整數之間的距離可以由 $2$-adic norm 來定義。
回到註釋， normalization 就是以 `r` 為基準，將 `a`, `b` 的 $2$-adic norm 變成跟 `r` 一樣，所以就是將 `a`, `b` 右移到三者的 LBS 位置一樣，此時的 `a`, `b`, `r` 三者的 $2$-adic norm 相同，最大的 $2$ 的公因數也相同。
:::info
Let $p$ be a prime number.
1. The $p$-adic order for $\mathbb{Z}$ is the function $\nu_p:\mathbb{Z}\rightarrow\mathbb{N}$ defined by
$\nu_p(n) =
\begin{cases}
\max\{k\in\mathbb{N}:p^k|n\} & \text{if } n\neq 0\\
\infty & \text{if } n=0
\end{cases}$
2. The $p$-adic absolute value on $\mathbb{Q}$ is the function $|\cdot|_p:\mathbb{Q}\rightarrow\mathbb{R}_{\geq 0}$ defined by
$|r|_p=p^{-\nu_p(r)}$
For example, $n = 45$ gives $\nu_3(n)=2$ and $|n|_3=3^{-\nu_3(n)}=\dfrac{1}{9}$ .
:::
（待續）
## 測驗 4
以下是 bit array (bitset) 的應用。
```graphviz
digraph {
node[shape=record;];
rankdir=LR;
expand [
label="{0|1|1|0|1|0|0|...|0|1|1|0|0|1|0|1|...|1|...}";
xlabel="expand the bitmap";
];
i [
label="i";
shape=plaintext;
];
bitmap_ent [ label="<e0>bitmap[0]|<e1>bitmap[1]|...|bitmap[bitmapsize - 1]"; ];
bitmap [ label="*bitmap"; ];
bitset1 [ label="{<e0>1|1|0|0|<e1>1|0|1|...|1}"; ];
bitset0 [ label="{<e0>0|1|1|0|1|0|0|...|0}"; ];
bitmap -> bitmap_ent:e0;
bitmap_ent:e1 -> bitset1:e0;
bitmap_ent:e0 -> bitset0:e0;
i -> bitset1:e1;
}
```
`bitmap` 是一個 64 位元無號整數陣列，而每個整數都代表著一個 `bitset` ，使用 `p + i` 來計算 bit 的 index ，是將 `bitmap` 所有 bits 都攤平，所以 `p` 會是每加一次就是加 64 ，代表一個整數有 64 位元，然後加上在 `bitset` 裡的位移 `i` ，傳到 `out` 裡面。
```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;
}
```
以下使用 GNU extension 的 [`__builtin_ctzll`](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html) 作等價改寫。先觀察程式碼，改變的地方從 `p = k * 64` 開始，變成 `while (bitset != 0)` ，然後計算 trailing zero 的數量，將之與 `k * 64` 取和放進 `out` ，應該就是對應上面程式碼的 index 計算，然後做 XOR 運算，繼續下一個迴圈。所以推測應該每次迴圈就是計算一個 index ，放進 `out` 之後就要將那個 bit 給反轉，所以每次 `r` 都可以更新到左邊一個 bit 的位置，進而到 `bitset` 變為 `0` ，算完所有的 bits ，結束迴圈。
```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;
}
```
這樣 EXP6 就要填那個準備反轉的 bit ，利用 `bitset & -bitset` 可以得到 LSB 位置，即該準備反轉的 bit 。
:::info
延伸問題:
1. 解釋上述程式運作原理，並舉出這樣的程式碼用在哪些真實案例中;
2. 設計實驗，檢驗 ctz/clz 改寫的程式碼相較原本的實作有多少改進？應考慮到不同的 [bitmap density](http://www.vki.com/2013/Support/docs/vgltools-3.html);
3. 思考進一步的改進空間;
4. 閱讀 [Data Structures in the Linux Kernel](https://0xax.gitbooks.io/linux-insides/content/DataStructures/linux-datastructures-3.html) 並舉出 Linux 核心使用 bitmap 的案例;
:::
新算法可以跳過 `0`-bit ，只對 `1` 作計算。根據不同的 bit density ，密度越高對新算法越不利，因每次迴圈的計算量較多。
（檢驗待測）
## 測驗 5
LeetCode - [Fraction to Recurring Decimal](https://leetcode.com/problems/fraction-to-recurring-decimal/) ，給一分數的分子與分母，回傳其以字串表示的小數。若有循環小數，則將循環部份用括號括起來。
```c=
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "list.h"
struct rem_node {
int key;
int index;
struct list_head link;
};
static int find(struct list_head *heads, int size, int key)
{
struct rem_node *node;
int hash = key % size;
list_for_each_entry (node, &heads[hash], link) {
if (key == node->key)
return node->index;
}
return -1;
}
```
這裡對循環小數的對策是使用 hash table ，將遇到的餘數儲存，每算個新的位數，就與之前儲存的餘數作比對，檢查是否出現計算重複的餘數，即出現循環小數。
剛開始先檢查有沒有 trivial solution ，是否有分子為 $0$ 或分母為 $0$ 的情況。再來做正負號檢測，儲存在 `sign` 中。在第 35 行做了無謂的檢查，因為 `n`, `p` 在第 24, 26 行調整為正數，所以 `division` 一定是正數，不需再做檢查，即該行應改為 `sprintf(p, "%ld", (long)division);` 較佳。
```c=23
char *fractionToDecimal(int numerator, int denominator)
{
int size = 1024;
char *result = malloc(size);
char *p = result;
if (denominator == 0) {
result[0] = '\0';
return result;
}
if (numerator == 0) {
result[0] = '0';
result[1] = '\0';
return result;
}
/* using long long type make sure there has no integer overflow */
long long n = numerator;
long long d = denominator;
/* deal with negtive cases */
if (n < 0)
n = -n;
if (d < 0)
d = -d;
bool sign = (float) numerator / denominator >= 0;
if (!sign)
*p++ = '-';
long long remainder = n % d;
long long division = n / d;
sprintf(p, "%ld", division > 0 ? (long) division : (long) -division);
if (remainder == 0)
return result;
p = result + strlen(result);
*p++ = '.';
/* Using a map to record all of reminders and their position.
* if the reminder appeared before, which means the repeated loop begin,
*/
char *decimal = malloc(size);
memset(decimal, 0, size);
char *q = decimal;
size = 1333;
struct list_head *heads = malloc(size * sizeof(*heads));
for (int i = 0; i < size; i++)
INIT_LIST_HEAD(&heads[i]);
for (int i = 0; remainder; i++) {
int pos = find(heads, size, remainder);
if (pos >= 0) {
while (PPP > 0)
*p++ = *decimal++;
*p++ = '(';
while (*decimal != '\0')
*p++ = *decimal++;
*p++ = ')';
*p = '\0';
return result;
}
struct rem_node *node = malloc(sizeof(*node));
node->key = remainder;
node->index = i;
MMM(&node->link, EEE);
*q++ = (remainder * 10) / d + '0';
remainder = (remainder * 10) % d;
}
strcpy(p, decimal);
return result;
}
```
看完前半部份再來到第 54 行的迴圈，第一次迴圈 hash table 還是空的，所以 `pos` 會是 `-1` ，再來將儲存餘數進 hash table ，這邊注意到 `key` 是 `int` ，而 `remainder` 是 `long long` ，所以存進去的餘數有可能存不完整， `find` 在比對 `key` (truncated `remainder`) 與 `node->key` (truncated `remainder`) 是否相同時會有問題，應將 `key` 改為 `long long` 。之後將位數 `i` 存進 `index` ，插入 hash table ，所以 MMM 應該為 `list_add` ，而位置要跟 `find` 裡的 hash function 對應，所以 EEE 應該為 `&heads[remainder % size]` 。最後將該位數的商算出，放進 `decimal` ，然後更新 `remainder` ，進入下一個迴圈。
回到剛跳過的第 78 行區塊，在若 `find` 發現該餘數已經出現過了，則回傳前一次出現的位數 `pos` ，這時我們知道有循環小數，準備將剛才放在 `decimal` 的小數外圍用括號括起來，放到我們要回傳的 `result` 。第 80 行將 `decimal` 前面部份複製到 `result` 上，然後複製完後給一個括號，所以推測停止條件裡的 PPP 為 `pos--` 。之後將循環部份複製，最後補個下括號和 null terminator 就可以回傳了。
若沒有遇到循環小數，也就是在整除的時候，結束迴圈，將 `decimal` 複製到 `result` 回傳結束。
:::info
延伸問題:
1. 解釋上述程式碼運作原理，指出其中不足，並予以改進
> 例如，判斷負號只要寫作 `bool isNegative = numerator < 0 ^ denominator < 0`;
> 搭配研讀 [The simple math behind decimal-binary conversion algorithms](https://indepth.dev/posts/1019/the-simple-math-behind-decimal-binary-conversion-algorithms)
2. 在 Linux 核心原始程式碼的 `mm/` 目錄 (即記憶體管理) 中，找到特定整數比值 (即 fraction) 和 bitwise 操作的程式碼案例，予以探討和解說其應用場景
:::
## 測驗 6
由於要從記憶體上快速找到某資料位置，一種加速方式是對資料做[對齊](https://en.wikipedia.org/wiki/Data_structure_alignment)，讓現代的 CPU 讀寫時更有效率。而使用 [`__alignof__`](https://gcc.gnu.org/onlinedocs/gcc/Alignment.html) 可以得知某資料的對齊大小，進而讓我們對資料位置的控制更方便。
以下 `t` 用 `int` 代入為例，中間的 anonymous structure 第一個成員為 `c` 是 `char` 型態，第二個成員為 `_h` 是 `int` 型態，有 4 bytes 。由對齊的特性， `char` 的地址會在 `1` 的倍數上，換句話說，任意地址皆可，反而 `int` 是 4 bytes ，所以地址會在 `4` 的倍數上。將這個 anonymous structure 放在 `0` 上， `0` 為任意數的倍數，用第一個成員 `c` 佔據，然後以 `char` 這個最小單位給 `_h` 做位移，強制我們要計算對齊的對象 `_h` 的地址在 `1` 以後（含）。
這時有趣的事發生了，假如 `t` 的對齊是 `1` 的倍數，則此時 `_h` 就會剛好靠在 `char` 的旁邊，也就是地址 `1` 上，但現在我們 `t` 的對齊為 `4` 的倍數，所以 `_h` 不會在地址 `1` 上，而是出現在下一個 `4` 的倍數地址上，即地址 `4` ，而中間空的部份稱為 padding 。
```graphviz
digraph {
node[shape=record;]
{
rank=same;
struct [
label="struct \{ char c; int _h; \}";
shape=plaintext;
];
a [ label="<c>char|<p0>____|<p1>____|<p2>____|<i>int 0|int 1|int 2|int 3"; ];
}
addr0 [ label="address = 0"; shape=plaintext; ];
p0 [ label="padding"; shape=plaintext; ];
addr2 [ label="address = 4"; shape=plaintext; ];
addr0 -> a:c;
p0 -> a:p0;
p0 -> a:p1;
p0 -> a:p2;
addr2 -> a:i;
}
```
根據以上的例子，對齊可以由 `_h` 與 `0` 的位移量得知，所以 M 應該指 `_h` 這個成員，得到前者 `_h` 的地址，而後者 X 對應 anonymous structure 的起始位置 `0` 當基準，計算位移量。
```c
/*
* ALIGNOF - get the alignment of a type
* @t: the type to test
*
* This returns a safe alignment for the given type.
*/
#define ALIGNOF(t) \
((char *)(&((struct { char c; t _h; } *)0)->M) - (char *)X)
```
這邊有兩個問題，為何要轉換成 `(char *)` ？還有為何要相減，直接拿 `_h` 的地址不就好了。第一個問題，因為要做指標的減法，所以需要兩個指標的型態相同。而對於前者的對齊可能很多種，不如直接選 `(char *)` ，因為任何數都是 `1` 的倍數，這樣可以避免某些錯誤和編譯器的警告。
第二個問題目前只是猜測，根據 N1256 (6.5.6-9) ，兩者相減得到的型態不會是 `(char *)` ，而是 `ptrdiff_t` ，也就是會有轉換型態發生。用數學舉例，在二維空間中，兩點的座標相減得到是向量，而向量跟座標是不一樣的概念。
> N1256 (6.5.6-7)
> For the purposes of these operators, a pointer to an object that is not an element of an array behaves the same as a pointer to the first element of an array of length one with the type of the object as its element type.
> N1256 (6.5.6-9)
> When two pointers are subtracted, both shall point to elements of the same array object, or one past the last element of the array object; the result is the difference of the subscripts of the two array elements. The size of the result is implementation-defined, and its type (a signed integer type) is `ptrdiff_t` defined in the `<stddef.h>` header. [...]
> N1256 (6.2.5-27)
> A pointer to void shall have the same representation and alignment requirements as a pointer to a character type. __Similarly, pointers to qualified or unqualified versions of compatible types shall have the same representation and alignment requirements__. All pointers to structure types shall have the same representation and alignment requirements as each other . All pointers to union types shall have the same representation and alignment requirements as each other. Pointers to other types need not have the same representation or alignment requirements.
:::info
延伸問題:
1. 解釋上述程式碼運作原理
2. 在 Linux 核心原始程式碼中找出 [`__alignof__`](https://gcc.gnu.org/onlinedocs/gcc/Alignment.html) 的使用案例 2 則，並針對其場景進行解說
3. 在 Linux 核心源使程式碼找出 `ALIGN`, `ALIGN_DOWN`, `ALIGN_UP` 等巨集，探討其實作機制和用途，並舉例探討 (可和上述第二點的案例重複)。思索能否對 Linux 核心提交貢獻，儘量共用相同的巨集
:::
## 測驗 7
[FizzBuzz](https://en.wikipedia.org/wiki/Fizz_buzz) 是一個數學遊戲，剛開始由 $1$ 開始往上數，每當數字可以被 $3$ 整除，這時玩家不該喊數字，而是代替那數字喊 Fizz 。同樣地，數字可以被 $5$ 整除，則是喊 Buzz 。當可以被 $3$, $5$ 同時整除，則是喊 FizzBuzz 兩個同時。
```c
/* naive.c */
#include <stdio.h>
int main() {
for (unsigned int i = 1; i < 100; i++) {
if (i % 3 == 0) printf("Fizz");
if (i % 5 == 0) printf("Buzz");
//if (i % 15 == 0) printf("FizzBuzz"); /* no need to print it again */
if ((i % 3) && (i % 5)) printf("%u", i);
printf("\n");
}
return 0;
}
```
上面為直覺的作法，當數字 `i` 可以被 $3$ 整除就印出 Fizz ，當數字可以被 $5$ 整除則印出 Buzz ，當可以同時被兩者整除時，則在兩個 if 檢查的時候就會依序印出。最後都不能整除的時候，則是印出該數字 `i` 。
```c
string literal: "FizzBuzz%u"
offset: 0 4 8
```
觀察上面同時整除的時候，依序印出就像是兩個字接在一起，換個方式想，那可以控制字串的起始位置和長度，來選擇印出的字串為何，就像上面的 `"FizzBuzz%u"` 那樣，選擇 offset 和字串長度來控制要印出的格式。
```c
/* bitwise.c */
static inline bool is_divisible(uint32_t n, uint64_t M)
{
return n * M <= M - 1;
}
static uint64_t M3 = UINT64_C(0xFFFFFFFFFFFFFFFF) / 3 + 1;
static uint64_t M5 = UINT64_C(0xFFFFFFFFFFFFFFFF) / 5 + 1;
int main(int argc, char **argv)
{
for (size_t i = 1; i <= 100; i++) {
uint8_t div3 = is_divisible(i, M3);
uint8_t div5 = is_divisible(i, M5);
unsigned int length = (2 << KK1) << KK2;
char fmt[9];
strncpy(fmt, &"FizzBuzz%u"[(9 >> div5) >> (KK3)], length);
fmt[length] = '\0';
printf(fmt, i);
printf("\n");
}
return 0;
}
```
先跳過 `is_divisible` 函式，到後面再來分析，從名字可以看出就是判斷 `n` 是否能被 bitmask `M` 對應的整數整除。那開始看主函式，使用我們剛才的想法，需要知道要印出的字串起始位置和長度，有的資訊則是可否被 $3$, $5$ 整除，下面用表格表示相對關係。
| `dev3` | `dev5` | start index | `length` |
|:------:|:------:| :----------:|:--------:|
| `0` | `0` | `8` | `2` |
| `0` | `1` | `4` | `4` |
| `1` | `0` | `0` | `4` |
| `1` | `1` | `0` | `8` |
從此表可以看出 `length` 是每多一個可整除的因數（ $3$ or $5$ ）就多增加一倍，所以可以直接 KK1 填 `dev3` ， KK2 填 `dev5` （可以互換）。
而起始位置原為 `&"FizzBuzz%u"[(9 >> div5) >> (KK3)]` 裡的 `9` 應改為 `8` 才對，不然不符合 KK3 的填空要求。從表格觀察，若 `div5` 為 `0` ，則起始位置為 `8` ，沒有右移，若 `div5` 為 `1` ，則起始位置為 `4` ，這兩個事件在 `(8 >> div5)` 就已經寫好了。那麼再來是 `dev3` ，若為 `0` 則保持原狀，若為 `1` 則歸零，所以 KK3 應該為右移位數為零或者右移位數大於 $4$ 的一個表示式（可以將 `8` 歸零），所以可以填 `dev3 << 2` 。
:::info
延伸問題:
1. 解釋上述程式運作原理並評估 `naive.c` 和 `bitwise.c` 效能落差
* 避免 stream I/O 帶來的影響，可將 `printf` 更換為 `sprintf`
2. 分析 [Faster remainders when the divisor is a constant: beating compilers and libdivide](https://lemire.me/blog/2019/02/08/faster-remainders-when-the-divisor-is-a-constant-beating-compilers-and-libdivide/) 一文的想法 (可參照同一篇網誌下方的評論)，並設計另一種 bitmask，如「可被 3 整除則末位設為 1」「可被 5 整除則倒數第二位設定為 1」，然後再改寫 `bitwise.c` 程式碼，試圖運用更少的指令來實作出 branchless;
* 參照 [fastmod: A header file for fast 32-bit division remainders on 64-bit hardware](https://github.com/lemire/fastmod)
3. 研讀 [The Fastest FizzBuzz Implementation](https://tech.marksblogg.com/fastest-fizz-buzz.html) 及其 [Hacker News](https://news.ycombinator.com/item?id=29413656) 討論，提出 throughput (吞吐量) 更高的 Fizzbuzz 實作
4. 解析 Linux 核心原始程式碼 [kernel/time/timekeeping.c](https://github.com/torvalds/linux/blob/master/kernel/time/timekeeping.c) 裡頭涉及到除法運算的機制，探討其快速除法的實作 (注意: 你可能要對照研讀 `kernel/time/` 目錄的標頭檔和程式碼)
> 過程中，你可能會發現可貢獻到 Linux 核心的空間，請充分討論
:::
Sign in with Wallet
Connect another wallet