2022q1 Homework2 (quiz2)

contributed by < ppodds >

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

EXP1

第二種寫法相當於

⌊ \frac{a}{2} ⌋ + ⌊ \frac{b}{2} ⌋ + x

很顯然我們要補的就是被floor捨去掉的值。

舉個簡單的例子

\frac{1 + 1}{2} = 1

⌊ \frac{1}{2} ⌋ + ⌊ \frac{1}{2} ⌋ = 0

我們必須在 a 和 b 同時是奇數時補 1 給最終結果。寫法如下:

(a & 1) & (b & 1)

\to

a & b & 1

故 EXP1 為 a & b & 1

EXP2 & EXP3

改寫為

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

限定使用 OR AND XOR

直覺想到用加法器的概念做

加法
$\to$ XOR
進位
$\to$ AND

先把 a + b 拆成用加法器的 bitwise 形式

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

補上除2答案就出來了

(a + b) / 2 = (a & b) + (a ^ b) >> 1

故 EXP2 為 a & b ， EXP3 為 a ^ b

延伸問題:

解釋上述程式碼運作的原理
比較上述實作在編譯器最佳化開啟的狀況，對應的組合語言輸出，並嘗試解讀 (可研讀 CS:APP 第 3 章)
研讀 Linux 核心原始程式碼 include/linux/average.h，探討其 Exponentially weighted moving average (EWMA) 實作，以及在 Linux 核心的應用

移動平均（Moving average)，又稱滾動平均值、滑動平均，在統計學中是種藉由建立整個資料集合中不同子集的一系列平均數，來分析資料點的計算方法。
移動平均通常與時間序列資料一起使用，以消除短期波動，突出長期趨勢或周期。短期和長期之間的閾值取決於應用，移動平均的參數將相應地設置。例如，它通常用於對財務數據進行技術分析，如股票價格、收益率或交易量。它也用於經濟學中研究國內生產總值、就業或其他宏觀經濟時間序列。

延伸2

比較上述實作在編譯器最佳化開啟的狀況，對應的組合語言輸出，並嘗試解讀 (可研讀 CS:APP 第 3 章)

avg.c

#include <stdint.h>
#include <stdio.h>

uint32_t average_1(uint32_t a, uint32_t b)
{
    return (a + b) / 2;
}

uint32_t average_2(uint32_t a, uint32_t b)
{
    return (a >> 1) + (b >> 1) + (a & b & 1);
}

uint32_t average_3(uint32_t a, uint32_t b)
{
    return (a & b) + ((a ^ b) >> 1);
}

int main()
{
    int a = 1, b = 1;
    printf("%d\n", average_1(a, b));
    printf("%d\n", average_2(a, b));
    printf("%d\n", average_3(a, b));
    return 0;
}

編譯指令 gcc avg.c -O2 -o avg

objdump -D -M intel avg > avg.dump

先備知識
rdi 和 rsi 是 gcc x64 calling convention 的前兩個參數

00000000000011c0 <average_1>:
    11c0:       f3 0f 1e fa             endbr64
    11c4:       8d 04 37                lea    eax,[rdi+rsi*1]
    11c7:       d1 e8                   shr    eax,1
    11c9:       c3                      ret
    11ca:       66 0f 1f 44 00 00       nop    WORD PTR [rax+rax*1+0x0]

lea 將後者 address 存的值放到前者之中，且允許後者可以進行暫存器運算

根據上面的解釋，可以看出第一行其實在執行 a + b 並把結果存入 eax 中。下一行的 shr eax,1 即是 eax / 2 最終做出 (a + b) / 2

00000000000011d0 <average_2>:
    11d0:       f3 0f 1e fa             endbr64
    11d4:       89 f8                   mov    eax,edi
    11d6:       89 f2                   mov    edx,esi
    11d8:       21 f7                   and    edi,esi
    11da:       d1 e8                   shr    eax,1
    11dc:       d1 ea                   shr    edx,1
    11de:       83 e7 01                and    edi,0x1
    11e1:       01 d0                   add    eax,edx
    11e3:       01 f8                   add    eax,edi
    11e5:       c3                      ret
    11e6:       66 2e 0f 1f 84 00 00    nop    WORD PTR cs:[rax+rax*1+0x0]
    11ed:       00 00 00

mov eax,edi 、 mov edx,esi 先各做了一份 a 和 b 變數的備份，待會算 a & b & 1 會用到。接下來直接做 a & b 存到 edi 中，並執行 a >> 1 與 b >> 1。此時才會再把 edi & 1 (在此可以看出 compiler 在優化時有重排我們的程式碼)，最後把前面計算出來的東西用 add 加起來。

00000000000011f0 <average_3>:
    11f0:       f3 0f 1e fa             endbr64
    11f4:       89 f8                   mov    eax,edi
    11f6:       21 f7                   and    edi,esi
    11f8:       31 f0                   xor    eax,esi
    11fa:       d1 e8                   shr    eax,1
    11fc:       01 f8                   add    eax,edi
    11fe:       c3                      ret
    11ff:       90                      nop

先備份 a 到 eax 後把 edi 設為 a & b ，再把 eax 設為 (a ^ b) / 2，最終把兩個結果加在一起。

延伸3

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

EWMA介紹

大致上可以理解為一個簡單的動態平均系統，一開始先設置一個權重，之後每次計算時根據權重和新進來的數值做加權平均。

