2022q1 Homework2 (quiz2)

contributed by < linjohnss >

測驗題目

測驗 1 無號整數取平均值

程式運作原理

兩個無號整數取平均，代表在實做過程中我們不需要考慮 sign bit 的問題，但仍需考慮 overflow 的問題。

第一個等價實作

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

從上述程式法可以看到它先將 a 與 b 右移一個位元，代表分別將 a 與 b 除 2。接著將兩者加在一起。

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

再來，我們需要考慮 a、b 是否為奇數相加，若果都是奇數則需考慮進位。
因此，EXP1 可以寫成 a & b & 1 來處理奇數的進位。

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

我們可以透過帶入數字來演示

a = 1001 (9)
b = 1111 (15)
a >> 1 = 0100 (4)
b >> 1 = 0111 (7)
0100 + 0111 = 1011 (11) 缺少 last bit 的進位
a & b & 1 = 0001 都為 1 (奇數)故需要進位
1011 + 0001 = 1100 (12)

第二個等價實做

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

題目這邊要求只能使用 AND、OR、XOR 來實作
若用 c 來實做半加器
SUM = (a ^ b)
CIN = (a & b) (CIN 需要 *2 達到進位)
a + b = (a ^ b) + ((a & b) << 1)

接著將 a + b 除以 2，可以推得
(a + b) >> 1= (a ^ b) >> 1 + (((a & b) << 1) >> 1) = (a ^ b) >> 1 + (a & b)
因此，得知 EXP2 = a & b 以及 EXP3 = a ^ b。

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

開啟編譯器最佳化

編譯器版本:gcc version 9.3.0
使用 -O3 最高級別優化
gcc -S -O3 interger_average.c

第一個等價實作

優化前：

average1:
.LFB0:
	.cfi_startproc
	endbr64
	pushq	%rbp
	.cfi_def_cfa_offset 16
	.cfi_offset 6, -16
	movq	%rsp, %rbp
	.cfi_def_cfa_register 6
	movl	%edi, -4(%rbp)
	movl	%esi, -8(%rbp)
	movl	-4(%rbp), %eax
	shrl	%eax
	movl	%eax, %edx
	movl	-8(%rbp), %eax
	shrl	%eax
	addl	%eax, %edx
	movl	-4(%rbp), %eax
	andl	-8(%rbp), %eax
	andl	$1, %eax
	addl	%edx, %eax
	popq	%rbp
	.cfi_def_cfa 7, 8
	ret
	.cfi_endproc

優化後：

average1:
.LFB0:
	.cfi_startproc
	endbr64
	movl	%edi, %eax
	movl	%esi, %edx
	andl	%esi, %edi
	shrl	%eax
	shrl	%edx
	andl	$1, %edi
	addl	%edx, %eax
	addl	%edi, %eax
	ret
	.cfi_endproc

第二個等價實作

優化前：

average2:
.LFB1:
	.cfi_startproc
	endbr64
	pushq	%rbp
	.cfi_def_cfa_offset 16
	.cfi_offset 6, -16
	movq	%rsp, %rbp
	.cfi_def_cfa_register 6
	movl	%edi, -4(%rbp)
	movl	%esi, -8(%rbp)
	movl	-4(%rbp), %eax
	andl	-8(%rbp), %eax
	movl	%eax, %edx
	movl	-4(%rbp), %eax
	xorl	-8(%rbp), %eax
	shrl	%eax
	addl	%edx, %eax
	popq	%rbp
	.cfi_def_cfa 7, 8
	ret
	.cfi_endproc

優化後：

average2:
.LFB1:
	.cfi_startproc
	endbr64
	movl	%edi, %eax
	andl	%esi, %edi
	xorl	%esi, %eax
	shrl	%eax
	addl	%edi, %eax
	ret
	.cfi_endproc

Linux 核心原始程式碼

指數移動平均（英語：exponential moving average，EMA或EWMA）是以指數式遞減加權的移動平均。 各數值的加權影響力隨時間而指數式遞減，越近期的數據加權影響力越重，但較舊的數據也給予一定的加權值。 -維基百科

在 include/linux/average.h 定義了一個 macro，其參數包含：

