2022q1 Homework2 (quiz2)

contributed by < Destiny0504 >

題目連結

測驗一

uint32_t average(uint32_t a, uint32_t b)
{
    // 判斷 a 和 b 都是奇數的情況下，才加一
    return (a >> 1) + (b >> 1) + (a & b & 1);
}

用 a & b 來實現 a - b 的操作，

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

兩種實作方式的組合語言比較

第一種實作方式的 x86 組合語言版本

average(unsigned int, unsigned int):
        // 將 edi 的值複製到 eax
        mov     eax, edi
        // 將 esi 的值複製到 edx
        mov     edx, esi
        // edi = edi & esi
        and     edi, esi
        // eax = edi >> 1
        shr     eax
        // edx = esi >> 1
        shr     edx
        // edi = edi & exi & 1
        and     edi, 1
        // eax = (esi >> 1) + (edi >> 1)
        add     eax, edx
        // eax = (esi >> 1) + (edi >> 1) + (edi & exi & 1)
        add     eax, edi
        ret

第二種實作方式的 x86 組合語言版本

average(unsigned int, unsigned int):
        // 將 edi 的值複製到 eax
        mov     eax, edi
        // edi = edi & esi
        and     edi, esi
        // eax = edi | esi
        or      eax, esi
        // eax = (edi | esi) >> 1
        shr     eax
        // eax = ((edi | esi) >> 1) + (edi & esi)
        add     eax, edi
        ret

	第一種實作方式	第二種實作方式
register 使用數量	4	3
add指令數	2	1
and指令數	2	1
mov指令數	2	1
shr指令數	2	1
or指令數	0	1

結論

第一種實作方式需要用到 a 、b 個別的值，所以需要花費額外的 register 將它存起來，所以在實作的要盡量避免這種所有要加起來的值都有關係的情況。

Exponentially weighted moving average ( EWMA )

優點

只要紀錄上一次的值就可以得到下一項應該是多少，能減少使用的記憶體空間。
可以動態的計算結果，不需要一次就看完所有的資料。

在 linux kernel 中的實作

linux kernel 中是以一個 pointer 指向一個 unsigned long 的變數，而這個變數就是用來儲存過去計算過得結果的。
考慮以下實作的方法，internal 先乘以 2 的 weight_rcp 次方減去原本的 internal，相當於
$internal \times \frac{1}{2^{w e i g h t_r c p}}$ ，在加上新的 val 。
- e.g. 如果 weight_rcp = 8 ，相當於每次把過去運算的結果乘以
  $\frac{7}{8}$ ，再加上新的 val 乘以
  $\frac{1}{8}$ 。

    WRITE_ONCE(e->internal, internal ?			\
        (((internal << weight_rcp) - internal) + 	\
            (val << precision)) >> weight_rcp :	\
        (val << precision));

測驗二

由 a <= b 來判斷我們應該回傳 a 還是 b

當 a > b 的時候 -(a <= b) 會得到零，所以可以直接無視 (b ^ a) 會得到的值，而零跟 a 做 XOR 則會得到 a ，所以回傳值為 a
當 a > b 的時候-(a <= b) 會得到 1 ，所以可以直接視為回傳的值為 a ^ (b ^ a)，而 a 跟 b ^ a 做 XOR 則會得到 b ，最後得到的回傳值為 b

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

正確答案是給 a ^ b 和 a < b, 但我想我的實作也可以達到同樣的效果。

Branch-free 版本

利用 a < b 時，a - b 會 overflow 的特性，所以只要判斷是否有 overflow 的情況即可。

因為 a 和 b 都是 unsigned int，所以不能只由最高位來判斷是否出現 overflow，還要接著判斷 a 的最高位是否大於
$2^{31} - 1$ ，如果兩者皆符合，才能判斷為遇到 overflow

uint32_t max(uint32_t a, uint32_t b)
{
    return a ^ ((b ^ a) & (-((a - b) >> 31) ^ -(a >> 31)));
}

測驗三

以下程式碼的運作原理為 Euclidean algorithm ，差別只在於每一步找都先把 v 換成奇數而已。而之所以可以只留下奇數的部份，是因為我們前面已經把 u 跟 v 所擁有的 2 的倍數的公因數都去除掉了。( 可以參考 reference 提到的演算法）

#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
    // 如果 u 跟 v 其中一個是 0，直接回傳不為 0 的那一個
    if (!u || !v) return u | v;
    int shift;
    // 用來找出 u 跟 v 公因數有 2 的幾次方
    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 (v);
    // 最後回傳的時候要把前面除去的 2 的次方的公因數乘回來
    return u << shift;
}
}

