2022q1 Homework2 (quiz2)

contributed by < ChenBoSyun >

測驗一

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

解題思路: (a >> 1) + (b >> 1) 等同於 a / 2 + b / 2
因此推斷 EXP1 是為了加回當 a, b 都是奇數時，被向下取整捨去的 1 (0.5 + 0.5)

EXP1 = a & b & 1

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

解題思路: average 可以視為相加後再除以 2
我們先做出一個加法器 add， a & b 是進位的部份，因此需要左移 1，a ^ b 是沒有進位的部份，兩者相加後就是完整的加法。

uint32_t add(uint32_t a, uint32_t b)
{
    return ((a & b) << 1) + (a ^ b);
}

現在已經跟 average 的答案很相近了，只差了一個除以 2，也就是右移 1。

EXP2 = a & b
EXP3 = a ^ b

解析組合語言

透過 gcc -S -O3 add.c 生成最佳化後的組合語言

version1

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

average:
.LFB0:
	.cfi_startproc
	endbr64
	movl	%esi, %eax
	subl	%edi, %eax
	shrl	%eax
	addl	%edi, %eax
	ret
	.cfi_endproc

%edi, %esi, %eax 都是常見的暫存器，分別是第一個參數，第二個參數，返回值
%edi 的首字母 e 代表該暫存器存取 32bit

movl, subl 等指令的尾字母 l 代表指令操作對象的大小為 4 byte

指令	指令說明	執行完的結果
`movl %esi, %eax`	把參數2 high 賦予給回傳值	`%eax: high`
`subl %edi, %eax`	%eax 減去參數1 low	`%eax: high - low`
`shrl %eax`	%eax 做 1 bit 的邏輯右移	`%eax: (high - low) / 2`
`addl %edi, %eax`	%eax 加上參數1 low	`%eax: low + (high - low) / 2`

version2

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

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

指令	指令說明	執行完的結果
`movl %edi, %eax`	把參數1 a 賦予 %eax	`%eax: a`
`movl %esi, %edx`	把參數2 b 賦予 %edx	`%edx: b`
`xorl %esi, %edi`	參數1 a 跟參數2 b 做 XOR 並把結果放進參數2	`%edi: a ^ b`
`shrl %eax`	%eax 做 1 bit 的邏輯右移	`%eax: a >> 1`
`shrl %edx`	%edx 做 1 bit 的邏輯右移	`%edx: b >> 1`
`xorl $1, %edi`	%edi 跟 1 做 XOR 並把結果放進 %edi	`%edi: a ^ b ^ 1`
`addl %edx, %eax`	%edx 加進 %eax	`%eax: (a >> 1) + (b >> 1)`
`addl %edi, %eax`	%edi 加進 %eax	`%eax: (a >> 1) + (b >> 1) + (a ^ b ^ 1)`

version3

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

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

指令	指令說明	執行完的結果
`movl %edi, %eax`	把參數1 a 賦予 %eax	`%eax: a`
`andl %esi, %edi`	參數1 a 跟參數2 b 做 and 並把結果放進參數1	`%edi: a & b`
`xorl %esi, %eax`	參數2 b 跟 %eax 做 XOR 並把結果放進 %eax	`%eax: a ^ b`
`shrl %eax`	%eax 做 1 bit 的邏輯右移	`%eax: (a ^ b) >> 1`
`addl %edi, %eax`	%eax 跟參數 1 相加	`%eax: (a & b) + ((a ^ b) >> 1)`

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

測驗二

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

解題思路: 這題會用到 XOR 的一個特性，相同的數字做 XOR 會抵銷掉，例如 x ^ x = 0

EXP5 是條件判斷式，結果為 false 或 true
當 EXP5 結果為 false 時
EXP4 & -0 = 0
a ^ 0 = a
因此推斷 EXP5 = a < b

當 EXP5 = a < b 為 true 時
EXP4 & -1 = EXP4
a ^ EXP4 需要等於 b，因為 a < b 為 true，要回傳 b
推斷 EXP4 = a ^ b

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

