---
tags: linux2022
---

# 2022q1 Homework2 (quiz2)

contributed by < `SmallHanley` >

> [linux2022-quiz2](https://hackmd.io/@sysprog/linux2022-quiz2)

## 測驗 `1`

:::success
解釋上述程式碼運作的原理
:::

考慮以下對二個無號整數取平均值的程式碼：

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

這個直覺的解法會有 overflow 的問題，若我們已知 a, b 數值的大小，可用下方程式避免 overflow：

#include <stdint.h>
uint32_t average(uint32_t low, uint32_t high)
{
    return low + (high - low) / 2;
}

接著我們可改寫為以下等價的實作：

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

**原理：**
如果是我們平常使用的兩實數取平均，可以寫成：

\begin{gather}
\frac{a + b}{2} = \frac{a}{2} + \frac{b}{2}
\end{gather}

而上述的程式碼考慮的是非負整數，因此我們可以使用右移運算子 `>>` 代替上式的除以 2。不過這還不是等價的程式碼，會忽略到兩數都是奇數時的進位問題，因此要加上 `a & b & 1`，使得當兩數都是奇數時會再加一。

我們再次改寫為以下等價的實作：

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

**原理：**

上述程式碼讓我們很容易聯想到加法器，我們先觀察一下半加器，以下是真值表 ({C, S} = A + B)：

| A | B | C | S |
| --- | --- | --- |:---:|
| 0 | 0 | 0 | 0 |
| 0 | 1 | 0 | 1 |
| 1 | 0 | 0 | 1 |
| 1 | 1 | 1 | 0 |

得到以下結果：

\begin{gather}
S = A \oplus B \\
C = A \times B
\end{gather}

而全加器則是除了兩個一位元的輸入訊號，還需要低一位的進位訊號 (即 {C~out~, S} = A + B + C~in~)。當我們要計算多位元的非負整數加法，可以利用多個全加器串接 (或稱級聯，cascade)。若要使用到半加器的結果，當前位數的運算結果，可以看作是當前位數的兩輸入訊號相加再加上前一位的進位訊號 (即 S + C~in~，注意，這裡的 S 為半加器的 Sum，C~in~ 為前一位的 C~out~)，也就是當前位數兩輸入訊號的 bitwise XOR 加上 前一位數兩輸入訊號的 bitwise AND。因此，兩多位數 (分別為 a 和 b) 相加的結果可以寫作：

((a & b) << 1) + (a ^ b)

兩數平均可以寫作：

(a & b) + ((a ^ b) >> 1)

:::success
比較上述實作在編譯器最佳化開啟的狀況，對應的組合語言輸出，並嘗試解讀 (可研讀 [CS:APP 第 3 章](https://hackmd.io/@sysprog/CSAPP-ch3))
:::

==暫略==，我目前沒有觀察到什麼有趣的現象。

:::success
研讀 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/)取決於應用，移動平均的參數將相應地設置。例如，它通常用於對財務數據進行技術分析，如股票價格、收益率或交易量。它也用於經濟學中研究國內生產總值、就業或其他宏觀經濟時間序列。
:::

在 [include/linux/average.h](https://github.com/torvalds/linux/blob/master/include/linux/average.h) 中使用巨集 `DECLARE_EWMA(name, _precision, _weight_rcp)` 用來生成相關的物件和函式。其中，`_precision` 是用來表示小數部份的精準度，也就是我們使用多少位元表示。而 `_weight_rcp` 則是權重。

而巨集內部包含以下結構體：

struct ewma_##name {						\
		unsigned long internal;					\
	};

以及以下函式：

* **void ewma_##name##_init(struct ewma_##name *e)***
將 `internal` 初始化為 `0`。
* **unsigned long ewma_##name##_read(struct ewma_##name *e)***
將 `internal` 右移 `_precision` 的位元數回傳。
* **void ewma_##name##_add(struct ewma_##name *e, unsigned long val)***
將 `val` 新增進來，新的 `internal` 會是舊的 `internal` 和 `val` 的線性組合 (`val` 佔 1 / _weight_rcp，舊的 `internal` 佔 1 - 1 / _weight_rcp)。
舊的 `internal` 佔比的部份，原本應寫作：
\begin{gather}
\mbox{(internal << precision) - (internal << precision >> weight_rcp)}
\end{gather}經化簡可寫成：
\begin{gather}
\mbox{(internal << weight_rcp) - internal}
\end{gather}

**應用：**

比方說 [drivers/net/wireless/realtek/rtw88/main.h](https://github.com/torvalds/linux/blob/8efd0d9c316af470377894a6a0f9ff63ce18c177/drivers/net/wireless/realtek/rtw88/main.h) 有使用到：

DECLARE_EWMA(tp, 10, 2);

struct ewma_tp tx_ewma_tp;

可以推測跟一些訊號處理有關係。

## 測驗 `2`

:::success
解釋上述程式碼運作的原理
:::

改寫〈[解讀計算機編碼](https://hackmd.io/@sysprog/binary-representation)〉一文的「不需要分支的設計」一節提供的程式碼 `min`，我們得到以下實作 (`max`)：

#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b)
{
    return a ^ ((a ^ b) & -(a < b));
}

**原理：**

我們知道 a 和 b 做 bitwise XOR，再和 a 做 bitwise XOR 會得到 b，因此我們可以運用這個原理讓 a 和 (a ^ b) 或是 0 做 bitwise XOR 而分別得到 b 和 a。至於要怎麼做到不使用分支，我們可以利用 a < b 的結果，加上負號，會產生兩種結果 (0x00...00 和 0xff...ff)，再加上 bitwise AND 可以巧妙的做到無分支判斷最大值。

:::success
針對 32 位元無號/有號整數，撰寫同樣 [branchless](https://en.wikipedia.org/wiki/Branch_(computer_science)) 的實作
:::

當我們要判斷一個 32 位元無號整數是否為 2 的冪，我們可以用以下程式碼判斷：

#include <stdint.h>
#include <stdbool.h>
bool is_power_of_2(uint32_t n)
{
    uint32_t t = 0x80000000;
    while (t) {
        if (n == t) {
            return true;
        }
        t >>= 1;
    }
    return false;
}

我們可以改寫成以下 branchless 的實作：

#include <stdint.h>
#include <stdbool.h>
bool is_power_of_2(uint32_t n)
{
    return !(n & (n - 1)) && n;
}

:::success
Linux 核心也有若干 branchless / branch-free 的實作，例如 [lib/sort.c](https://github.com/torvalds/linux/blob/master/lib/sort.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` 檢索。
:::

在 [Replace int2float() with an optimized version.](https://git.kernel.org/pub/scm/linux/kernel/git/torvalds/linux.git/commit/?id=747f49ba67b8895a5831ab539de551b916f3738c) (commit 747f49ba) 中，將 32 位元整數的 bit pattern 轉成浮點數的 bit pattern，有使用到類似的實作手法 (不過目前 r600_blit 相關程式已背棄用)。以下是 [IEEE 754-2008](https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4610935) 定義的 `binary32` 格式，其中藍、綠、紅分別代表 sign、exponent 和 fraction，分別佔用 1、8、23 位元：

digraph G                     
{ 
	n
    [
    shape = none
	label = <<table border="0" cellspacing="0"><tr>
            <td border="1" bgcolor="blue">0</td>
            <td border="1" bgcolor="green">1</td>
            <td border="1" bgcolor="green">1</td>
            <td border="1" bgcolor="green">0</td>
            <td border="1" bgcolor="green">1</td>
            <td border="1" bgcolor="green">1</td>
            <td border="1" bgcolor="green">0</td>
            <td border="1" bgcolor="green">0</td>
            <td border="1" bgcolor="green">0</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">0</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">0</td>
            <td border="1" bgcolor="red">0</td>
            <td border="1" bgcolor="red">0</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">0</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">0</td>
            <td border="1" bgcolor="red">0</td>
            <td border="1" bgcolor="red">0</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">0</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">1</td>
            <td border="1" bgcolor="red">0</td>
            <td border="1" bgcolor="red">0</td>
            </tr></table>>
    ]
}

完整實作如下：

/* 23 bits of float fractional data */
#define I2F_FRAC_BITS  23
#define I2F_MASK ((1 << I2F_FRAC_BITS) - 1)

/*
 * Converts unsigned integer into 32-bit IEEE floating point representation.
 * Will be exact from 0 to 2^24.  Above that, we round towards zero
 * as the fractional bits will not fit in a float.  (It would be better to
 * round towards even as the fpu does, but that is slower.)
 */
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;
}

其中用到一個 `ror32` 的函式，是被定義在 [include/linux/bitops.h](https://github.com/torvalds/linux/blob/a48b0872e69428d3d02994dcfad3519f01def7fa/include/linux/bitops.h#L116)，程式碼如下：

static inline __u32 ror32(__u32 word, unsigned int shift)
{
	return (word >> (shift & 31)) | (word << ((-shift) & 31));
}

有些指令集架構有支援 rotate 指令，比方說 [x86_64 指令集](https://en.wikipedia.org/wiki/X86_instruction_listings) 具備 `ror`、`rol` 指令，上述程式碼經過編譯器最佳化，確實會產生出 `ror` 指令。

新版 `int2float` 的程式碼先使用 `__fls` 來找出該整數的指數部份。原本的實作 fraction 會用左移加上條件判斷來求值，新版使用 `(msb - I2F_FRAC_BITS) & 0x1f)`，搭配 `ror32` 來實作，減少檢查是否大於 $2^{24}$ 所需的分支，前者利用 unsigned integer 減法的溢位來決定輸入小於 $2^{24}$ 時的位移數目，搭配後者 rotate 指令，不管輸入大於或是小於 $2^{24}$，皆可以做到保留 MSB 之後的 23 個位元給 `fraction` 的效果。

## 測驗 `3`

:::success
解釋上述程式運作原理;
:::

考慮以下實作：

#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 (v);
    return u << shift;
}

性質：

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(\frac{x}{2}, \frac{y}{2})$ because $2$ is a common divisor.
4. If x is even and y is odd, $gcd(x, y) = gcd(\frac{x}{2}, y)$.
5. On similar lines, if x is odd and y is even, then $gcd(x, y) = gcd(x, \frac{y}{2})$. It is because $2$ is not a common divisor.

一開始的 if 判斷式對應到性質 1、2，如果 x、y 皆為 0，回傳 0；如果 y 為 0，回傳 x；如果 x 為 0，回傳 y。接著，for 迴圈對應到性質 3，並用 `shift` 記錄乘了多少次 2。while 迴圈對應到性質 4。do while 迴圈則是用到性質 5 以及 [Euclid's_algorithm](https://en.wikipedia.org/wiki/Greatest_common_divisor#Euclid's_algorithm)。做到最後 `u` 和 `v` 會相等，會跳到 else 區塊，做完該區塊後 `u` 為一非 0 數值，`v` 為 0，跳出迴圈，再將 `u` 左移 `shift`。

## 測驗 `4`

:::success
解釋上述程式運作原理，並舉出這樣的程式碼用在哪些真實案例中;
:::

在影像處理中，[bit array](https://en.wikipedia.org/wiki/Bit_array) (也稱 `bitset`) 廣泛使用，考慮以下程式碼:

#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) 的行為是回傳由低位往高位遇上連續多少個 `0` 才碰到 `1`。

> 範例: 當 `a = 16`
`16` 這個十進位數值的二進位表示法為 00000000 00000000 00000000 00010000
從低位元 (即右側) 往高位元，我們可發現 $0 \rightarrow 0 \rightarrow 0 \rightarrow 0 \rightarrow 1$，於是 ctz 就為 4，表示最低位元往高位元有 4 個 `0`

用以改寫的程式碼如下:

#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 = bitset & -bitset;
            int r = __builtin_ctzll(bitset);
            out[pos++] = k * 64 + r;
            bitset ^= t;                           
        }
    }
    return pos;
}

在 `naive()` 這個函式中，需要將 bitmap 中的所有 `1` 的位置記錄到 `out` 中，並回傳 `pos` 代表 `out` 的 size。其中，我們可以使用 `__builtin_ctzll` 優化，如 `improved()` 函式所示。使用 `ctz` 找到目前 bitmap 中的 least significant `1`，並將該 `1` clear，進行下一次迭代。我們可以觀察到 `bitset ^= t` 這一行程式就是 clear bit 的操作。至於 `t = bitset & -bitset` 則是利用二補數的特性 (可觀察 [解讀計算機編碼](https://hackmd.io/@sysprog/binary-representation) 的圖)，取二補數相當於對原數作 bitwise NOT 後加 1。二補數再和原數作 bitwise AND 可以得到一個數，該數的第 r 個 bit 為 set，其餘為 clear (等價於 `1 << r`)。

:::success
設計實驗，檢驗 ctz/clz 改寫的程式碼相較原本的實作有多少改進？應考慮到不同的 [bitmap density](http://www.vki.com/2013/Support/docs/vgltools-3.html);
:::

:::success
思考進一步的改進空間;
:::

可改寫成以下等價實作：

#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) {
            int r = __builtin_ctzll(bitset);
            out[pos++] = k * 64 + r;
            bitset ^= 1UL << r;                           
        }
    }
    return pos;
}

:::success
閱讀 [Data Structures in the Linux Kernel](https://0xax.gitbooks.io/linux-insides/content/DataStructures/linux-datastructures-3.html) 並舉出 Linux 核心使用 bitmap 的案例;
:::

Linux 的 [kernel/sched/sched.h](https://elixir.bootlin.com/linux/latest/source/kernel/sched/sched.h#L254)，有關 RT scheduling 就有使用到 bitmap：

/*
 * This is the priority-queue data structure of the RT scheduling class:
 */
struct rt_prio_array {
	DECLARE_BITMAP(bitmap, MAX_RT_PRIO+1); /* include 1 bit for delimiter */
	struct list_head queue[MAX_RT_PRIO];
};

搭配下述程式碼可以找出 bitmap 中的 first set：

/*
 * Every architecture must define this function. It's the fastest
 * way of searching a 100-bit bitmap.  It's guaranteed that at least
 * one of the 100 bits is cleared.
 */
static inline int sched_find_first_bit(const unsigned long *b)
{
#if BITS_PER_LONG == 64
	if (b[0])
		return __ffs(b[0]);
	return __ffs(b[1]) + 64;
#elif BITS_PER_LONG == 32
	if (b[0])
		return __ffs(b[0]);
	if (b[1])
		return __ffs(b[1]) + 32;
	if (b[2])
		return __ffs(b[2]) + 64;
	return __ffs(b[3]) + 96;
#else
#error BITS_PER_LONG not defined
#endif
}