用 `__builtin_ctz` 改寫 GCD 的版本

uint64_t gcd64(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v;
    int shift = min(__builtin_ctz(u) , __builtin_ctz(v));

    u >>= shift, v >>= shift;
    u >>= __builtin_ctz(u);
    do {
        v >>= __builtin_ctz(v);
        if (u < v) {
            v -= u;
            printf("%ld %ld\n",u, v);
        } else {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    } while (v);
    return u << shift;
}

Reference

GCD 演算法

測驗四

因為 2 complement 的編碼特性，所以一個數字 a 在跟 -a 取 bitwise and 會剛好找出目前最低位元的 1 。
- e.g. a = 12 轉換為二進位 = 01100 , -a 轉換為二進位 = 10100，取 bitwise and 之後剛好等於 00100

#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 = bitset & -bitset;
            int r = __builtin_ctzll(bitset);
            out[pos++] = k * 64 + r;
            bitset ^= t;                           
        }
    }
    return pos;
}

測驗五

#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 (pos-- > 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;

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

        *q++ = (remainder * 10) / d + '0';
        remainder = (remainder * 10) % d;
    }

    strcpy(p, decimal);
    return result;
}

測驗六

我們在測驗中所作的版本能夠直接回傳 _h 的 datatype 在 memory 中所佔有的 byte 數量。
- 運作的方式是先拿到 _h 的記憶體位置，減去 char 的 pointer 的位置。
作法有點類似 container_of，差別在於一個是加上 offset 得到指向某一個變數的 pointer，而 alignof 是算出 offset
- 注意：這邊 char c 要擺在 t _h 前面，否則無法算出 t 的記憶體大小。（因為 struct 為連續的記憶體區塊的關係）

/*
 * 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)->_h) - (char *)0)

align function 說明

ALIGN(x, a) 要拿 x 以 a 作為對齊的基準進行向上對齊
- e.g. a = 128, x = 12 會回傳 128
- 注意：這個 function 只有 a 為 2 的整數次方算出來的對齊結果才是對的。因為這個 function 通常是用來計算記憶體的位置的，所以只有在上述的情況下算出來的結果才是對的也沒問題。

#define ALIGN(x, a)		__ALIGN_KERNEL((x), (a))

ALIGN_DOWN(x, a) 要拿 x 以 a 作為對齊的基準進行向下對齊
- e.g. a = 128, x = 12 會回傳 0
- 注意：這個 function 跟 ALIGN(x, a) 的使用條件一樣

#define ALIGN_DOWN(x, a)	__ALIGN_KERNEL((x) - ((a) - 1), (a))

Reference

linux 實作的 align
linux 實作的 align_kernel

#define __ALIGN_KERNEL_MASK(x, mask)	(((x) + (mask)) & ~(mask))
#define __ALIGN_KERNEL(x, a)		__ALIGN_KERNEL_MASK(x, (typeof(x))(a) - 1)

測驗七

因為在 compile 的時候 M3 、M5 都已經被計算好，而且乘法的運算速度又比除法快，所以會比另一種實作方法快。

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

        char fmt[9];
        strncpy(fmt, &"FizzBuzz%u"[(9 >> div5) >> (div3 << 2)], length);
        fmt[length] = '\0';

        printf(fmt, i);
        printf("\n");
    }
    return 0;
}

兩種方法的時間比較

用 linux 的 time 指令進行測量，且迴圈執行 100000 次
方法一為題目中的 naive.c ，方法二為上面的程式碼
從下方的結果來看，方法二確實有更好的 performance

	方法一	方法二
實驗一	0.438	0.330
實驗二	0.442	0.445
實驗三	0.478	0.447
實驗四	0.427	0.464
實驗五	0.457	0.478
實驗六	0.494	0.451
實驗七	0.518	0.303
實驗八	0.480	0.375
實驗九	0.394	0.447
實驗十	0.586	0.269
average	0.4714	0.4009
max	0.586	0.478
min	0.394	0.269

2022q1 Homework2 (quiz2)

測驗一

兩種實作方式的組合語言比較

結論

Exponentially weighted moving average ( EWMA )

優點

在 linux kernel 中的實作

測驗二

Branch-free 版本

測驗三

用 __builtin_ctz 改寫 GCD 的版本

Reference

測驗四

測驗五

測驗六

align function 說明

Reference

測驗七

兩種方法的時間比較

Read more

2022q1 Homework3 (quiz3)

2022q1 Homework3 (fibdrv)

2022q1 Homework4 (quiz04)

2022q1 Homework1 (quiz1)

用 `__builtin_ctz` 改寫 GCD 的版本