• name（用於定義 struct 與 helper functions 名稱)
• _precision（用於定義固定精度值的小數位數)
• _weight_rcp（權重，用於確定新資料與舊資料的加權方式)

首先，定義一個開頭為 ewma_ 的結構，內部包含一個 unsigned long internal。

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

使用前的初始化，檢查各項參數，並令 e->internal = 0。這能使其在 compile 時期就作檢查。

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

e->internal >> (_precision) 也就是讀取 internal 的整數值。

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

此函數為輸入 unsigned long val 並計算 e->internal 的值。
其中 weight_rcp 為

$lo{g}_2$(_weight_rcp)，因為 EWMA 的權值是利用指只函數計算。

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

根據公式 internal = 1/weight_rcp * val + (1 - 1/weight_rcp) * internal 用位移的方式去實做。

internal = (((internal << weight_rcp) - internal) + (val << precision)) >> weight_rcp

測驗 2 不需要分支的 Max

程式運作原理

分析不需要分支的設計的 Min

找到兩個有號整數的最校值

int32_t min(int32_t a, int32_t b) {
    int32_t diff = (a - b);
    return b + (diff & (diff >> 31));
}

1. 計算差值 diff

nt32_t diff = (a - b);

2. 找到差值得 sign bit

(diff >> 31)

3. 判斷 sign bit 正負
   • 正（a > b）：b 即為最小值。
   • 負（b > a）：b 要加上 diff，才會變成 a。

b + (diff & (diff >> 31))

算術位移 (Arithmetic shift) : 補上號數 (sign bit) 也就是最高有效位元的值在左側

兩個無號整數的最大值

找到兩個無號整數的最大值，達到 branchless 的實做。

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

首先，我們可以根據 a ^ ((EXP4) & -(EXP5)) 得知當 a > b 時 ((EXP4) & -(EXP5)) 應該要是 0。
反之，當 a < b 時，((EXP4) & -(EXP5)) 應該要是 a ^ b。

• a >= b ⇒ a ^ 0
• a < b ⇒ a ^ (a ^ b)

題目提示 EXP4 為 bitwise 操作，所以 EXP4 應該是 a ^ b。

a ^ ((a ^ b) & -(EXP5))	// EXP4 = a ^ b

題目提示 EXP5 為 logical operator 操作，代表結果只可能為 0 或 1。
然後我們可以看到他是先作負號運算

• -0 = 0x00000000
• -1 = 0xFFFFFFFF

接著與 a ^ b 作 AND 運算，因此當 a >= b 時，EXP5 要等於 0x00000000，反之要等於 0xFFFFFFFF。將此寫法在精簡，故 EXP5 = a < b。

a ^ ((a ^ b) & -(a > b))	// EXP5 = a > b

32 位元有號整數的最大值

int32_t max(int32_t a, int32_t b)
{
    return a ^ ((a ^ b) & -(a < b));
}

實做後發現，因為 -(a < b) 不論 sign 或 unsign 其結果都是 0x00000000 或 0xFFFFFFFF 其中之一，故可以用相同方法實做。

找出Linux 核心原始程式碼中更多類似的實作手法

測驗 3 64位元 GCD

程式運作原理

GCD 演算法運作原理

總是拿比較大的數值取餘數(相減)比較小的數值，直到其中一個數值為 0 時，另一個數值就是最大公因數。

1. 如果 u, v 都是 0 的話，其最大公因數必為 0。

• $gcd(0,0)=0$

2. 如果 u, v 其中一個為 0 的話，其最大公因數必為不為 0 的那個數。

• $gcd(u,0)=u$
• $gcd(0,v)=v$

3. 接著將 v 當作除數; u 當作被除數，作 u % v。
4. 令餘數(u % v)為除數v，令原本的除數 t 為 u。
5. 重複 2, 3 直到 v = 0，u 即為最大公因數。

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

改寫等價實做

這裡將原本的 Modulo Operator，改為用其他 Operator 實作。

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

先計算兩個數字的 2 指數最大公因數。

• if u, v are both even,
  $gcd(u,v)=2\ast gcd(\frac{u}{2},\frac{v}{2})$

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

因為我們已經找到公因數

