2024q1 Homework4 (quiz3+4)

contributed by < st10740 >

第三週測驗題

測驗三

這段程式碼的目的是計算以 2 為底的對數，其做法是找到該數字的二進位表示中，最高位元為 1 的索引值，而該索引值即為以 2 為底的對數值。

然而，若是透過線性搜尋 (Linear search) 的方式會找太慢，因此以下方法採取二分搜尋法 (Binary search) 的方式進行搜索。首先會先利用 while (i >= ...) 確認左半邊的位元是否存在 1，若存在則搜尋左半邊，並將 result 加上一半的位元數，若左半邊不存在為 1 的位元，則搜尋右半邊，並不斷縮小範圍直到找到為止。

static size_t ilog2(size_t i)
{
    size_t result = 0;
    while (i >= 65536) {
        result += 16;
        i >>= 16;
    }
    while (i >= 256) {
        result += 8;
        i >>= 8;
    }
    while (i >= 16) {
        result += 4;
        i >>= 4;
    }
    while (i >= 2) {
        result += 1;
        i >>= 1;
    }
    return result;
}

而另外一個作法的概念也相同，只是利用了 GCC 提供的內置函式 __builtin_clz ，其作用是計算無號整數 v 的二進位表示式中，從最高有效位開始的連續零位的数量。因此透過將 31 減去該值就可以得到最高位元為 1 的索引值。

int ilog32(uint32_t v)
{
    return (31 - __builtin_clz(v));
}

在 Linux 核心原始程式碼中有使用 log2 的相關程式碼

在 tools/include/linux/log2.h 中定義了巨集 ilog2。

這個方法會先利用 __builtin_constant_p 判斷 n 在編譯時期是否為常數，若是則會在編譯時期完成如下面 (n) & (1ULL << 63) ? 63 的運算來求得最高位為 1 的索引值，因此 ilog2 的運算結果可以在編譯時期就確定；但若不是常數則需要等到執行時期再進行運算，並且會根據 n 的大小決定呼叫 __ilog2_u32 或是 __ilog2_u64。因此這個作法可以節省輸入為常數時的執行時間。

#define ilog2(n)				\
(						\
	__builtin_constant_p(n) ? (		\
		(n) < 2 ? 0 :			\
		(n) & (1ULL << 63) ? 63 :	\
		(n) & (1ULL << 62) ? 62 :	\
		(n) & (1ULL << 61) ? 61 :	\
		(n) & (1ULL << 60) ? 60 :	\
		(n) & (1ULL << 59) ? 59 :	\
		(n) & (1ULL << 58) ? 58 :	\
		(n) & (1ULL << 57) ? 57 :	\
		(n) & (1ULL << 56) ? 56 :	\
		(n) & (1ULL << 55) ? 55 :	\

        // ...
		(n) & (1ULL << 10) ? 10 :	\
		(n) & (1ULL <<  9) ?  9 :	\
		(n) & (1ULL <<  8) ?  8 :	\
		(n) & (1ULL <<  7) ?  7 :	\
		(n) & (1ULL <<  6) ?  6 :	\
		(n) & (1ULL <<  5) ?  5 :	\
		(n) & (1ULL <<  4) ?  4 :	\
		(n) & (1ULL <<  3) ?  3 :	\
		(n) & (1ULL <<  2) ?  2 :	\
		1 ) :				\
	(sizeof(n) <= 4) ?			\
	__ilog2_u32(n) :			\
	__ilog2_u64(n)				\
 )

我進行了簡單的測試確認以上巨集處理常數的部份是否真的會在編譯時期取得運算結果，因此我拿取了巨集中處理常數的程式碼建立了一個新的巨集，並在 main 中使用該巨集，將常數 1022 作為參數輸入，並將結果印出來。

int main() {
	printf("log_2_%d: %d\n", 1022, ilog2(1022));
	return 0;
}

接著我將程式碼進行編譯，取得其編譯後的組合語言，可以觀察到第 10 行的 9 即為 ilog2(1022) 的結果，在這邊被當作常數來使用。

main:
.LFB0:
	.cfi_startproc
	endbr64
	pushq	%rbp
	.cfi_def_cfa_offset 16
	.cfi_offset 6, -16
	movq	%rsp, %rbp
	.cfi_def_cfa_register 6
	movl	$9, %edx
	movl	$1022, %esi
	leaq	.LC0(%rip), %rax
	movq	%rax, %rdi
	movl	$0, %eax
	call	printf@PLT
	movl	$0, %eax
	popq	%rbp
	.cfi_def_cfa 7, 8
	ret
	.cfi_endproc

測驗四

Exponentially Weighted Moving Average (EWMA; 指數加權移動平均) 是一種取平均的作法，在訊號處理中可以作為濾波器使用，其公式如下：

