2022q1 Homework2 (quiz2)

contributed by < laneser >

作業要求

2022q1 第 2 週測驗題

測驗 1

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

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

改寫為

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

想法 & 思考

看到 a >> 1 就會直覺想到是結果相等於把 a 除以 2. 但會無條件捨去最後一個 bit.
因此 EXP1 應該就是彌補了 a, b 都是無條件捨去的情況下，結果跟 (a + b) / 2 的差異。
也就是 EXP1 需要利用 a, b 的最後一個 bit 產生一個彌補值，相當於 (a + b) / 2 跟 (a >> 1) + (b >> 1) 之間的差異。
所以可以列出 a, b 最後一個 bit 跟我們期望的答案的真值表。

a lowest bit	b lowest bit	EXP1 result
0	0	0
0	1	0
1	0	0
1	1	1

那就很明顯是這兩個 bits 的 AND 結果。
所以答案就是 (a & 1) & (b & 1), 簡化為 a & b & 1。

改寫為

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

題目要求只能用 AND, OR, XOR 三種 bitwise 運算子，也就意味著這三種方法跟 +,- 有一定的關聯。只好求救 Google, 找到 Bitwise XORing two numbers results in sum or difference of the numbers, 提到
x + y = (x ^ y) + ((x & y) << 1) , 這個等式就可以聯想到我們要求的答案是 (a + b) >> 1 , 所以將這個等式兩邊都除以 2, 也就是 bitwise 右移一位。(x + y) >> 1 = (x ^ y) >> 1 + (x & y), 這樣就跟我們要的答案非常接近了。EXP2 應該就是 a & b, 而 EXP3 就是 a ^ b.

延伸問題

嘗試解讀編譯器最佳化開啟的狀況，對應的組合語言輸出。

利用 gcc -S -O3 test1.c 可以得到每個 average 實作函式對應的組合語言輸出。
以下都移除 assembly directives, 因為那是給編譯器看的，按題目需求應該是專注於給 CPU 執行的指令 (instructions).

(a + b) / 2

; uint32_t average(uint32_t a, uint32_t b) {return (a + b) / 2;}
average:
	endbr64
	leal	(%rdi,%rsi), %eax
	shrl	%eax
	ret

endbr64 根據 What does the endbr64 instruction actually do? 解釋為 End for branch 64 bit. 在後續的其他函式也都是一開始就執行這個指令，在這裡也可以解釋為 nop, 沒有任何動作。

leal 這個指令是 load effactive address 的 64bits 版本，但常常用來取代 add 指令, What's the purpose of the LEA instruction? 可以理解為，為了支持高階語言的陣列功能，x86 指令集特別弄了這個，但被發現跟 add 有異曲同工之妙，但又比 add 好用。
總之在這裡，就意味著把傳進來的兩個參數 a,b, 分別放在 esi,edi 這兩個暫存器內，直接加起來把結果存入 eax 暫存器。

shrl 很好理解，就是右移一個 bit. 把 eax 暫存器內容除以 2 的意思。
最後回傳值放在 eax.

low + (high - low) / 2

;uint32_t average2(uint32_t a, uint32_t b) {return a + (b-a) / 2;}
average2:
	endbr64
	movl	%esi, %eax
	subl	%edi, %eax
	shrl	%eax
	addl	%edi, %eax
	ret

前兩行就是把 b-a 的結果放在 eax, 再右移一個 bit. 最後加上 a. 所以如果 b-a 的結果會 overflow, 那最後的結果也就不正確了。

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












;uint32_t average_EXP1(uint32_t a, uint32_t b){return (a >> 1) + (b >> 1) + (a & b & 1);}
average_EXP1:
	endbr64
	movl	%edi, %eax
	movl	%esi, %edx
	andl	%esi, %edi
	shrl	%eax
	shrl	%edx
	andl	$1, %edi
	addl	%edx, %eax
	addl	%edi, %eax
	ret

第 4,7 行實作了 b >> 1 , 第 5,8 行實作了 a >> 1, 而這兩個結果分別放在 eax,edx.
第 6,9 行實作了 a & b & 1, 結果存放於 edi.
最後 10,11 行把這三個結果加起來回傳。

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

