# 2022q1 Homework2 (quiz2)

contributed by < [`tinyynoob`](https://github.com/tinyynoob) >

> [作業要求](https://hackmd.io/@sysprog/BybKYf5xc)

## 測驗 1

### 解釋程式碼

這題的目標是對兩數取平均，由於直覺的寫法會有 overflow 的問題，因此這裡試圖用 bitwise operation 來達成。

#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
    return (a >> 1) + (b >> 1) + (EXP1);
}

這裡的 `>> 1` 顯然是除以 2 的意思，那麼這個 `EXP1` 應該就是來解決兩個奇數分開運算造成的錯誤，例如：

(1 >> 1) + (1 >> 1) = 0 neq 1

要確認兩個都是奇數，我想到的最簡易方法是

`EXP1` = `a & b & 1`

---

取平均的另外一種等價的寫法

uint32_t average(uint32_t a, uint32_t b)
{
    return (EXP2) + ((EXP3) >> 1);
}

根據[你所不知道的 C 語言：數值系統篇](https://hackmd.io/@sysprog/c-numerics#%E9%81%8B%E7%94%A8-bit-wise-operator)裡的提點，我們了解到可以直接在 C 語言上應用數位邏輯設計，嘗試再次推導：

(x + y) / 2
=> (x + y) >> 1
=> ((x ^ y) + cin) >> 1
=> ((x ^ y) + ((x & y) << 1)) >> 1
=> ((x ^ y) >> 1) + (x & y)

> [C Operator Precedence](https://en.cppreference.com/w/c/language/operator_precedence)

注意到不要把這個數學推導過程的左右移想成 C 語言內的左右移，否則會誤以為最後一個步驟不相等。

:::warning
這裡的加法長得是如此簡單又沒有 ripple 的問題，使我納悶在邏輯設計時為何 carry lookahead adder 要用那麼複雜的方式，或許 shift register 影響此法的硬體可行性？
:::

:::info
已解決。
前面推導的過程將 `+` 誤改成 `^` 導致我想錯了
:::

結論得

- `EXP2` = `a & b`
- `EXP3` = `a ^ b`

### 編譯器最佳化下的比較

開啟編譯器最佳化 `-O2` ，對於上述兩個函式得到以下結果：

...
average1:
	endbr64
	movl	%edi, %eax
	movl	%esi, %edx
	andl	%esi, %edi
	shrl	%eax
	shrl	%edx
	andl	$1, %edi
	addl	%edx, %eax
	addl	%edi, %eax
	ret
...
average2:
	endbr64
	movl	%edi, %eax
	andl	%esi, %edi
	xorl	%esi, %eax
	shrl	%eax
	addl	%edi, %eax
	ret
...

兩種 C 語言寫法在 assembly 層級會相差 3 行指令，主要的差異來自於第 5, 8, 10 行。

第一種算法 average1 意圖在 `%eax` 和 `%edx` 分別算出 `a >> 1` 和 `b >> 1` ，並在 `%edi` 算出 `a & b & 1` ，最後加總到返回值 `%eax` 。
可以想見若我們先做 `a & b & 1` 則可以不必用到 4 個 register ，也有機會減少 1 行指令，變成：

movl	%edi, %eax
andl    %esi, %eax
andl    $1, %eax
shrl    %edi
shrl    %esi
addl    %edi, %eax
addl    %esi, %eax
ret

不過我認為編譯器原先的作法可以降低前後指令的相依關係，帶來更好的平行性。假如以指令的行數來表示，則可以有：

$$\begin{split}
4\to 7&\to 10\\
5\to 8&\to 11\\
6\to 9&\underbrace{\to}_{同步}
\end{split}
$$

在多處理器的情況下有機會做 3 步就得出結果

第二種算法 average2 雖然生成的組合語言指令較少，但前後的相依性大，在多處理器的情況下可能需要執行更久：

$$\begin{split}
16\to 18\to 19&\to 20\\
17&\underbrace{\to}_{同步}
\end{split}
$$

不過如此短的計算應該還是會在單一處理器上完成，結論而言我認為仍是法二較佳。

:::warning
不確定我這裡對 multicore 相關的認知是否正確
:::

### 研讀 kernel 中的 EWMA

> [linux/average.h](https://github.com/torvalds/linux/blob/master/include/linux/average.h)

從定義

#define DECLARE_EWMA(name, _precision, _weight_rcp)	

可以看出這個 macro 接收了 3 個參數，並且可以從註解中得知

- `precision` 是定點數在小數點以下的精度
- `weight_rcp` 是加權平均的權重，令它為 $w$ 則
$$
y_n = \frac{1}{w}\cdot x_n +\frac{w-1}{w}\cdot y_{n-1} \tag{1}
$$
其中 $x_n$ 為原始輸入， $y_n$ 是移動加權後的值

$w$ 被規定須選為 power of 2 ，如此一來除法便可以化為簡單的右移運算。

我們主要需要解釋的是這一段程式碼：

unsigned long weight_rcp = ilog2(_weight_rcp);	
...
WRITE_ONCE(e->internal, internal ?			\
	(((internal << weight_rcp) - internal) +	\
		(val << precision)) >> weight_rcp :	\
	(val << precision));				\

令 weight_rcp = $\log_2$_weight_rcp

若 `internal` 非 0 ，則依據 Equation $(1)$ 更新 `internal` 值。若為 0 ，則可以直接省略加號右邊的部份。

`WRITE_ONCE()` 的功能是在多執行緒的情境下進行正確寫入，把第二個參數的值 assign 給第一個參數。

> [Linux内核中的READ_ONCE和WRITE_ONCE宏](https://zhuanlan.zhihu.com/p/344256943)
:::warning
不知道哪裡可以找到 `WRITE_ONCE()` 的第一手資料，它不在 **compiler.h** 裡
:::

:::warning
我推測這裡是使用定點數運算，不過 `>> precision` 不是等同於把低位元丟棄嗎，這是如何把 unsigned long 變成定點數的呢？
:::

moving average 的功能主要是讓資料變得更平滑，不因突發的波動而隨之起舞，其實它就是離散時間的低通濾波器。以這題而言，我們可以嘗試由 Z transform 計算 frequency response ：

$$\begin{aligned}
\mathcal{Z}\{y_n\} &= \mathcal{Z}\left\{\frac{1}{w}\cdot x_n +\frac{w-1}{w}\cdot y_{n-1}\right\} \\
Y(z) &= \frac{1}{w}X(z) + \frac{w-1}{w}\cdot z^{-1}Y(z) \\
H(z) = \frac{Y(z)}{X(z)} &= \frac{z}{1-w+zw}
\end{aligned}
$$

換成頻域得

$$
H(e^{j\Omega}) = \frac{e^{j\Omega}}{1-w+w\,e^{j\Omega}}
= \frac{1}{1-(1-\frac{1}{w})\,e^{-j\Omega}}
$$

對於頻率響應形式為 $\dfrac{1}{1-a\,e^{-j\Omega}}$ 的離散濾波器，有：

$$
\text{filter type }
\begin{cases}
\text{highpass}, &-1<a<0 \\
\text{lowpass}, &0<a<1 \\
\text{unstable}, &|a|\geq 1
\end{cases}
$$

因為 $0<(1-\frac{1}{w})<1$ ，我們的 moving average 是一個低通濾波器。

### 探討 kernel/sched/fair.c 中的 Low-pass filter

:::info
studying...
:::

## 測驗 3

### 解釋程式碼

#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v;
    int shift;
    for (shift = 0; !((u | v) & 1); shift++) {
        u /= 2, v /= 2;
    }
    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;
}

該程式計算 gcd ，首先處理 `u`, `v` 有 0 的情形。

第 6, 7 行基於 $\mathrm{gcd}(x,y) = 2\cdot\mathrm{gcd}(\frac{x}{2},\frac{y}{2})$ ，先將答案中 2 的成份 (以質因數分解的角度) 之 power 存到 `shift` 。

此時必有其一為奇數，接著要利用 $\mathrm{gcd}(x,y) = \mathrm{gcd}(\frac{x}{2},y)$ 的性質，先嘗試將 `u` 的大小盡量降低。

最後 11 到 21 行進行 [Euclidean algorithm](https://en.wikipedia.org/wiki/Euclidean_algorithm) ，不停做相減直到出現 1 ：
1. 由於 `v` 是每次相減後新出現的值，有可能又是偶數， 12, 13 行先確認能不能再用性質把值降下來。
2. 將大減小放到 `v` 並將原先較小的數放到 `u` 。

因此 `COND` 為 Euclidean algorithm 繼續的條件，即是 `v` 不等於 `1` ，考量在計算過程中 `u`, `v` 可能變相同，例如： $\mathrm{gcd}(3,3)=3$ ，因此 `v` 的中止條件也可能為 `0` ，總結來說 `COND` 為 `(v != 1) || (v != 0)` 。

答案便是留下來的 `u` 值再補回早先存的 `shift` ，因此 `RET` 為 `u << shift` 。

拿去測試之後發現有錯，練習使用 gdb 除蟲：
1. 以 `u = 3`, `v = 5` 呼叫
2. break 在 while (上方程式的 21 行)，依序得
3.
|`u`|`v`|
|:-:|:-:|
|3|2|
|1|2|
|1|0|
4. 接著程式卻沒有正確結束

判斷 while 條件有問題，將 `COND` 更正為 `(v != 1) && (v != 0)` 後得到理想結果。

### 透過 `__builtin_ctz` 改寫並比較效能

由於目前程式碼中存在大量判斷偶數與右移的操作，因此可以用 `ctz()` 來處理，從 [gcc doc](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html) 找到說明：

> 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.

improved 版實作程式碼：

#define min(x, y)          \
    ({                     \
        typeof(x) _x = x;  \
        typeof(y) _y = y;  \
        _x < _y ? _x : _y; \
    })

static uint64_t gcd64_imp(uint64_t u, uint64_t v)
{
    if (!u || !v)
        return u | v;
    int shift = min(__builtin_ctz(u), __builtin_ctz(v));

    u >>= __builtin_ctz(u);
    do {
        v >>= __builtin_ctz(v);
        if (u < v) {
            v -= u;
        } else {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while ((v != 1) && (v != 0));
    return u << shift;
}

為了進行實驗，我需要一些隨機的整數輸入 `gcd()` ， C 常用的 `rand()` 最大只到 32767 。考量兩種計算 gcd 的方法可能在大數 cases 有更顯著的差距，因此我參考了 [stack overflow](https://stackoverflow.com/questions/28115724/getting-big-random-numbers-in-c-c) ，使用 C++ 來產生亂數，[程式碼](https://github.com/tinyynoob/2022linux-quiz2/blob/master/RNG.cpp)。

我打算分成兩組輸入數值範圍做比較：

- 0 ~ 65535
- 65536 ~ 4294967295

:::warning
這個測試方法好像怪怪的，因為隨機出來的 pair 大部分都會互質，但是我不確定什麼樣的測試方式更好
:::

參考[共筆](https://hackmd.io/@KYWeng/rkGdultSU#%E9%87%8F%E6%B8%AC%E6%99%82%E9%96%93%E7%9A%84%E6%96%B9%E5%BC%8F)在 user space 量測時間的方式，撰寫用來測量的[程式碼](https://github.com/tinyynoob/2022linux-quiz2/blob/master/3.c)。

輸入

$ sudo taskset -c 0 ./3

進行量測

#### 結果

每組 100000 個數據進行測量，得到結果：

|| small numbers | big numbers |
|:-:|:-:|:-:|
|original|27944393|44117441 (ns)|
|improved|16084115|23787207 (ns)|
| the ratio | 1.73739 | 1.85467 |

small number 組：


big number 組：


顯示出利用 `__builtin_ctz()` 優化後的版本比原版加快將近一倍，看來 shift operation 的運算量在原版程式中佔了很大的比例；而 big number 組的差距並不比 small number 組來得顯著，與我原先猜想的不同。

### 解釋 kernel 中的實作

> [linux/lib/math/gcd.c](https://github.com/torvalds/linux/blob/master/lib/math/gcd.c)

## 測驗 4

### 解釋程式碼

首先看 naive 版

#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;
}

因為看不出這到底在做什麼，直接丟數字進去 run ：

    uint64_t bitmap[2] = {0xFF22, 0x033F};
    uint32_t out[16];
    printf("%d\n", naive(bitmap, 2, out));

得到

(gdb) run
Starting program: /home/tinyynoob/code/linux2022-quiz2/4 
18

Breakpoint 1, main () at 4.c:25
25          return 0;
(gdb) p out
$1 = {1, 5, 8, 9, 10, 11, 12, 13, 14, 15, 64, 65, 66, 67, 68, 69}

與輸入值 `{0000 0011 0011 1111, 1111 1111 0010 0010}` (左邊為高位) 一起觀察，發現 `out` 算出了每一個 set bit 的 $\log_2()$ 。而 gdb 第 3 行的 18 則是算出 set bit 的總數。理解程式功能後往下繼續看 improved 版：

#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;
}

由於要使用 `ctz()` 函式，因此每次計算完低位後就要把它設為 0 ，讓更高位可以繼續使用 `ctz()` ，以 `00001010` 舉例來說，首先 `out[0]` 會被設為 1 ，接著為了讓 `out[1]` 被設為 3 ，必須在算完 `out[0]` 之後將 `bitset` 改為 `00001000` ，因此 `EXP6` 應為 `1u << __builtin_ctzll(bitset)` ，如果這裡能把第 10 行往前顯然更好。

### 實驗比較

[實驗程式碼](https://github.com/tinyynoob/2022linux-quiz2/blob/master/4.c) (配合 [Makefile](https://github.com/tinyynoob/2022linux-quiz2/blob/master/Makefile) 會輕鬆許多)

我用了四種 bitmap density 進行實驗， set bit 的個數分別是 8, 24, 40, 56 ，每一組都跑 100000 個數據。得到如下結果：

#### 8

#### 24

#### 40

#### 56

暫且用 [Hamming weight](https://en.wikipedia.org/wiki/Hamming_weight) 來稱呼 bitmap 中 set bit 的個數。

可以看到在 weight 8 的情境下， improved 比 naive 快了將近 10 倍，不過當 Hamming weight 變大時，兩條線就越來越靠近，此結果佐證了作業說明中的：

> 若 bitmap 越鬆散 (即 1 越少)，於是 improved 的效益就更高。

可以見到當 weight = 56 時， improved 雖然仍比 naive 更好，但已經沒有顯著的改善了。

`improved()` 的效能隨著 Hamming weight 增大而隨之掉落符合我們的預期，值得注意的是在另一方面，我們可以發現 `naive()` 在 weight 56 的效能好於 weight 40 ，甚至好於 weight 8 ，目前不曉得該如何解釋這個現象

:::info
待
:::

### 進一步改進

## 測驗 7

### 探討 faster remainder

> https://lemire.me/blog/2019/02/08/faster-remainders-when-the-divisor-is-a-constant-beating-compilers-and-libdivide/

Let $n, d$ be constants.
Suppose that we want to do $n \div d$. we may consider

$$
n\div d = n \times \frac{2^N}{d} \times \frac{1}{2^N}
$$

If we can precompute $2^N/d$, then the division can be replaced by multiplication and bit-shift.

由於整數運算到小數點後就捨去，所以我嘗試用一些亂數跑跑看 accuracy 和 precision，整體看起來 precision 並不高，甚至在少部份的情況會出現低 accuracy 的結果：

> [accuracy and precision](https://en.wikipedia.org/wiki/Accuracy_and_precision)


---

接著我們嘗試分析 `is_divisible()` 函式

uint64_t c = 1 + UINT64_C(0xffffffffffffffff) / d;

// given precomputed c, checks whether n % d == 0
bool is_divisible(uint32_t n) {
  return n * c <= c - 1; 
}

簡記 `UINT64_MAX` 為 MAX

RHS = c - 1
    = 1 + MAX / d - 1
    = MAX / d

LHS = n * c
    = n * (1 + MAX / d)
    = n + MAX / d * n
    = n + MAX / d * (n / d * d + n % d)
    = n + (MAX * (n / d)) + (MAX / d) * (n % d)
    = n - (n / d) + (MAX / d) * (n % d)

因此，若 $n, d$ 滿足
- $n - n/d\leq \mathrm{MAX}/d$

則有：

$$
\begin{cases}
LHS\leq RHS &\text{if } d \text{ divides } n\\
LHS > RHS &\text{if } \text{not}
\end{cases}
$$

假設 $n, d$ 的型態皆為 `uint32_t`，該條件可被滿足。

### 補完程式碼

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"[(8 >> div5) >> (KK3)], length);
        fmt[length] = '\0';

        printf(fmt, i);
        printf("\n");
    }
    return 0;
}

直接列出全部四種情況：

- `(div3 && div5)`
- start = 0
- length = 8
- `(div3)`
- start = 0
- length = 4
- `(div5)`
- start = 4
- length = 4
- Otherwise
- start = 8
- length = 2

於是可以輕易推得：`KK1` = `div3`, `KK2` = `div5`, `KK3` = `4 * div3`