$2^{shift}$ ，後續不會在遇到公因數 2，因此可以先將 2 除掉。
!(u & 1) 是用於判斷是否為偶數。

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

進入迴圈後，同樣先將 v 變為奇數。

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

根據輾轉相除法的演算法，這裡將取餘數改成用連續相減來達成。

• u, v 相減直到 u > v 時交換(此時 v 為餘數)。
• 迴圈條件應該要是餘數為 0 時候跳出，故 COND = v。

do {
    while (!(v & 1))
    	v /= 2;
    if (u < v) {
        v -= u;
    } else {
        uint64_t t = u - v;
        u = v;
        v = t;
    }
} while (COND);

而最後回傳的值要在乘上先前除掉的

$2^{shift}$，所以 RET = u << shift

return RET;

透過 __builtin_ctz 改寫 GCD

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.

回傳 interger x 最末端 0 的個數。
因此，我們可以從上述程式法發現__builtin_ctz 可以用於偶數除法的部份。

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

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

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

等價實作

uint64_t gcd64(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v;
    int ctz_u = __builtin_ctz(u);
    int ctz_v = __builtin_ctz(v);
    int shift = ctz_u > ctz_v ? ctz_v : ctz_u;
    u >> shift;
    v >> shift;

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

分析效能差異

Linux 核心的 GCD 實作

測驗 4 bit array

程式運作原理

此函式的功能是將一個長度為 bitmapsize 的 64 位元無號數陣列的所有 1 記錄下來，並用陣列 out 將這些 1 依序記錄下來，最後回傳 1 的數量。

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

為了找到陣列中的每一個 1，它利用 for 迴圈一個一個慢慢找，然而我們可以用 __builtin_ctzll 來找到最低位元 1 的位置來改寫。

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

首先，我們看到 t 會在後面與 bitset作 XOR 運算，代表其作用可能是清掉 bitset 的 1。因此，我們可以得知 EXP6 應該是取得最低位元 1 的位置，故 EXP6 = bitset & -bitset。

測驗 5 Fraction to Recurring Decimal

leetcode 166. Fraction to Recurring Decimal

程式運作原理

程式規則

計算一個分數的十進制成果。

• 若相除後有循環小數，則需 enclose 重複的部份

Input: numerator = 4, denominator = 333
Output: "0.(012)"

用 linux list.h 實作

使用 linux kernel 的 list API 來實做 hash table。

struct rem_node {
    int key;
    int index;
    struct list_head link;
};  

find 這個 function 會根據 key 找到 hash table 中對應的 key 是否存在，若不存在則回傳 -1。

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

若分子或分母任意一個為 0 則回傳 0。

if (denominator == 0) {
	result[0] = '\0';
    return result;
}

if (numerator == 0) {
    result[0] = '0';
    result[1] = '\0';
    return result;
}

改用 long long 型態來處理，避免 overflow。

/* 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++ = '-';

這裡處理除法以及取餘數的部份，若不是整除 (remainder !＝0)，則代表要處理小數部份，反之可以直接回傳相除的成果。

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++ = '.';

接來就是確認餘數（remainder）是否為循環小數

• 先從 hash table 找尋 remainder 是否存在
• 若不存在（pos = -1），則將 remainder 加到 hash table 中，並紀錄 i(位數)。
  • MMM = list_add
  • EEE = &heads[remainder % size]
• 若存在（pos >= 0），代表有循環小數，往回尋找重複的位數，加入括號
  • PPP = pos--
• 計算 (remainder * 10) / d，並將商加到 q 中
• 令 remainder 為餘數，重複上述動作直到 remainder 為 0

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

測驗 6 __alignof__

程式運作原理

alignof 功能

__alignof__ 的功能是能夠指定物件在記憶體中的對齊，將資料儲存在資料大小倍數的位址時，可以改善 cache 的效能。

等價實做

我們可以利用 struct 的對齊性質來實作 __alignof__。首先，可以看到 struct 內部包含一個 char c 以及一個 t _h，目的是找到型別 t 對齊時所需的位元大小，因此M = _h, X = 0 計算頭尾相減即為答案。

/*
 * 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)

測驗 7 FizzBuzz