;uint32_t average_EXP2_EXP3(uint32_t a, uint32_t b){return (a & b) + ((a ^ b) >> 1);}
average_EXP2_EXP3:
	endbr64
	movl	%edi, %eax
	andl	%esi, %edi
	xorl	%esi, %eax
	shrl	%eax
	addl	%edi, %eax
	ret

前三行實作了 a & b, a ^ b, 分別放在 edi,eax.
接著再把 a ^ b 右移一個 bit.
最後加起來回傳。

研讀 Linux 核心原始程式碼 `include/linux/average.h`，探討其 Exponentially weighted moving average (EWMA) 實作，以及在 Linux 核心的應用

先了解 EWMA 的定義:

S_{t} = {\begin{cases} Y_{0}, & t = 0 \\ α Y_{t} + (1 - α) \cdot S_{t - 1}, & t > 0 \end{cases}

average.h 定義了一個大 macro, 共有三個參數, 分別是

name: 表示整個 macro 定義出來的一組 macro 以及 struct 的名稱。
_percision: 表示用多少個 bits 來表示內部儲存的小數部分。
_weight_rcp: 表示權重
$α$ 的倒數,
$α = 1 / w e i g h t_r c p$ 。

接著就是拆解整個 macro 的每一部分:

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

表示這個 macro 會定義出一個 ewma_ 開頭, name 參數結尾的 struct, 內部只有一個 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;
}

在使用之前，需要先呼叫 ewma_##name##_init 這個 incline function.
會在 compiler 時期就檢查所有該檢查的參數，而真正的 code 只有 e->internal = 0, 看來 linux kernel 傾向一切能夠在 compiler 時就能做的就盡量做，這樣可以最早發現問題，而且又確保了執行期的速度與正確性。

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

看到這裡我有個小疑問，難道這四行檢查碼會因為放在不同的 macro 當中而有所不同嗎?
看起來都一樣的程式為什麼要寫第二遍呢?
我想像中，只要寫一次，應該會在 compiler 時期就會發生錯誤警告了。
不過，做完檢查就只剩下一行, 也就是讀取值的時候直接把 internal 右移 _precision 這麼多個 bits 後回傳。

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

這個 macro 就是最主要的功能。
其中利用了 READ_ONCE 來讀取 e->internal, 也利用 WRITE_ONCE 來寫入 e->internal, 看來是為了避免在不同的 thread 當中存取所發生的問題。
另外 ilog2 函式看來也是個 macro 可以在 compiler time 就可以運算出 log2 的值。
算出 log2 的值就可以運用 bits 位移加速乘法以及除法運算。
為了避免 double evaluation, _preision 先指派給 precision 取值後再使用， macro 的世界都要小心翼翼的。

所以最後一行的意思是 (用 _val 取代 val << precision)

internal 沒有值的時候: e->internal = _val
internal 有值的時候: e->internal = (internal*(_weight_rcp - 1) + _val)/_weight_rcp
如果把
$_w e i g h t_r c p = 1 / α$
e->internal = (
$i n t e r n a l * (1 / α - 1) +_v a l) * α$
e->internal =
$i n t e r n a l * (1 - α) +_v a l * α$
這正是 Exponentially weighted moving average (EWMA) 的定義。

探究這個 DECLARE_EWMA 在 linux kernel 當中的應用，最快的方法就是搜尋它被使用在哪裡。

用在 network or network driver, 推論跟 ethernet 有關
- net/mac80211/sta_info.h
- drivers/net/wireless/mediatek/mt76/mt76x02_mac.h
- drivers/net/wireless/realteck/rtw88/main.h
- drivers/net/wireless/mediatek/mt7601u/mt7601u.h
- …還有十幾處都是 net or drivers/net 下方的實作
gpu
- drivers/gpu/drm/i915/gt/intel_context_types.h
- drivers/gpu/drm/drm_self_refresh_helper.c

以網路來說, 參考 sta_info.h 的註解說明, 這個是 average ACK signal 的意思, 在每次 ack 的時候都會更新一次.