找出 Linux 核心原始程式碼中，branchless / branch-free 的實作

測驗三

解釋下面這段 gcd 程式碼運作原理























#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 both x and y are 0, gcd is zero
$g c d (0, 0) = 0$

$g c d (x, 0) = x$ and
$g c d (0, y) = y$ because everything divides
$0$

If x and y are both even,
$g c d (x, y) = 2 * g c d (\frac{x}{2}, \frac{y}{2})$ because
$2$ is a common divisor.
Multiplication with
$2$ can be done with bitwise shift operator.

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

On similar lines, if x is odd and y is even, then
$g c d (x, y) = g c d (x, \frac{y}{2})$ It is because
$2$ is not a common divisor.

if (!u || !v) return u | v; 對應到演算法 1, 2 點，當 u, v 有任一個數字為 0 時，回傳 u | v
若 u == 0, v != 0 回傳 v，同理 u ！= 0, v == 0 回傳 u，當 u, v 皆為 0 時回傳 0
6 ~ 8 行對應到演算法第 3 點，提出公因數為 2 的倍數

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

9 ~ 10 行對應到演算法第 4 點


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

11 ~ 22 行，求非 2 的倍數的公因數，並乘上先前計算的 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;

算法核心是輾轉相減法，結束條件是當其中一方減至 0
COND = v > 0
結束時，剩餘的那方就是公因數，乘上之前計算的 2 倍數的公因數
RET = u << shift

以下用兩個實際例子說明演算法運作
表格的每列是每次迴圈執行完的結果

以 u = 3, v = 5 為例

u	v
3	5
3	2
3	1
1	2
1	1
1	0

回傳 u << shift 即 1 << 0，符合

g c d (3, 5) = 1

以 u = 21, v = 14 為例

u	v
21	14
21	7
7	14
7	7
7	0

回傳 u << shift 即 7 << 0，符合

g c d (21, 14) = 7

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

Built-in Function: 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.

__builtin_ctz 可以用來把

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

取代成

v = v >> __builtin_ctz(v);

改寫後的 gcd 程式碼

