---
tags: linux2022
---
# 2022q1 Homework2 (quiz2)
contributed by < `cwl0429` >
## 測驗 `1`
觀察下列程式碼
- `a >> 1` 和 `b >> 1` 的功用都是將其除 2，此時會遇到最低位元遺失問題
- 需要將最低位元補上，方法是 a 和 b 作 `AND`，接著和 1 作 `AND`，便可取得遺失的位元，因此 `EXP1` 為 `a & b & 1`
```c
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
return (a >> 1) + (b >> 1) + (EXP1);
}
```
再觀察改寫後的程式碼，發現完全想不到該如何下手，於是嘗試用湊數字的方式，最後得出:
- `EXP2` 為 `a & b`
- `EXP3` 為 `a ^ b`
```c
uint32_t average(uint32_t a, uint32_t b)
{
return (EXP2) + ((EXP3) >> 1);
}
```
後來，參照了 [laneser](https://hackmd.io/@laneser/linux2022-quiz2) 的筆記，得知其實 `x + y = (x ^ y) + ((x & y) << 1)`，而 `(x + y) / 2` 就等於 `((x ^ y) + ((x & y) << 1)) >> 1` 也就是 `(x & y) + ((x ^ y) >> 1)`
因此得出:
- `EXP2` 為 `a & b`
- `EXP3` 為 `a ^ b`
### 比較上述在編譯器最佳化開啟的狀況
先放一張來自 [x64](https://cs.brown.edu/courses/cs033/docs/guides/x64_cheatsheet.pdf) 文件的圖片，方便檢視暫存器的大小
接著嘗試理解下列四種 average 的組合語言，使用 `gcc -S 1-1.c -O3 -fno-asynchronous-unwind-tables` 生成
- `-O3` 開啟編譯器最佳化第四層
- `-fno-asynchronous-unwind-tables` 過濾不必要的 `.cfi`
```c
uint32_t average1(uint32_t a, uint32_t b)
{
return (a + b) / 2;
}
average1:
endbr64
leal (%rdi,%rsi), %eax
shrl %eax
ret
```
根據 [Introduction to the leal instruction](https://www.cs.swarthmore.edu/~newhall/cs31/resources/addrleal.php) 的描述
> Load effective address: leal S,D # D<--&S, where D must be a register, and S is a Memory operand.
>
>leal looks like a mov instr, but does not access Memory. Instead, it takes advantage of the addressing circuitry and uses it to do arithmetic (as opposed to generating multiple arithmetic instructions to do arithmetic).
>
> (ex) if edx holds the value of x:
leal (%eax),%ecx # R[%ecx]<--&(M[R[%eax]])
>
> \# this moves the value stored in %eax to %ecx
>
> The key is that the address of (M[ at address x ]) is x, so this is moving the value stored in %eax to %ecx; there is no memory access in this instruction's execution.
- `leal` 會將 `(%rdi,%rsi)` 的值放至 `%eax`，其間沒有 memory access 發生
- 接著利用 `shrl` 將 `%eax` 右移一位元，達成除 2 的效果
```c
uint32_t average2(uint32_t low, uint32_t high)
{
return low + (high - low) / 2;
}
average2:
endbr64
movl %esi, %eax
subl %edi, %eax
shrl %eax
addl %edi, %eax
ret
```
- `movl` 將 `high` 存入 `%eax`
- `subl` 使 `%eax` 減 `low` 後，將 `%eax` 往右移一個 bit
- `addl` 將結果相加並存放於 `%eax`
```c
uint32_t average3(uint32_t a, uint32_t b)
{
return (a >> 1) + (b >> 1) + (a & b & 1);
}
average3:
endbr64
movl %edi, %eax
movl %esi, %edx
andl %esi, %edi
shrl %eax
shrl %edx
andl $1, %edi
addl %edx, %eax
addl %edi, %eax
ret
```
- 將 `x` 及 `y` 用 `movl` 讀至 `%eax` 及 `%edx`
- 接著做對應的位元運算 `andl` 及 `shrl`
- 最後便將三個結果 `%edx`, `%edi` 及 `%eax` 相加，存放在 `%eax`
```c
uint32_t average4(uint32_t a, uint32_t b)
{
return (a & b) + ((a ^ b) >> 1);
}
average4:
endbr64
movl %edi, %eax
andl %esi, %edi
xorl %esi, %eax
shrl %eax
addl %edi, %eax
ret
```
- 將 `a` 存入 `%eax` ，此時 `%edi` 及 `%eax` 為 a `%esi` 為 b
- 做對應的位元運算，`andl` 及 `xorl`
- 完成後將 `%eax` 右移
- 最後將 `a & b` 結果加至 `eax`
### 研讀 Linux 核心原始程式碼 include/linux/average.h，探討其 Exponentially weighted moving average (EWMA) 實作，以及在 Linux 核心的應用
TODO..
## 測驗 `2`
```c
#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b)
{
return a ^ ((EXP4) & -(EXP5));
}
```
觀察題目，若 a 較大，則需 `return a ^ 0` 才能成立，反之則需要 `return a ^ (a ^ b)` 才能回傳 b
到此根據老師的提示，便可猜出 `EXP4` 應為 `a ^ b`
剩下 `EXP5` 要決定回傳 `a ^ 0` 還是 `a ^ (a ^ b)`
- `a ^ 0`
- 代表 `a ^ ((a ^ b) & 0x0)`
- `a ^ (a ^ b)`
- 代表 `a ^ ((a ^ b) & 0xFFFFFFFF)`
所以得出 `EXP5` 為 `a < b`
### 針對 32 位元無號/有號整數，撰寫同樣 branchless 的實作
無號數版本能和有號數幾乎一樣，更改回傳型態及參數型態即可
```c
#include <stdint.h>
int32_t max(int32_t a, int32_t b)
{
return a ^ ((a ^ b) & -(a < b));
}
```
### Linux 核心也有若干 branchless / branch-free 的實作，請在 Linux 核心原始程式碼中，找出更多類似的實作手法。請善用 git log 檢索。
在 `drivers/gpu/drm/radeon/r600_blit.c`
> Replace int2float() with an optimized version. We use
> __fls() to find the most significant bit. Using that, the
> loop can be avoided. A second trick is to use the
> behaviour of the rotate instructions to expand
> the range of the unsigned int to float
> conversion to the full 32 bits in a branchless way.
>
>The routine is now exact up to 2^24. Above that, we truncate which is equivalent to rounding towards zero.
>
>Signed-off-by: Steven Fuerst <svfuerst@gmail.com>
這份 commit 使用 branchless 的技巧優化了 `int2float`
主要優化的位置在第 12 行
```c=
uint32_t int2float(uint32_t input)
{
u32 result, i, exponent, fraction;
if ((input & 0x3fff) == 0)
result = 0; /* 0 is a special case */
else {
exponent = 140; /* exponent biased by 127; */
fraction = (input & 0x3fff) << 10; /* cheat and only
handle numbers below 2^^15 */
for (i = 0; i < 14; i++) {
if (fraction & 0x800000)
break;
else {
fraction = fraction << 1; /* keep
shifting left until top bit = 1 */
exponent = exponent - 1;
}
}
result = exponent << 23 | (fraction & 0x7fffff); /* mask
off top bit; assumed 1 */
}
return result;
}
```
新的 `int2float` 利用 `__fls` 找到 `msb` 的位置，再利用 `ror` rotate 到正確位置，使我們不需判斷是否低於 2 ^ 24
```c
uint32_t int2float(uint32_t x)
{
uint32_t msb, exponent, fraction;
/* Zero is special */
if (!x) return 0;
/* Get location of the most significant bit */
msb = __fls(x);
/*
* Use a rotate instead of a shift because that works both leftwards
* and rightwards due to the mod(32) behaviour. This means we don't
* need to check to see if we are above 2^24 or not.
*/
fraction = ror32(x, (msb - I2F_FRAC_BITS) & 0x1f) & I2F_MASK;
exponent = (127 + msb) << I2F_FRAC_BITS;
return fraction + exponent;
}
```
## 測驗 `3`
__GCD 演算法__
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$ brcause everything divides 0.
3. If x and y are both even, $gcd(x, y) = 2 * gcd(x/2, 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(x/2, y)$.
5. On similar lines, if x is odd and y is even, then $gcd(x, y) = gcd(x, y/2)$. It is because 2 is not a common divisor.
```c
if (!u || !v) return u | v;
```
參考老師提供的 GCD 演算法，這段程式碼實作了 1. 及 2. ，也就是 gcd(0, 0), gcd(x, 0) 及 gcd(0, y) 的情況
```c
for (shift = 0; !((u | v) & 1); shift++) {
u /= 2, v /= 2;
}
```
這是 3. 的情況，當二數的最低位元皆為 0 時，便將二數都除 2 並且紀錄除的次數，也就是 `shift` 的數值
```c
while (!(u & 1))
u /= 2;
```
這是 4. 的情況，當 x 為偶數時，gcd(x ,y) = gcd(x / 2, y)
```c
while (!(u & 1))
u /= 2;
do {
while (!(v & 1))
v /= 2;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (COND);
return RET;
```
最後這段則是 5. 的情況及 GCD 的原始算法
這邊需要找出 `COND` 及 `RET` 為何，觀察 `COND` ，得知其作用為中止 while 的判斷式，那根據 `GCD` 的原始定義，當兩數相減至其中一數為 0 時則中止，而且得知
- `COND` 為 `v`
最後 `RET` 需要回傳二數中剩餘非 0 的數，也就是 `u`，但要注意在前面的 `shift` 代表了二數整除於 2 多少次數，所以需要將 `u` 左移 `shift` 次
- `RET` 為 `u << shift`
### 在 x86_64 上透過 __builtin_ctz 改寫 GCD，分析對效能的提升
> 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. 取自 [`GNU` 官方網站](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html)
`__buildin_ctz` 能取得 `unsigned int x` 從最低位元算起的連續 0-bits 個數，再回頭看 `gcd_64` 程式碼，得知主要想改善的部份為：
```c
while (!(u & 1))
u /= 2;
```
於是將相關操作改成 `__buildin_ctz` 後，得到以下實作：
```c
uint64_t gcd64_ctz(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
int shift;
int u_ctz =__builtin_ctz(u);
int v_ctz =__builtin_ctz(v);
shift = u_ctz > v_ctz ? v_ctz : u_ctz;
u >> u_ctz;
v >> v_ctz;
do {
v = v >> __builtin_ctz(v);
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u << shift;
}
```
接著利用 perf 進行效能測試
用亂數進行測試， N 的大小為 10000
```c
int main()
{
srand(time(NULL));
for(int i = 0; i < N; i++) {
//gcd64(rand(), rand());
gcd64_ctz(rand(), rand());
}
return 0;
}
```
使用以下命令進行測試
`perf stat --repeat 10 ./1-3`
`gdb64`
```shell
0.0055975 +- 0.0000618 seconds time elapsed ( +- 1.10% )
```
`gdb64_ctz`
```shell
0.002960 +- 0.000433 seconds time elapsed ( +- 14.63% )
```
可看出 ctz 的版本約快了 1.8910 倍
### Linux 核心中也內建 GCD (而且還不只一種實作)，例如 [lib/math/gcd.c](https://github.com/torvalds/linux/blob/master/lib/math/gcd.c)，請解釋其實作手法和探討在核心內的應用場景。
研究 `ffs` 的用途
> `ffs` — find first set bit in word
得知 `ffs` 能找到從右邊數起的第一個 1 位置，再根據[維基百科](https://en.wikipedia.org/wiki/Find_first_set)的解釋， `ffs` 和之前題目用到的 `ctz` 的關係能用數學式子表示
在 x 不等於 0 的情況下 $ctz(x) = ffs(x) - 1$
舉例來說， `x = 0x11000000` 則 `ctz(x) = 6` 而 `ffs(x) = 7`
但實際在 trace 此函式時，若是根據上述的定義，這個函式有機會進入無窮迴圈，才發現 `__ffs` 並不等於 `ffs`
於是重新翻了資料，找到 `__ffs` 實際定義的地方 [asm-generic/bitops/__ffs.h](https://docs.huihoo.com/doxygen/linux/kernel/3.7/include_2asm-generic_2bitops_2____ffs_8h_source.html)
```c
/**
* __ffs - find first bit in word.
* @word: The word to search
*
* Undefined if no bit exists, so code should check against 0 first.
*/
```
因此 `__ffs` 和 `ctz` 會回傳一樣的結果
另一個需要注意的地方是 `r & -r` 會保留從右往左數第一個為 1 的位元，其餘位元皆設定為 0
> 有趣的是，若將 `r & -r` 取 $log2$ ，其實就會是 `ctz(r)` 的結果
> 也學到教訓，找資料請盡量找第一手資料
解決上述兩點後，能觀察到這個版本的 `GCD` 和我們的版本差異為使用了 `__ffs` 進行加速，演算法相同
```c
unsigned long gcd(unsigned long a, unsigned long b)
{
unsigned long r = a | b;
if (!a || !b) /*1. If both x and y are 0, gcd is zero $gcd(0, 0) = 0$.*/
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;
}
}
```
接著來探討 lib/math/gcd.c 實際的應用情境
利用 `grep -r "#include <linux/gcd.h>"` 列出有哪些檔案使用了 `lib/math/gcd.c` ，發現在 sound, net, drivers 等等都有用到 gcd，這邊我嘗試理解 `linux/sound/soc/codecs/arizona.c` 內 gcd 的用途
觀察目錄，發現 `arizona.c` 和音訊解碼有關
Sign in with Wallet
Connect another wallet