---
tags: linux2023
---

# 2023q1 Homework3 (fibdrv)
contributed by < `YSRossi` >

## 開發環境

$ gcc --version
gcc (Ubuntu 11.3.0-1ubuntu1~22.04) 11.3.0

$ lscpu
Architecture:            x86_64
  CPU op-mode(s):        32-bit, 64-bit
  Address sizes:         48 bits physical, 48 bits virtual
  Byte Order:            Little Endian
CPU(s):                  8
  On-line CPU(s) list:   0-7
Vendor ID:               AuthenticAMD
  Model name:            AMD FX-8320E Eight-Core Processor
    CPU family:          21
    Model:               2
    Thread(s) per core:  2
    Core(s) per socket:  4
    Socket(s):           1
    Stepping:            0
    Frequency boost:     enabled
    CPU max MHz:         3200.0000
    CPU min MHz:         1400.0000
    BogoMIPS:            6429.64

## 開發紀錄
### 改成 Fast Doubling 版本（成功計算到 F~92~）
參考作業說明所提供 Fast Doubling 的虛擬碼，搭配 bitwise 實作。

- `n` 的型態為 `long long`，使用 `sizeof(long long) << 3` 計算 `long long` 所佔有的位元數，減去 `n` 的 `leading 0-bits` 數量，得到 `n` 實際上有用到的位元數。其中 `n` 的 `leading 0-bits` 數量，使用 `__builtin_clzll(n)` 獲得。
- 在迴圈中，目前的位元以 `n & (1 << (i - 1))` 表示，將 1 左移 `i - 1` 把對應目前位元的位置設成 1 ，與 `n` 做 `and` 判斷目前位元是否為 1 。

Fast_Fib(n)
    a = 0; b = 1;       // m = 0
    for i = (number of binary digit in n) to 1
        t1 = a*(2*b - a);
        t2 = b^2 + a^2;
        a = t1; b = t2; // m *= 2
        if (current binary digit == 1)
            t1 = a + b; // m++
            a = b; b = t1;
    return a;

- 實驗比較 Fast Doubling 版本與原本方法在 kernel space 所花費的時間（還未排除干擾效能分析的因素）
:::warning
TODO : 嘗試解決在 GRUB 的設定檔中加入 ```isolcpus=7``` 重新開機後仍無法起作用的問題。

GRUB_CMDLINE_LINUX_DEFAULT="quiet splash isolcpus=7"

> 已解決，使用 ```sudo update-grub``` 更新 GRUB 設定檔
:::
- 實驗比較 Fast Doubling 與原始方法在 kernel space 所花費的時間(已排除干擾效能分析的因素)

    static long long fib_sequence_fast_doubling(long long k)
    {
        long long a = 0, b = 1;

        for (int i = 63; i >= 0; i--) {
            long long t1 = a * ((2 * b) - a);
            long long t2 = b * b + a * a;
            a = t1;
            b = t2;
            if ((k >> i) & 1) {
                t1 = a + b;
                a = b;
                b = t1;
            }
        }

        return a;
    }