$$S_t=\{\begin{array}{ll}{Y}_0,& t=0\\ \alpha Y_t+(1-\alpha )S_{t-1},& t>0\end{array}$$
可以透過

$\alpha$ 決定新資料與歷史資料的權重，介於 0 到 1 之間，

$\alpha$ 值越大代表歷史資料的權重減少地越快。

在實作方面，首先宣告一個結構體 ewma ，internal 用來儲存當前算出來的平均值；weight 用來存新資料與歷史資料的的權重，

$\frac{1}{weight}$ 等同於上面公式中的

$\alpha$ ；factor 的目的是保留平均數的小數部份，其值存的是欲保留的小數點位數，因為 internal 的型別為 unsigned long，進行平均值運算時可能會丟失小數部份，因此透過 factor 來避免誤差。

struct ewma {
    unsigned long internal;
    unsigned long factor;
    unsigned long weight;
};

ewma 的初始化確保傳入的 weight 和 factor 都是 2 的冪，原因是為了提升效能，會使用位元運算子進行乘除運算，接著分別取得 weight 和 factor 欲左移或右移的位元數，並將 interval 初始化為 0。

void ewma_init(struct ewma *avg, unsigned long factor, unsigned long weight)
{
    if (!is_power_of_2(weight) || !is_power_of_2(factor))
        assert(0 && "weight and factor have to be a power of two!");

    avg->weight = ilog2(weight);
    avg->factor = ilog2(factor);
    avg->internal = 0;
}

ewma_add 用以計算加入新資料後的整體平均值，當資料為第一筆時則直接把值當作平均值，即

$S_0=Y_0$ ，對應到程式碼中的第 6 行，注意到這邊有對 val 進行左移 avg->factor 個位元的操作，目的是預先保留 avg->factor 個位元當作小數部分，避免之後進行平均值運算時丟失小數資訊。

而當資料不為第一筆時，則需計算新資料與歷史資料的權重，即

$S_t=\alpha Y_t+(1-\alpha )S_{t-1}$，對應到程式碼的第 4, 5 行，這邊針對新資料 val 也有先左移 avg->factor 個位元，另外觀察這行表示式 (((avg->internal << avg->weight) - avg->internal) + (val << avg->factor)) >> avg->weight 可以拆成：((avg->internal << avg->weight) - avg->internal) >> avg->weight 和 (val << avg->factor) >> avg->weight 相加，分別對應到

$(1-\alpha )S_{t-1}$ 和

$\alpha Y_t$ 的運算。

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

因為 EWMA 的計算過程中為了保留小數部分有將新資料都先左移 avg->factor 個位元再進行運算，因此若要讀取計算完成的 EWMA，需要將 avg->interval 右移 avg->factor 個位元校正回來。

unsigned long ewma_read(const struct ewma *avg)
{
    return avg->internal >> avg->factor;
}

在 Linux 核心原始程式碼中有使用 EWMA 的相關程式碼

在 /include/linux/average.h 中定義了巨集 DECLARE_EWMA 用以進行 EWMA 的運算。裡面定義的函式 ewma_##name##_init、ewma_##name##_read 和 ewma_##name##_add 都有進行以下的確認才進行後續的操作：

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

首先透過 BUILD_BUG_ON(!__builtin_constant_p(...)) 可以確保 _precision 和 _weight_rcp 在編譯時期都是常數。

BUILD_BUG_ON((_precision) > 30) 用以確保 _precision 的值不會大於 30，根據註解說明其目的是為了避免過多位元都被當作小數部分，不禁好奇為甚麼是使用 30 作為分界呢？ 但是翻過歷史 commit 紀錄，都沒有說明使用 30 的原因。

另外我發現，在 Linux 核心原始程式碼中的 _precision 功能與測驗題的 factor 相同，但是比起 factor，_precision 的命名方式更可以讓我明白它的用途，而 commit eb1e011 做的變更就是為了這個目的把 _factor 改成 _precision。

測驗五

測驗三的 ilog2 計算的是

$\lfloor lo{g}_2(x)\rfloor$ ，而測驗五的 ceil_ilog2 計算的是

$\lceil lo{g}_2(x)\rceil$。在作法上，ilog2 需要求該數字的二進位表示中，最高位元為 1 的索引值，其索引從 0 開始；而 ceil_ilog2 的作法則需要分成兩種情況探討，一種是其值為 2 的冪的情況，在這個情況下作法同 ilog2，另外一種則是值不為 2 的冪，其作法同 ilog2 但是最後需要再加上 1 。

一開始 x-- 的操作用以處理 x 為 2 的冪的情況，使其不再是 2 的冪，因此便可以採取處理不是 2 的冪的方式進行計算，且最終也會得到正確的結果，因此接下來的操作適用於任何大於 0 的整數。

接下來的操作與 ilog2 找最高位元為 1 的方式大同小異，都是採取二分搜尋法，r 用來存最高位的索引值，不過與 ilog2 不同的地方在於，這邊都是使用位元運算子來進行操作，可以達到 branchless ，提高效能。

int ceil_ilog2(uint32_t x)
{
    uint32_t r, shift;

    x--;
    r = (x > 0xFFFF) << 4;                               
    x >>= r;
    shift = (x > 0xFF) << 3;
    x >>= shift;
    r |= shift;
    shift = (x > 0xF) << 2;
    x >>= shift;
    r |= shift;
    shift = (x > 0x3) << 1;
    x >>= shift;
    return (r | shift | x > 0x1) + 1;       
}