include/linux/average.h 定義了一個大micro DECLARE_EWMA(name, _precision, _weight_rcp) 。可讓 user 根據傳入的參數產生 EWMA 的 struct 和操作函式。分別為 ewma_##name (struct) 、ewma_##name##_init 、 ewma_##name##_read 、 ewma_##name##_add 。

#define DECLARE_EWMA(name, _precision, _weight_rcp)			\
	struct ewma_##name {						\
		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;					\
	}								\
	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);			\
	}								\
	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));				\
	}

從上面的 code 可觀察出 linux 在 compile time 時就有對巨集傳入的參數做檢查，來避免 runtime 時的錯誤。

##name##會被代換成巨集的 name 參數

ewma_##name##_init

初始化 struct 的內容(把 internal 設成 0)

ewma_##name##_read

讀取 EWMA 當前的數值

ewma_##name##_add

加入新的一筆資料進到 EWMA 中

read 和 add 剛好對應到 EWMA 的基本操作

測驗2

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

由於題目有一點難度，需要慢慢把思考過程寫下來。

首先從EXP5開始猜測， EXP5 是一個 a 和 b 的比較，意思是只會傳回 true 或 false，即 1 或 0 。經過負號轉換後 1 變為 -1 ，即 0xFFFFFFFF 。到此為止可以稍微有點頭緒， (EXP4) & -(EXP5) 在 EXP5 的條件滿足時，會把 EXP4 的數值完整留下來，否則結果會是 0 。當此敘述為 0 時，經過和 a 的 XOR 運算，結果會是 a。固可先猜測 EXP5 在 a < b 時會滿足。

接下來考慮 EXP5 滿足的情形。當 EXP5 滿足時，我們需要把 b 作為最終結果回傳回去。簡化後可以看成 a ^ (EXP4) 要等於 b 。回想 XOR 的特性，某數對同樣數字 XOR 兩次以後會得到自己，因此可猜測 EXP4 為 a ^ b 。

EXP4 = a ^ b
EXP5 = a < b

延伸問題:

解釋上述程式碼運作的原理
針對 32 位元無號/有號整數，撰寫同樣 branchless 的實作

Linux 核心也有若干 branchless / branch-free 的實作，例如 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 檢索。

延伸2

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

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

在這裡只要 a 和 b 同時是有號數或無號數結果就會正確(不一樣的話會造成 a < b 的比較與預期不同)。

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

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$

剛好對應上

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

if x, y are both even,
$g c d (x, y) = 2 * g c d (\frac{x}{2}, \frac{y}{2})$

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

只要 u 和 v 都是偶數，就除二並把 shift 加 1 (此處用 last bit 是否等於 1 來判斷奇偶)。當迴圈結束時，會剛好達成

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

(此處的 u 和 v 指的是函式剛傳進來時的值)。

If x is even and y is odd, then
$g c d (x, y) = g c d (\frac{x}{2}, y)$

If x is odd and y is even, then
$g c d (x, y) = g c d (x, \frac{y}{2})$

此時已經不能再提出 2 的公因數了，因此可以把 u 和 v 都除 2 到變成奇數為止。

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

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

對到最後的剩下的兩條。

最後的 do while 迴圈中止條件和一般的 gcd 一樣是 v=0，故 COND = v，而回傳值如上面所述， RET = u << shift

COND = v
RET = u << shift

延伸2

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

延伸3

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

測驗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

-bitset 會將 bitset 最右邊的 1 前方的所有 bits 都翻轉(二補數)，和 bitset 做 AND 後就會只剩下最右邊的 1 的 bit。故 EXP6=bitset & -bitset 。

EXP6=bitset & -bitset

延伸1

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

延伸2

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

延伸3

思考進一步的改進空間

延伸4

閱讀 Data Structures in the Linux Kernel 並舉出 Linux 核心使用 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 = ? (表示式)

PPP=
MMM=
EEE=

測驗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 等價。

這題的答案我不確定

M=_h
X=0

將 0 轉型成 struct { char c; t _h; } 的 pointer，並取得 -h 的 address，會拿到 -h 在結構中的 offset，再扣掉一開始 struct base address 的位置。

不解的是為什麼要扣掉 base address ，如果 base address 已經是 0 了應該不需要再扣一次。

延伸2

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

延伸3

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

測驗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 = ? (表示式，不包含小括號)

KK1=div3
KK2=div5
KK3=div3 << 2

原理是透過位元運算湊出 fmt (要輸出的字串)數值

fmt 的決定方法是由 &"FizzBuzz%u"[(9 >> div5) >> (KK3)] 和 length 共同控制。第一項是字串的開頭，後面是持續長度，列出表以後可以求出上面三條式子。

9 怪怪的，要改成 8，不然吃不到 %u ，最後印出來的是 u

延伸1

解釋上述程式運作原理並評估 naive.c 和 bitwise.c 效能落差

避免 stream I/O 帶來的影響，可將 printf 更換為 sprintf

延伸2

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

延伸3

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

延伸4

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

2022q1 Homework2 (quiz2)

2022q1 第 2 週測驗題

測驗 1

EXP1

EXP2 & EXP3

延伸2

延伸3

測驗2

延伸2

測驗3

延伸2

延伸3

測驗4

延伸1

延伸2

延伸3

延伸4

測驗5

測驗6

延伸2

延伸3

測驗7

延伸1

延伸2

延伸3

延伸4

Read more

Linux Kernel Project 1

2022q1 Homework1 (lab0)