- 若將上方程式碼的第 5 行改成下方程式碼，計算 fib(0) 在 user space 會花費大量時間。

    for (int i = 63 - __builtin_clzll(k); i >= 0; i--) {

`write` 回傳 Fast Doubling 版本的 kernel space time
`read` 回傳 Fast Doubling 版本計算的 Fibonacci 值

    /* client.c */
    for (int i = 0; i <= offset; i++) {
        struct timespec start, end;
        lseek(fd, i, SEEK_SET);
        clock_gettime(CLOCK_MONOTONIC, &start);
        sz = read(fd, buf, 1);
        clock_gettime(CLOCK_MONOTONIC, &end);
        printf("Reading from " FIB_DEV " at offset %d, returned the sequence ",
               i);
        printf("%lld.\n", sz);
        kernel_time = write(fd, write_buf, 0);
        user_time = (long long)(end.tv_sec * 1e9 + end.tv_nsec)
                 - (start.tv_sec * 1e9 + start.tv_nsec);
    }

- 在計算 fib(0) 會有異常的原因是 [__builtin_clzll](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html#Other-Builtins) 的參數不可以為0，因此在計算 fib(0) 時，直接回傳 0 即可解決問題。
> Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.


### 排除干擾效能分析的因素
* 限定 CPU 給特定的程式使用
使用 ```isolcpus``` Linux 核心起始參數，讓特定的 CPU 核在開機時就被保留下來，設定完成後使用 ```sudo update-grub``` 更新 grub 設定檔

    GRUB_CMDLINE_LINUX_DEFAULT="quiet splash isolcpus=7"

* 設定 [NO_HZ](https://www.kernel.org/doc/Documentation/timers/NO_HZ.txt) 減少影響程式執行的中斷，設定完成後使用 ```sudo update-grub``` 更新 grub 設定檔

    GRUB_CMDLINE_LINUX_DEFAULT="quiet splash isolcpus=7  nohz_full=7"

* 抑制 [address space layout randomization](https://en.wikipedia.org/wiki/Address_space_layout_randomization) (ASLR)

    $ sudo sh -c "echo 0 > /proc/sys/kernel/randomize_va_space"

* 設定 scaling_governor 為 performance。準備以下 shell script，檔名為 `performance.sh`:

    for i in /sys/devices/system/cpu/cpu*/cpufreq/scaling_governor
    do
        echo performance > ${i}
    done

之後再用 `$ sudo sh performance.sh` 執行

* 關閉 [Processor boosting](https://www.kernel.org/doc/Documentation/cpu-freq/boost.txt)

    sudo bash -c "echo 0 > /sys/devices/system/cpu/cpufreq/boost"

* 設定 [SMP IRQ Affinity](https://www.kernel.org/doc/html/next/core-api/irq/irq-affinity.html) ，減少目標 cpu 處理 interrupt request ，使用 cpu7 來做此實驗，所以 mask 為 0x7f

    sudo bash -c "echo 7f > /proc/irq/default_smp_affinity"

設定完成後，雖然還是有些浮動，但可以測量到相對穩定的實驗結果

### 支援大數運算（成功計算到 F~186~）
#### 大數資料結構
`uppper` 紀錄大數較高的 `64` 位元，`lower` 紀錄大數較低的 `64` 位元

struct BigN {
    unsigned long long upper, lower;
};

依照前面實作， Fast-doubling 需要加法、減法、乘法、左移一位運算，其中大數乘法運算還會使用右移一位運算。

##### `bn_add` 加法
兩數 `upper` 與 `lower` 分別相加，判斷 `lower` 相加後是否溢位，溢位的話將 `upper` 加一。

struct BigN bn_add(struct BigN a, struct BigN b)
{
    struct BigN result;
    result.lower = a.lower + b.lower;
    result.upper = a.upper + b.upper;
    if (a.lower > ~b.lower)
        result.upper++;

    return result;
}

##### `bn_sub` 減法
兩數 `upper` 與 `lower` 分別相減，若 `a.lower < b.lower`，代表需要跟 `upper` 借位，所以 `upper` 減一。

struct BigN bn_sub(struct BigN a, struct BigN b)
{
    struct BigN result;
    result.upper = a.upper - b.upper;
    if (a.lower < b.lower)
        result.upper--;
    result.lower = a.lower - b.lower;

    return result;
}

##### `bn_lshift` 左移一位 / `bn_rshift` 右移一位
實作左移一位，將兩數 `upper` 與 `lower` 分別左移一位，若 `lower` 在位移前的最高有效位元為 `1`，則 `upper` 的最低有效位元設為 `1`。實作右移一位與左移相似，兩數 `upper` 與 `lower` 分別右移一位，若 `upper` 在位移前的最低有效位元為 `1`，則 `lower` 的最高有效位元設為 `1`。

struct BigN bn_lshift(struct BigN a)
{
    a.upper <<= 1;
    if (a.lower & MSB_MASK)
        a.upper |= 1;
    a.lower <<= 1;
    
    return a;
}

struct BigN bn_rshift(struct BigN a)
{
    a.lower >>= 1;
    if (a.upper & 1)
        a.lower |= MSB_MASK;
    a.upper >>= 1;
    
    return a;
}

##### `bn_mul` 乘法
每次將 `a` 左移一位，並將 `b` 右移一位，找出 `b` 中值為 `1` 的位元，若找到為 `1` 的位元，則將目前 `a` 的值累加到 `result` ，達到類似乘法直式的運算方式。
例如 : a = 0b1010, b = 0b101
a * b = 0b1010 * 0b101 = 0b1010 + 0b101000 = 0b110010

struct BigN bn_mul(struct BigN a, struct BigN b)
{
    struct BigN result = bn_init(0, 0);
    while (b.upper || b.lower) {
        if (b.lower & 1)
            result = bn_add(result, a);
        a = bn_lshift(a);
        b = bn_rshift(b);
    }
    
    return result;
}

##### `bn_to_str` 將大數轉成字串
`input` 為要轉換成字串的大數， `num` 為最靠近且小於等於 `input` 的 `10` 的冪，將 `input` 不停減掉 `num` 直到 `input < num` ，減去的次數即目前 `input` 最高位的十進位數。
例如 : `input = 3021`, 最靠近且小於等於 `input` 的 `10` 的冪為 `num = 1000` ， `input` 需要減掉 `num` `3` 次才會滿足 `input < num` 的條件，因此得到最高位的十進位數為 `3` ，後面的位數以此類推。

int bn_to_str(char *buf, char *str)
{
    struct BigN input =
        bn_init(*(unsigned long long *) buf, *((unsigned long long *) buf + 1));
    struct BigN num = bn_init(0, 1);
    struct BigN ten = bn_init(0, 10);
    /* maximum value of BigN is 340,282,366,920,938,463,463,374,607,431,768,211,455 (39 digits)*/
    int pos = 39; /* highest possible digit index*/
    int dig_idx = 0;
    int dig_num;

    while (pos--) {
        for (dig_num = 0, num = bn_init(0, 1); dig_num < pos; dig_num++)
            num = bn_mul(num, ten);

        if (dig_idx || !bn_less(input, num) || !pos) {
            for (dig_num = 0; !bn_less(input, num); dig_num++)
                input = bn_sub(input, num);
            str[dig_idx++] = dig_num + '0';
        }
    }
    str[dig_idx] = 0;
    return dig_idx;
}

#### 實驗結果

### 支援大數運算(可動態改變大數長度)
#### 大數資料結構
`digits` 紀錄大數的數值，`size` 紀錄大數使用幾個 uint64_t，`alloc` 配置幾個 uint64_t 空間。

typedef u64 bn_digit;
typedef u64 bn_size;
struct {
    bn_digit *digits; 
    bn_size size;     
    bn_size alloc;    
}

##### 重新配置空間

若 `size` 超過 `alloc` 時，代表配置的空間不夠，需要 realloc，`BN_MIN_ALLOC(n, s)` 將大數 `n` 重新配置成 `s` 個 uint64_t 的大小，並確保配置的大小為 4 的倍數個 uint64_t。
`BN_SIZE(n, s)` 的功能類似，會額外將大數 `n` 的 `size` 設為 `s`。

#define BN_MIN_ALLOC(n, s)                                               \
    do {                                                                 \
        bn *const __n = (n);                                             \
        const bn_size __s = (s);                                        \
        if (__n->alloc < __s) {                                          \
            __n->digits =                                                \
                bn_resize(__n->digits, __n->alloc = ((__s + 3) & ~3U)); \
        }                                                                \
    } while (0)

#define BN_SIZE(n, s)                                                          \
    do {                                                                       \
        bn *const __n = (n);                                                   \
        __n->size = (s);                                                       \
        if (__n->alloc < __n->size) {                                          \
            __n->digits =                                                      \
                bn_resize(__n->digits, __n->alloc = ((__n->size + 3) & ~3U));  \
        }                                                                      \
    } while (0)

##### `bn_add` 加法

將大數 `a` 與大數 `b` 對應的位數相加，若有進位會累加到下一位。為了確保 a = c 與 b = c 時能正確運算不改到 `a` 或 `b` 的值，使用 `tmp` 紀錄結果，最後再將 `tmp` 複製到 `c`。

/* a + b = c */
void bn_add(bn *a, bn *b, bn *c)
{
    /* a + 0 or 0 + b */
    if (a->size == 0) {
        if (b->size == 0)
            c->size = 0;
        else
            bn_set(c, b);
        return;
    } else if (b->size == 0) {
        bn_set(c, a);
        return;
    }

    /* a + b = 2 * a */
    if (a == b) {
        bn_lshift(a, c);
        return;
    }

    bn_size min_size, max_size;
    bn_digit *tmp;
    bn *max_bn;
    if (a->size <= b->size) {
        min_size = a->size;
        max_size = b->size;
        max_bn = b;
    } else {
        min_size = b->size;
        max_size = a->size;  
        max_bn = a;
    }
    tmp = kmalloc((max_size + 1) * BN_DIGITS_SIZE, GFP_KERNEL);
    for (int i = 0; i < max_size + 1; i++)
        tmp[i] = 0;
    
    bn_digit cy = 0;
    for (int i = 0; i < min_size; i++) {
        bn_digit u = a->digits[i];
        bn_digit v = b->digits[i];
        cy = (u += cy) < cy;
        tmp[i] = u + v;
        cy += tmp[i] < v;
    }
    if (cy) {
        for (int i = min_size; i < max_size; i++) {
            bn_digit u = max_bn->digits[i];
            cy = (u += cy) < cy;
            tmp[i] = u;
        }
        if (cy)
            tmp[max_size++] = cy;
    } else
        bn_copy(max_bn->digits + min_size, max_size - min_size, tmp + min_size);
    BN_SIZE(c, max_size);
    bn_copy(tmp, max_size, c->digits);
    kfree(tmp);
}

##### `bn_lshift` 左移一位

從低位到高位尋訪大數的 `digits`，每次左移 1 位，藉由將 `digits` 右移 63 位獲取該 `digits` 的最高位，並與下一個 `digits` 做 OR 運算，將下一個 `digits` 的最低位設成目前 `digits` 的最高位。

/* c = a << 1 */
void bn_lshift(bn *a, bn *c)
{
    if (a->size == 0) {
        bn_set(c, a);
        return;
    }

    if(a != c)
        BN_SIZE(c, a->size);

    const unsigned int subp = 63;
    bn_digit cy = 0;
    bn_size size = a->size;
    for (int i = 0; i < size; i++) {
        bn_digit p = *(a->digits + i);
        *(c->digits + i) = (p << 1) | cy;
        cy = p >> subp;
    }

    if (cy) {
        BN_SIZE(c, c->size + 1);
        c->digits[c->size - 1] = cy;
    } 
}

#### `bn_mul` 乘法

大數 `a` 的每個位數與大數 `b` 的每個位數相乘，將結果紀錄在 `tmp` ，二個 64-bit 整數相乘會需要 128-bit ，所以使用 GCC Extension `__uint128_t`，透過右移 64 位，將較高的 64-bit 存在 `tmp` 目前位數的下一位，最後再將 `tmp` 複製到 `c`。

void bn_mul(bn *a, bn *b, bn *c)
{
    if (a->size == 0 || b->size == 0)
        c->size = 0;
    
    bn_size csize = a->size + b->size;
    bn_digit *tmp = kmalloc(csize * BN_DIGITS_SIZE);
    for (int i = 0; i < csize; i++)
        tmp[i] = 0;
   
    for (int i = 0; i < a->size; i++) {
        for (int j = 0; j < b->size; j++) {
            __int128 product = (__int128)a->digits[i] * (__int128)b->digits[j];
            bn_size low =  product;
            bn_size high =  product >> BN_DIGITS_BITS;
            tmp[i + j] += low;
            uint32_t cy = tmp[i + j] < low;
            tmp[i + j + 1] += high + cy;
            if (cy)
                cy = tmp[i + j + 1] <= high;
            else
                cy = tmp[i + j + 1] < high;

            for (int k = i + j + 2; cy; k++) {
                tmp[k] += cy;
                cy = tmp[k] < cy;
            }
        }
    }

    while (csize > 0 && tmp[csize - 1] == 0)
        csize--;
    BN_SIZE(c, csize);
    bn_copy(tmp, csize, c->digits);
    kfree(tmp);
}

##### `bn_to_str` 將大數轉成字串

基本想法為對大數除以 10 取餘數，反覆對其商除以 10 取餘數，將取得的餘數逐一串接起來，最後所得到的就是該數的十進位表示法。

123 / 10 = 12...3 <-   
 12 / 10 =  1...2 <-
  1 / 10 =  0...1 <-

- 方法一(對應到圖中的 Not divided by max_radix)
每次迴圈將大數除以 10 取餘數，這個餘數就是該數的十進位最低位的數值，保留商到下一次迴圈再除以 10，直到商等於 0 時停止迴圈。
方法一對大數一次只除以 10，需要跑很多次迴圈才可以讓商等於 0，因此有了方法二來加速。
- 方法二(對應到圖中 Divided by max_radix)
`max_radix` 為 10^19^，相當於 uint64_t 可表示值的範圍中，最大的 10 的冪的值，`max_power` 為 19，紀錄 10 的冪次。
每次迴圈藉由將大數除以 `max_radix`，獲得一個 uint64_t 可表示的餘數，再用方法一來取得最終的餘數，重複上述操作直到商等於 0 時停止迴圈。

static char *bn_to_str(const bn_digit *u,
                         bn_size size,
                         char *out)
{
    while (size && !u[size -1]) /* 取得 u 的有效 size(byte) */
        --size;
    /* u = 0 ~ 9 */
    if (size == 0 || (size == 1 && u[0] < 10)) {
        if (!out)
            out = malloc(2);
        out[0] = size ? radix_chars[u[0]] : '0';
        out[1] = '\0';
        return out;
    }

    const bn_digit max_radix = radix_table.max_radix;
    const unsigned int max_power = radix_table.max_power;

    if (!out)
        out = malloc(string_size(size) + 1);
    char *outp = out;

    bn_digit *tmp = BN_TMP_COPY(u, size);
    bn_size tsize = size;

    do {
        bn_digit r = bn_ddivi(tmp, tsize, max_radix);
        tsize -= (tmp[tsize - 1] == 0);
        unsigned int i = 0;
        do {
            bn_digit rq = r / 10;
            bn_digit rr = r % 10;
            *outp++ = radix_chars[rr];
            r = rq;
            if (tsize == 0 && r == 0)
                break;
        } while (++i < max_power);
    } while (tsize != 0);

    free(tmp);
    
    char *f = outp - 1;
    while (*f == '0')
        --f;
    f[1] = '\0';
    for (char *s = out; s < f; ++s, --f)
        SWAP(*s, *f);
    return out;
}


#### 降低 realloc 呼叫次數

在 [sysprog21/bignum](https://github.com/sysprog21/bignum) 中，初始化大數所 malloc 的大小皆相同，在計算 fib(k) 時，k 越大，所需的空間越多，realloc 的次數也越多。
透過預估 fib(k) 所需的位數，在初始化時就分配足夠的空間，降低 realloc 次數。
根據[作業說明](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-b)，當 $n$ 為足夠大的正整數時，
$$
F(n)\approx0.4472135955\times1.61803398875^n
$$

將近似值取 2 為底的對數，可得到需要使用的位元數
$$
log_2(0.4472135955\times1.61803398875^n)\approx-1.160964 + n\times0.694242
$$

除以 64 可得到需要使用的 `uint64` 數量
$$
(-1.160964 + n\times0.694242)\div64 = -0.0181400625 + n\times0.01084753125
$$

$$\lfloor\frac{-181400625 + n\times108475313}{10^{10}}\rfloor+1
$$

int bn_fib_digits_n(int n)
{
    if (n < 2)
        return 1;
    int digits = ((long)n * 108475313 - 181400625) / 10000000000 + 1;
    return digits;
}

由下圖能觀察到隨著數字增加，預先分配足夠空間的效果越明顯。

### 關於 [sysprog21/bignum](https://github.com/sysprog21/bignum) 的問題

:::info
在 [sysprog21/bignum/bignum.c](https://github.com/sysprog21/bignum/blob/master/bignum.c) 中，macro 使用 do while(0) 將多行內容框起來，我想原因是減少語法錯誤，像是下方在 if 中呼叫 macro，且沒有用大括號括住的情形依舊能正確執行，若 macro 沒有 do while(0) 則會出錯。
在 macro 中不用 if(1) { ... } 而是使用 do while(0) 將內容框起來的原因是單純 do while(0) 後可以加分號，因此使用該 macro 時可以較直觀的在後面加分號，還是有其他考量？

#define BN_MIN_ALLOC(n, s)        \
    do {                          \
        ...                       \
    } while (0) 

if (1)
    BN_MIN_ALLOC(n, s);

:::

:::info
在 [sysprog21/bignum/format.c](https://github.com/sysprog21/bignum/blob/master/format.c) 中，為什麼當 radix 不是 2 的冪時，結果需要額外加 2，而不是加 1？

static size_t apm_string_size(apm_size size, unsigned int radix)
{
    ASSERT(radix >= 2);
    ASSERT(radix <= 36);

    if ((radix & (radix - 1)) == 0) { 
        unsigned int lg = apm_digit_lsb_shift(radix);
        return ((size * APM_DIGIT_BITS + lg - 1) / lg) + 1;
    }
    /* round up to second next largest integer */
    return (size_t)(radix_sizes[radix] * (size * APM_DIGIT_SIZE)) + 2;
}

:::