#define MIN(x, y) (((x) < (y)) ? (x) : (y))
uint64_t gcd64(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v;
    int shift;

    int ctz_u = __builtin_ctz(u);
    int ctz_v = __builtin_ctz(v);
    shift = MIN(ctz_u, ctz_v);
    u = u >> ctz_u;
    v = v >> ctz_v;
    
    do {
        v = v >> __builtin_ctz(v);
        if (u < v) {
            v -= u;
        } else {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while (v > 0);
    return u << shift;
}

執行 100 萬次，觀察兩者的執行時間

int main(){
    clock_t begin = clock();
    int i = 0;
    srand(time(0));
    for(i = 0; i < 1000000; i++){
        gcd64((uint64_t)rand(), (uint64_t)rand());
    }
    clock_t end = clock();
    double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    printf("gcd64 spend time for naive is %f second\n", time_spent);
}

gcd64 spend time for naive is 0.309016 second
gcd64 spend time for __builtin_ctz is 0.116132 second

__builtin_ctz 方式性能提昇 2.66 倍

另外嘗試用 objdump 觀察組合語言跟原始程式碼的對照

$gcc -o gcd -g gcd.c
$objdump -S gcd

紀錄一下使用 perf record 也可以看到組合語言跟原始程式碼的對照

$gcc -o gcd -g gcd.c
$perf record ./gcd
$perf report

原始版本的組合語言

do {
        while (!(v & 1))
    1240:	eb 0b                	jmp    124d <gcd64+0x84>
            v /= 2;
    1242:	48 8b 45 e0          	mov    -0x20(%rbp),%rax
    1246:	48 d1 e8             	shr    %rax
    1249:	48 89 45 e0          	mov    %rax,-0x20(%rbp)
        while (!(v & 1))
    124d:	48 8b 45 e0          	mov    -0x20(%rbp),%rax
    1251:	83 e0 01             	and    $0x1,%eax
    1254:	48 85 c0             	test   %rax,%rax
    1257:	74 e9                	je     1242 <gcd64+0x79>

__builtin_ctz 版本的組合語言

do {
        v = v >> __builtin_ctz(v);
    1229:	48 8b 45 d0          	mov    -0x30(%rbp),%rax
    122d:	f3 0f bc c0          	tzcnt  %eax,%eax
    1231:	89 c1                	mov    %eax,%ecx
    1233:	48 d3 6d d0          	shrq   %cl,-0x30(%rbp)

從組合語言可以看出
原始版本的 gcd 需要 7 個 instruction
__builtin_ctz 需要 3 個 instruction，雖然 tzcnt 需要 2 個 clock cycle，但總和 clock cycle 還是比較少

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

測驗四

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

上述程式碼，用來檢查 bitmap 中有哪些 bit 數值為 1，並且回傳他們在 butmap 中的位置
可以看到裡面的 for 迴圈，會逐一檢查 bit 是否為 1

若是使用 __builtin_ctzll 就不需要逐一檢查，它會回傳從尾端數過來會經過幾個 0，才會遇到第一個 1

使用 __builtin_ctzll 的程式碼
















#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 跟 12 行，推測這兩行的目的是把 rightmost 1 變成 0，其他部份保留原樣
再看第 12 行
bitset 想要保留非 rightmost 1
XOR 有一個特性是 x ^ 0 = x
從這兩點來看 EXP6 要做的事情就是保留 rightmost 1，把其他的部份都變成 0

EXP6 = bitset & -bitset

原理: 因為二補數的特性 -bitset = ~bitset + 1
將 bitset 拆成兩個部份觀察，rightmost 1 左邊的部份，與右邊的部份
左邊的部份不會因為加上 1 而改變
因此左邊的部份就是 x & ~x = 0，全變成 0
右邊的部份，因為原本都是 0，反轉之後加上 1，就變回跟原本一樣 x & x = x
因此 rightmost 1 會被保留，而其他 1 都變成 0

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

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

測驗五

測驗五的目的是給定一個分子跟分母，計算該分數的數值並用小數形式表示，最後以字串的形式回傳，該題目最關鍵的地方在於若是循環小數，需要用小括號包住循環的部份
例如: numerator = 7, denominator = 10
7/10 =1.42857142857…
回傳結果是 "1.(428541)"

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

程式碼分成兩部份，小數點前跟小數點後的部份
小數點前的部份是比較基本的數值運算`

小數點後的部份，需要判斷循環小數
為了紀錄何時開始循環，用一個 hash table 來紀錄小數點後的每個數字
當找到 hash table 中有重複的數字時，代表找到循環小數

小數點部份的程式碼
























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;

13行到20行是當 remainder 不在 hash table 中
因此17行要做的事情是把 node 加進 linked list

MMM 是把 node 加進 list 的 define macro
MMM = list_add

EEE 是要加進的 list_head 指標
EEE = &heads[remainder % size]

19, 20行做除法運算
注意這邊移動的是指標 q，並沒有移動指標 decimal

再看第四行 while (PPP > 0) 就很明顯了，當第三行 if (pos >= 0) 偵測到循環小數，先取出非循環的部份

PPP = (q - decimal)

改寫程式碼

原本的程式碼有一些地方可以修改的更好

處理負數的情況，可以用 bitwise 來實作 abs 來移除 branch

    /* deal with negtive cases */
    if (n < 0)
        n = -n;
    if (d < 0)
        d = -d;

abs with bitwise

    n = ((n >> 63) ^ n) - (n >> 63);
    d = ((d >> 63) ^ d) - (d >> 63);

sign 判斷也可以用 bitwise 來做

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

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

測驗六

/*
 * 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 macro 定義了結構 struct { char c; t _h; }
並將 0 轉型成該結構的指標，只要計算 t _h 成員與結構開頭的距離就可以知道 alignment 多少

M = _h
X = 0

查看 linux 原始碼實作Alignof(https://github.com/coreutils/gnulib/blob/master/lib/stdalign.in.h#L71)