Google 文章說是用在估算 Round Trip Time (RTT), 相關資料在 RFC793 page 40 當中提到 Retransmission Timeout 需要動態被估算, 而且需要得到 RTT 這個資訊, 根據這份簡報提到

\begin{array}{r} E s t i m a t e d R T T = (1 - g) (E s t i m a t e d R T T) + g (S a m p l e R T T) \\ 0 \leq g \leq 1, u s u a l l y g = .1 o r .2 \end{array}

這個公式就是 EWMA.

而 gpu 似乎是跟 DRM (Direct Rendering Manager) 有關.

測驗 2

改寫〈解讀計算機編碼〉一文的「不需要分支的設計」一節提供的程式碼 min，我們得到以下實作 (max):
#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b)
{
    return a ^ ((EXP4) & -(EXP5));
}

想法 & 思考

既然是取比較大的值, 結果只會有兩種: a , b.
若是 a 比較大, 則結果應該類似 a ^ 0 = a, 若是 b 比較大, 結果應該類似 a ^ (a ^ b) = b.
所以 ((EXP4 & -(EXP5)) 的結果應該會如下:

if a >= b then ((EXP4 & -(EXP5)) = 0
if a < b then ((EXP4 & -(EXP5)) = (a ^ b)

這時可以想到 EXP4 應該是 a ^ b , 而 -(EXP5) 應該如下:

if a >= b then -(EXP5) = 0b 使得 (a ^ b) & 0 = 0
if a < b then -(EXP5) = 0xFFFFFFFF = -1 使得 (a ^ b) & 0xFFFFFFFF = (a ^ b)

因此 EXP5 就是 a < b

延伸問題

針對 32 位元無號/有號整數，撰寫同樣 branchless 的實作

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

不管 a , b 是否為 signed 或 unsigned, -(a < b) 的結果就是 signed int, 且為 0x000000 或者 0xFFFFFFFF 這兩種結果之一, 所以實作的方法是一樣的。

請在 Linux 核心原始程式碼中，找出更多類似 branchless 的實作手法。請善用 git log 檢索。

kvm/mmu.h from commits

/*
 * If CPL < 3, SMAP prevention are disabled if EFLAGS.AC = 1.
 *
 * If CPL = 3, SMAP applies to all supervisor-mode data accesses
 * (these are implicit supervisor accesses) regardless of the value
 * of EFLAGS.AC.
 *
 * This computes (cpl < 3) && (rflags & X86_EFLAGS_AC), leaving
 * the result in X86_EFLAGS_AC. We then insert it in place of
 * the PFERR_RSVD_MASK bit; this bit will always be zero in pfec,
 * but it will be one in index if SMAP checks are being overridden.
 * It is important to keep this branchless.
 */
unsigned long smap = (cpl - 3) & (rflags & X86_EFLAGS_AC);

net/xfrm.h from commits

static inline bool addr4_match(__be32 a1, __be32 a2, u8 prefixlen)
{
	/* C99 6.5.7 (3): u32 << 32 is undefined behaviour */
	if (sizeof(long) == 4 && prefixlen == 0)
		return true;
	return !((a1 ^ a2) & htonl(~0UL << (32 - prefixlen)));
}

這裡的 if 是為了避免產生 u32 << 32 這種未定義的行為。
而且這次 commit 是從組合語言角度說明了 Branch savings: 1.

Rework 64 bits shifts to avoid tests and branches - 但這個我真的看不懂

還有以前的一些紀錄在 git logs 裡面, 但目前的 linux kernel 似乎已經不存在了:

測驗 3

考慮以下 64 位元 GCD (greatest common divisor, 最大公因數) 求值函式:

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

改寫為以下等價實作:

#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 演算法，可以得到很好的理解。

if (!u || !v) return u | v;

就實作了兩個步驟:

if both x, y are 0,
$g c d (0, 0) = 0$
$g c d (x, 0) = x$ and
$g c d (0, y) = y$

int shift;
for (shift = 0; !((u | v) & 1); shift++) {
    u /= 2, v /= 2;
}

這個步驟就是
3. if x, y are both even,

g c d (x, y) = 2 * g c d (\frac{x}{2}, \frac{y}{2})

在迴圈執行完畢的時候, !((u | v) & 1) 是 false, 意味著 (u | v) 的最後一個 bit 是 1, 也就是 u 或 v 其中一個是奇數. 或者兩者都是奇數。
而此時的

2^{s h i f t} * g c d (u, v)

就是原本要求的 GCD。

接著的

while (!(u & 1))
    u /= 2;

以及

while (!(v & 1))
    v /= 2;

因為此時必定不會存在 u, v 同時為偶數的狀況，所以其中一個是偶數就可以立刻把 2 這個非公因數 給剔除。
相當於演算法當中的
4. If x is even and y is odd, then

g c d (x, y) = g c d (\frac{x}{2}, y)

以及
5. If x is odd and y is even, then

g c d (x, y) = g c d (x, \frac{y}{2})

最後迴圈中的

if (u < v) {
    v -= u;
} else {
    uint64_t t = u - v;
    u = v;
    v = t;
}

是利用輾轉相除法的原理，只是用了連續減法來找出相除的餘數結果。
輾轉相除法就是除到餘數為 0 就是答案。
最後一步 u = v 的時候, 就是餘數為 0 的時候，此時 uint64_t t = u - v 會讓 t 為 0, 然後被指派給 v.
所以 COND 就是餘數 v , 而 u 在這裡就是最後的那個最大公因數。
但要記得真正的答案是

2^{s h i f t} * g c d (u, v)

.
因此 RET 就是 u << shift.

延伸問題

在 x86_64 上透過 __builtin_ctz 改寫 GCD，分析對效能的提升;

在理解了 __builtin_ctz 的意義後，顯然是針對

while (!(u & 1))
    u /= 2;

這樣的程式碼，可以直接等於

u >> __builtin_ctz(u)

所以實作的程式碼為

uint64_t gcd64_v3(uint64_t u, uint64_t v)
{
    if (!u || !v)
        return u | v;
    int u_ctz = __builtin_ctzl(u);
    int v_ctz = __builtin_ctzl(v);
    u = u >> u_ctz;
    v = v >> v_ctz;
    int shift = u_ctz > v_ctz ? v_ctz : u_ctz;
    do
    {
        v = v >> __builtin_ctzl(v);
        if (u < v)
        {
            v -= u;
        }
        else
        {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while (v);
    return u << shift;
}

可以連續呼叫後，計算耗費的時間來分析效能:

float execute_cost(gcd_function func, int times)
{
    clock_t start = clock();
    while (times--)
    {
        func(rand64(), rand64());
    }
    clock_t end = clock();
    return (float)(end - start) / CLOCKS_PER_SEC;
}

在我的環境中執行 1000000 次的結果為

gcd64_v2 cost 0.474353 seconds, __builtin_ctz cost 0.203850 seconds,
Uplift performance 132.697067%

大約快了 1.3 倍左右。

Linux 核心中也內建 GCD (而且還不只一種實作)，例如 lib/math/gcd.c，請解釋其實作手法和探討在核心內的應用場景。

lib/math/gcd.c 當中主要的 gcd 實作如下:

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

這段程式碼主要由 Zeng Zhaoxiu 所提的 commit 貢獻的.
他在這段 commit log 當中直接把 benchmark 的程式跟結果展示出來，比原本的速度快了很多。

提到了加速的主因，就是 % 這個求餘數在 CPU 指令集是用除法做的，即使是 x86 這樣的已經針對除法優化的 CPU, 用減法搭配檢查 2 的因數依舊快了 25%, 更別提 Alpha 或者 ARMv6 這種僅僅是利用 emulation code 來實作除法的 CPU. 加速效果更是明顯。

There are two variants of the code here, depending on whether a fast
__ffs (find least significant set bit) instruction is available.

也提到了 __ffs 是利用 CPU 提供的指令，意義上等同 __builtin_ctz 的實作。

其中的 r & -r 讓我思考蠻久的，是把 r 跟 r 的二進位補數做 bitwise AND 運作，這會使得 r 僅保留從低位元數過來的第一個 1 bit，意義上就是 1 << __builtin_ctz(r).

剩下的部分就幾乎跟本次的測驗題目一樣了。

核心內對於 gcd 的應用場景

對於頻率相關運算，用在影音或無線網路等等，這應該是主要的應用了
- pcm_timer.c
- adau-utils.c
- aptina-pll.c
- iio-rescale.c
- amdgpu_afmt.c
- …
數學
- 最小公倍數 lcm.c 可以利用 gcd 實作
Load balancer IPVS
- ip_vs_mh.c
檔案系統中的 blk
- blk-settings.c
tty
- vt.c 捲動文字螢幕也要用 gcd ?!
其他我不懂的
- sparx5_calendar.c
- ecc.c 這個註解提到的 link 已經 404 了, 可以考慮用 From Euclid’s GCD to Montgomery
  Multiplication to the Great Divide 取代嗎，可憐的 sun 被 oracle 買下來後連文件都保留不了，不過這裡只是提到 gcd, 似乎沒使用 gcd.

測驗 4

考慮以下程式碼:

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

改寫的程式碼如下:
















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

其中第 9 行的作用是找出目前最低位元的 1，並紀錄到 t 變數。舉例來說，若目前 bitset 為 000101b，那 t 就會變成 000001b，接著就可以透過 XOR 把該位的 1 清掉，其他保留 (此為 XOR 的特性)。
請補完 EXP6

想法 & 思考

看到這裡剛好跟前一題在 linux kernel 中的 gcd 用到的 r & -r 好像。
所以就是 bitset & -bitset.

延伸問題

舉出這樣的程式碼用在哪些真實案例中

Linux kernel gcd implement.

設計實驗，檢驗 ctz/clz 改寫的程式碼相較原本的實作有多少改進？應考慮到不同的 bitmap density;

int main()
{
    srand(time(NULL));
    int bitmapsize = 100000;
    int maxoutsize = bitmapsize * sizeof(uint64_t) * 8;
    uint64_t *bitmap = malloc(sizeof(uint64_t) * bitmapsize);
    uint32_t *naive_out = malloc(sizeof(uint32_t) * maxoutsize);
    uint32_t *improved_out = malloc(sizeof(uint32_t) * maxoutsize);
    for (int i = 0; i < bitmapsize; i++)
        bitmap[i] = rand64();
    for (int i = 0; i < maxoutsize; i++)
        naive_out[i] = improved_out[i] = 0;
    clock_t t1 = clock();
    naive(bitmap, bitmapsize, naive_out);
    clock_t t2 = clock();
    improved(bitmap, bitmapsize, improved_out);
    clock_t t3 = clock();
    float naive_cost = (float)(t2 - t1) / CLOCKS_PER_SEC;
    float improved_cost = (float)(t3 - t2) / CLOCKS_PER_SEC;
    int isPass = 1;
    for (int i = 0; i < maxoutsize; i++)
    {
        if (naive_out[i] != improved_out[i])
        {
            printf("Mismatched at index %d, %u != %u\r\n", i, naive_out[i], improved_out[i]);
            isPass = 0;
            break;
        }
    }
    printf(isPass ? "Passed\r\n" : "Failed\r\n");
    printf("naive cost %f seconds, improved cost %f seconds,\r\nUplift performance %f%%\r\n",
           naive_cost, improved_cost, (1 / improved_cost - 1 / naive_cost) * 100 / (1 / naive_cost));
    free(bitmap);
    free(naive_out);
    free(improved_out);
    free(improved2_out);
}

執行結果

naive cost 0.020146 seconds, improved cost 0.002941 seconds,
Uplift performance 585.005066%

思考進一步的改進空間;

稍微寫了一個簡單的 mapping 機制來取代計算 bits, 但果然程式碼變大又變慢。
用 out[pos++] = (k << 6) + r; 來取代 out[pos++] = k * 64 + r; 也沒有比較快。
還沒有想到更多的改進空間。

閱讀 Data Structures in the Linux Kernel 並舉出 Linux 核心使用 bitmap 的案例;

lib/iommu-helper.c
tools/perf/util/pmu.y
drivers/firmware/smccc/kvm_guest.c
linux/linkmode.h
非常多且雜, 但基本上只要想要省記憶體, bitmap 是非常好用的資料結構。

測驗 5

以下是 LeetCode 166. Fraction to Recurring Decimal 的可能實作:

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

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

請補完
PPP = ? (表示式)
MMM = ? (巨集)
EEE = ? (表示式)

想法 & 思考

find 函式實作就像 Hash Map, 在 heads 當中用 key 去找到對應的 index.
所以 int pos = find(heads, size, remainder); 當 pos >= 0 就是找到了循環的小數起點, 也就是之前出現過相同的餘數.
而 decimal 放的是之前算出來的餘數除法結果, 所以

while (PPP > 0)
    *p++ = *decimal++;

就是要把非循環的數字放進 p, 直到之前找到的餘數為止, 所以 PPP 就跟 pos 這個數字有關, 這個迴圈應該要執行 pos 這麼多次, PPP 應該是 pos--

而 MMM(&node->link, EEE); 顯然就是沒找到的時候要建立 Hash Map, 此時就是把 &node->link 塞入 heads 這個 Hash Map, Hash function 就在 find 裡面的 int hash = key % size;.
因此 MMM 就是 list_add,
EEE 就是 &heads[remainder % size]

延伸問題

解釋上述程式碼運作原理，指出其中不足，並予以改進

第一眼看 fractionToDecimal 最不高興的就是只有 malloc 卻沒有 free.
所以我們可以寫一個 free head_head Array 的函式:

void freeListArray(struct list_head *heads, size_t size)
{
    for (int i = 0; i < size; i++)
    {
        struct list_head *n, *s;
        list_for_each_safe(n, s, &heads[i])
        {
            free(container_of(n, struct rem_node, link));
        }
    }
    free(heads);
}

老師說判斷正負號也可以用更好的方法:

bool sign = (float)numerator / denominator >= 0;
if (!sign)
    *p++ = '-';

改為

if ((numerator < 0) ^ (denominator < 0))
    *p++ = '-';

既然都已經 /* deal with negtive cases */ , 那麼 long long division = n / d; 就不會有負的問題

sprintf(p, "%ld", division > 0 ? (long)division : (long)-division);

改為

sprintf(p, "%ld", (long)division);

利用 llabs 簡化:

long long n = llabs(numerator);
long long d = llabs(denominator);

可以考慮不需要先輸出到 decimal, 而是直接輸出給結果, 找到重複的 reminder 的時候, 再來插入 ( 跟 ) 就好. 最後如下:

char *fractionToDecimal2(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 = llabs(numerator);
    long long d = llabs(denominator);
    if ((numerator < 0) ^ (denominator < 0))
        *p++ = '-';

    long long remainder = n % d;
    long long division = n / d;

    sprintf(p, "%ld", (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,
     */
    int listsize = 1333;
    struct list_head *heads = malloc(listsize * sizeof(*heads));
    for (int i = 0; i < listsize; i++)
        INIT_LIST_HEAD(&heads[i]);
    char *idx_begin = p;
    int i = 0;
    for (; remainder; i++)
    {
        int pos = find(heads, listsize, remainder);
        if (pos >= 0)
        {
            size_t len = i - pos;
            memmove(idx_begin + pos + 1, idx_begin + pos, len);
            idx_begin[pos] = '(';
            idx_begin[pos + len + 1] = ')';
            idx_begin[pos + len + 2] = '\0';
            freeListArray(heads, listsize);
            return result;
        }
        struct rem_node *node = malloc(sizeof(*node));
        node->key = remainder;
        node->index = i;

        list_add(&node->link, &heads[remainder % listsize]);

        idx_begin[i] = (remainder * 10) / d + '0';
        remainder = (remainder * 10) % d;
    }
    idx_begin[i] = '\0';
    freeListArray(heads, listsize);
    return result;
}

在 Linux 核心原始程式碼的 mm/ 目錄 (即記憶體管理) 中，找到特定整數比值 (即 fraction) 和 bitwise 操作的程式碼案例，予以探討和解說其應用場景

vmscan.c - 看不懂

但這裡讓我學到 github 搜尋可以指定目錄, 像是搜尋 mm/ 目錄中的 fraction

測驗 6

alignof 是 GNU extension，以下是其可能的實作方式:
/*
 * 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)
請補完上述程式碼，使得功能與 alignof 等價。

想法 & 思考

在只有定義一種 type t 的情況下, alignof 就是 t 的長度.
既然要算出 type t 的長度, 頭尾相減就是答案.
所以 M 表示的就是 type t 的變數: _h,
而 X 表現的就是頭, 就是 0.

延伸問題

在 Linux 核心原始程式碼中找出 alignof 的使用案例 2 則，並針對其場景進行解說

socket.h 利用 __alignof__ 找出 struct sockaddr, 在定義 socket storage 的時候使用.
signal_compat.c 利用 __alignof__ 強制在 compiler 時期就拒絕做出不合規格的 siginfo_t

在 Linux 核心源使程式碼找出 ALIGN, ALIGN_DOWN, ALIGN_UP 等巨集，探討其實作機制和用途，並舉例探討 (可和上述第二點的案例重複)。思索能否對 Linux 核心提交貢獻，儘量共用相同的巨集

ALIGN, ALIGN_DOWN, ALIGN_UP 應該是在 align.h

測驗 7

考慮貌似簡單卻蘊含實作深度的 FizzBuzz 題目:

從 1 數到 n，如果是 3的倍數，印出 “Fizz”
如果是 5 的倍數，印出 “Buzz”
如果是 15 的倍數，印出 “FizzBuzz”
如果都不是，就印出數字本身

以下是利用 bitwise 和上述技巧實作的 FizzBuzz 程式碼: (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;
}

作答區
KK1 = ? (變數名稱)
KK2 = ? (變數名稱)
KK3 = ? (表示式，不包含小括號)

想法 & 思考

按題意搭配 strncpy(fmt, &"FizzBuzz%u"[(9 >> div5) >> (KK3)], length);
length 應該

div3 = 1, div5 = 0: length = 4
div3 = 0, div5 = 1: length = 4
div3 = 1, div5 = 1: length = 8
div3 = 0, div5 = 0: length = 2

因此 KK1 給 div3, KK2 給 div5 就會符合 length 的需求.
而 (9 >> div5) >> (KK3) 按需求則是

div3 = 1, div5 = 0 : (9 >> div5) >> (KK3) = 0
div3 = 0, div5 = 1 : (9 >> div5) >> (KK3) = 4
div3 = 1, div5 = 1 : (9 >> div5) >> (KK3) = 0
div3 = 0, div5 = 0 : (9 >> div5) >> (KK3) = 8

因此 9 要改為 8, 然後 KK3 等於 div3 << 2 就會符合需求…

延伸問題

解釋上述程式運作原理並評估 naive.c 和 bitwise.c 效能落差
避免 stream I/O 帶來的影響，可將 printf 更換為 sprintf

分析 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

研讀 The Fastest FizzBuzz Implementation 及其 Hacker News 討論，提出 throughput (吞吐量) 更高的 Fizzbuzz 實作

解析 Linux 核心原始程式碼 kernel/time/timekeeping.c 裡頭涉及到除法運算的機制，探討其快速除法的實作 (注意: 你可能要對照研讀 kernel/time/ 目錄的標頭檔和程式碼)
過程中，你可能會發現可貢獻到 Linux 核心的空間，請充分討論

2022q1 Homework2 (quiz2)

2022q1 第 2 週測驗題

測驗 1

想法 & 思考

延伸問題

(a + b) / 2

low + (high - low) / 2

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

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

研讀 Linux 核心原始程式碼 include/linux/average.h，探討其 Exponentially weighted moving average (EWMA) 實作，以及在 Linux 核心的應用

測驗 2

想法 & 思考

延伸問題

測驗 3

想法 & 思考

延伸問題

測驗 4

想法 & 思考

延伸問題

測驗 5

想法 & 思考

延伸問題

測驗 6

想法 & 思考

延伸問題

測驗 7

想法 & 思考

延伸問題

Read more

2022q1 Homework3 (fibdrv)

2022q1 Homework3 (quiz3)

2022q1 Homework1 (lab0)

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

研讀 Linux 核心原始程式碼 `include/linux/average.h`，探討其 Exponentially weighted moving average (EWMA) 實作，以及在 Linux 核心的應用