---
tags: linux2022
---

# 2022q2 Homework (fibdrv)

contributed by <`ray90514`>

## 實驗環境

$ gcc --version
gcc (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
$ lscpu
Architecture:                    x86_64
CPU op-mode(s):                  32-bit, 64-bit
Byte Order:                      Little Endian
Address sizes:                   43 bits physical, 48 bits virtual
CPU(s):                          8
On-line CPU(s) list:             0-7
Thread(s) per core:              2
Core(s) per socket:              4
Socket(s):                       1
NUMA node(s):                    1
Vendor ID:                       AuthenticAMD
CPU family:                      23
Model:                           24
Model name:                      AMD Ryzen 5 PRO 3500U w/ Radeon Vega Mobile Gfx
Stepping:                        1
Frequency boost:                 enabled
CPU MHz:                         1675.446
CPU max MHz:                     2100.0000
CPU min MHz:                     1400.0000
BogoMIPS:                        4192.05
Virtualization:                  AMD-V
L1d cache:                       128 KiB
L1i cache:                       256 KiB
L2 cache:                        2 MiB
L3 cache:                        4 MiB
NUMA node0 CPU(s):               0-7

### 排除干擾效能分析的因素

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

執行的時候使用 `taskset -c 1`

### 用統計手法去除極端值

每一組測試總共執行 50 次，去掉最大值和最小值後剩下的取平均。

## 大數運算

使用 `unsigned long long` 的動態陣列儲存大數，並且有紀錄陣列長度跟最大長度的變數。

struct BigN {
    int len;
    int max_len;
    unsigned long long *digits;
};

將大數傳遞給 user 的時候只要根據 `len` 傳遞對應大小的 `digits` 。

### init

static struct BigN *init_BigN(int digits_num)
{
    struct BigN *n = kmalloc(sizeof(struct BigN), GFP_KERNEL);
    if (!n)
        return NULL;
    
    n->digits = kmalloc(digits_num * sizeof(unsigned long long), GFP_KERNEL);
    if (!n->digits) {
        kfree(n);
        return NULL;
    }

    n->len = 1;
    n->max_len = digits_num;    
    for (int i = 0; i < digits_num; i++)
        n->digits[i] = 0ULL;
    return n;
}

對 `BigN` 初始化， `digits_num` 是陣列的大小。

### free

static void free_BigN(struct BigN *n)
{
    kfree(n->digits);
    kfree(n);
}

釋放 `BigN` 的空間。

### swap

static inline void swap_BigN(struct BigN *dst, struct BigN *src)
{
    struct BigN temp = *dst;
    *dst = *src;
    *src = temp;
}

參考 [bignum](https://github.com/0xff07/bignum/blob/fibdrv/bignum.c) ，可以省下複製大數的時間。

### add

static void add_BigN(struct BigN *output,
                     const struct BigN *x,
                     const struct BigN *y)
{
    int max_len = max(x->len, y->len);
    unsigned long long carry = 0;

    for (int i = 0; i < max_len; i++) {
        unsigned long long result = x->digits[i] + carry;
        carry = x->digits[i] > ~carry || result > ~y->digits[i];
        result += y->digits[i];
        output->digits[i] = result;
    }

    if (output->len < output->max_len) {
        output->digits[max_len] = carry;
        output->len = max_len + carry;
    }
}

兩個大數相加，將結果存進 `output` 。具體作法是由低位將兩數相加，如果有 overflow 就將 `carry` 加至下一位的結果中。 `output` 可以與要做運算的兩數是相同的位址。
這裡有個細節是因為還有 `carry` 要補上去，可能造成額外的 overflow ，所以要做兩次判斷。

### sub

static void sub_BigN(struct BigN *output,
                     const struct BigN *x,
                     const struct BigN *y)
{
    unsigned long long borrow = 0;

    for (int i = 0; i < x->len; i++) {
        unsigned long long result = x->digits[i] - borrow;
        borrow = borrow > x->digits[i] || result < y->digits[i];
        result -= y->digits[i];
        output->digits[i] = result;
    }

    output->len = x->len;
    while (output->len > 1 && output->digits[output->len - 1] == 0)
        output->len--;
}

兩個大數相減，將結果存進 `output` ， `output` 可以與要計算的兩數位址相同。具體作法是從低位開始將兩數相減，如果有需要借位的情況，就將下一位的結果減一。這裡跟加法一樣，前一位的借位可能造成額外的借位，所以要判斷兩次，最後調整 `output` 的長度至第一個不為 0 的位數。有個要注意的地方是 `x` 的值必須大於等於 `y` 。

### lshift

static void lshift_BigN(struct BigN *output)
{
    unsigned long long right = 0;
    for (int i = 0; i < output->len; i++) {
        unsigned long long left = (output->digits[i] & ~(~0ULL >> 1)) > 0;
        output->digits[i] = output->digits[i] << 1 | right;
        right = left;
    }

    if (output->len < output->max_len) {
        output->digits[output->len] = right;
        output->len += right;
    }
}

將大數向左位移一位，相當於乘以二，不保證超出 `max_len` 的結果。具體作法是對陣列每個數做位移，將最高位的值放到下一個數的最低位。

### mul

static void mul_BigN(struct BigN *output,
                     const struct BigN *x,
                     const struct BigN *y,
                     struct BigN *carry)
{
    output->len = x->len + y->len;
    carry->len = output->len;

    if (output->len > output->max_len) {
        output->digits[0] = 0;
        output->len = 0;
        return;
    }

    for (int i = 0; i < output->len; i++) {
        output->digits[i] = 0;
        carry->digits[i] = 0;
    }

    for (int i = 0; i < x->len; i++) {
        for (int j = 0; j < y->len; j++) {
            unsigned __int128 result =
                (unsigned __int128) x->digits[i] * y->digits[j];

            /* add the lower of result */
            carry->digits[i + j + 1] +=
                output->digits[i + j] > (unsigned long long) (~result & ~0ULL);
            output->digits[i + j] += result & ~0ULL;

            /* add the upper of result */
            result >>= 64;
            if (i + j + 2 < output->max_len)
                carry->digits[i + j + 2] +=
                    output->digits[i + j + 1] > (unsigned long long) ~result;
            output->digits[i + j + 1] += result;
        }
    }
    add_BigN(output, output, carry);
    output->len -= output->digits[output->len - 1] == 0;
}

兩數相乘，將結果存入 `output` ，這裡 `output` 不能與兩數的位址相同。具體作法使用最慢但比較直觀的直式乘法，`x[i] * y[j]` 的結果累加至 `output[i + j]` ，為了避免補 `carry` 時要一直判斷有沒有 overflow ，將 overflow 的 `carry` 存至別的數，最後再一次加起來。

### add_constant

static void add_constant_BigN(struct BigN *output, unsigned long long c)
{
    unsigned long long carry = output->digits[0] > ~c;
    output->digits[0] += c;

    for (int i = 1; carry && i < output->len; i++) {
        carry = output->digits[i] > ~1ULL;
        output->digits[i]++;
    }

    if (output->len < output->max_len) {
        output->digits[output->len] = carry;
        output->len += carry;
    }
}

加一個常數至大數裡，處理進位即可。

### sub_constant

static void sub_constant_BigN(struct BigN *output, unsigned long long c)
{
    unsigned long long borrow = c > output->digits[0];
    output->digits[0] -= c;

    for (int i = 1; borrow && i < output->len; i++) {
        borrow = 1ULL > output->digits[i];
        output->digits[i]--;
    }

    while (output->len > 1 && output->digits[output->len - 1] == 0)
        output->len--;
}

將大數減一個常數，處理借位即可。

### 轉換成十進位

user 取得大數後，因為 `printf` 無法直接輸出，我的想法是先將陣列轉成方便輸出的形式，原本儲存的值改儲存`unsigned long long` 範圍內最大的 $10$ 的冪次方也就是 $10^{19}$ ，這樣就可以直接輸出。

#define MAX_P10 10000000000000000000ULL
void print_fib(int n, unsigned long long *digits, int len)
{
    unsigned __int128 value = 0;

    printf("Reading from " FIB_DEV " at offset %d, returned the sequence ", n);
    for (int i = 0; i < len; i++) {
        value = 0;
        digits[len] = 0;
        for (int j = len - 1; j >= i; j--) {
            value = value << 64 | digits[j];
            digits[j + 1] = value / MAX_P10;
            value = value % MAX_10P;
        }
        digits[i] = value;
        len += digits[len] > 0;
    }
    printf("%llu", digits[len - 1]);
    for (int i = len - 2; i >= 0; i--) {
        printf("%019llu", digits[i]);
    }
    printf("\n");
}

中間那段程式碼在做進位轉換，先是對大數除以 $10^{19}$ ，第 n 次除法的餘數為轉換後第 n 位的值，然後商成為下一輪的被除數，直到商為 0 。對大數的除法是採用長除法，在這段程式碼，餘數直接放在當前的位置，商放在高一位的位置，這樣不用額外的空間。

以上這些大數的運算比較針對 Fibonacci Number 的計算，所以會有些限制，但都符合之後計算的需求。

## Fibonacci Number

### 所需空間

Fibonnaci Number 有一個近似值的公式 $0.4472135955 * 1.61803398875^n$
利用這個公式我們可以近似出所需的空間。

- 二進位用 `unsigne long long` 儲存 : `2 + 7 * n / 640`
- 十進位用 `unsigne long long` 儲存 : `2 + 21 * n / 1900`

為了運算方便有多加一，如果要更準確一點的話使用以下的值 $log(F(n)) = 0.20898764025 * n, log_2(F(n)) = 0.694241914 * n$

### 迴圈

比較直觀的作法是根據定義 `F(n) = F(n - 1) + F(n - 2), F(0) = 0, F(1) = 1` 寫出一個迴圈的版本。

static struct BigN *fib_sequence_iterative(long long k)
{
    int digits_num = 2 + k * 7 / 640;
    struct BigN *f_n_prev = init_BigN(digits_num);
    struct BigN *f_n = init_BigN(digits_num);

    if(!(f_n && f_n_next))
        return NULL;

    f_n_prev->digits[0] = 0ULL;
    f_n->digits[0] = 1ULL;
    
    if(k == 0)
        f_n->digits[0] = 0ULL;
    for (long long i = 2; i <= k; i++) {
        add_BigN(f_n_prev, f_n_prev, f_n);
        swap_BigN(f_n_prev, f_n);
    }

    free_BigN(f_n_prev);
    return f_n;
}

因為 `add_BigN` 可以輸出至與輸入相同的大數，將結果儲存在算 F(n + 1) 時不需要用到的原本的 `f_n_prev` ，最後將已經是 F(n) 的 `f_n_prev` 跟 F(n - 1) 的 `f_n` 交換，這樣就不需要額外的空間。

### fast doubling

利用 F(n) 與 F(n + 1) 可以算出 F(2n) 與 F(2n + 1) 的值，公式如下。
$F(2n) = F(n) * [2F(n + 1) - F(N)]$
$F(2n + 1) = F(n + 1)^2 + F(N)^2$

static struct BigN *fib_sequence_test(long long k)
{
    int digits_num = 2 + k * 7 / 640;
    struct BigN *a = init_BigN(digits_num);
    struct BigN *b = init_BigN(digits_num);
    struct BigN *aa = init_BigN(digits_num);
    struct BigN *bb = init_BigN(digits_num);
    struct BigN *carry = init_BigN(digits_num);
    unsigned long long i = 1ULL << fls(k) >> 1;

    if (!(a && b && aa && bb && carry))
        return NULL;
    a->digits[0] = 0ULL;
    b->digits[0] = 1ULL;

    while (i) {
        mul_BigN(aa, a, a, carry);
        mul_BigN(bb, a, b, carry);
        lshift_BigN(bb);
        sub_BigN(a, bb, aa);

        mul_BigN(bb, b, b, carry);
        add_BigN(b, bb, aa);
        if (k & i) {
            add_BigN(a, a, b);  // m++
            swap_BigN(a, b);
        }
        i >>= 1;
    }

    free_BigN(b);
    free_BigN(aa);
    free_BigN(bb);
    free_BigN(carry);
    return a;
}

`a` 為 F(n) ， `b` 為 F(n + 1)

### fast doubling 與迴圈版本的比較

每隔 60 測試一次，總共測試 100 次
由圖可以很明顯看出迴圈的執行時間遠大於 fast doubling ，接下來會對 fast doubling 進行修改。

### 改進 1

在之前的程式碼中每一輪都會算出 F(n) 和 F(n + 1) ，就連最後一輪也是，但實際上只需要 F(n) 或是 F(n + 1) 就好，因此對最後一輪做特別的處理。

while (i > 1) {
    mul_BigN(aa, a, a, carry);
    mul_BigN(bb, a, b, carry);
    lshift_BigN(bb);
    sub_BigN(a, bb, aa);

    mul_BigN(bb, b, b, carry);
    add_BigN(b, bb, aa);
    if (k & i) {
        add_BigN(a, a, b);  // m++
        swap_BigN(a, b);
    }
    i >>= 1;
}
/* last round */
if (k & i) {
    mul_BigN(aa, a, a, carry);
    mul_BigN(bb, b, b, carry);
    add_BigN(a, aa, bb);
} else {
    lshift_BigN(b);
    sub_BigN(b, b, a);
    mul_BigN(aa, b, a, carry);
    swap_BigN(aa, a);
}

接下來比較兩者的差異，每隔 99 測試一次，總共測試 100 次。


從上面的程式碼可以看出 k 為奇數的時候可以省下一次的乘法運算，為偶數時可以省下兩次乘法運算，省下的乘法運算也是計算中最耗時的最後一輪，因此對比原本的寫法可以改進約 20% 和 43% 。奇數與偶數差了一次乘法運算造成了時間呈鋸齒狀。

### 改進 2

參考自 [斐波那契数列第一千万项怎么求](https://zhuanlan.zhihu.com/p/98064307) ，使用 [Catalan's identity](https://en.wikipedia.org/wiki/Cassini_and_Catalan_identities) 改寫原本的公式
$F(n) * F(n + 1) = F(n + 1)^2 - F(n)^2 - (-1)^n$
$F(2n) = 2[F(n + 1)^2 - F(n)^2] - F(n)^2 -2(-1)^n$
這樣就可以省下一次乘法，程式碼如下

while (i > 1) {
    mul_BigN(aa, a, a, carry);
    mul_BigN(bb, b, b, carry);
    sub_BigN(a, bb, aa);
    lshift_BigN(a);
    (k & (i << 1)) ? add_constant_BigN(a, 2) : sub_constant_BigN(a, 2);
    sub_BigN(a, a, aa);
    add_BigN(b, aa, bb);
    if (k & i) {
        add_BigN(a, a, b);  // m++
        swap_BigN(a, b);
    }
    i >>= 1;
}
/* last round */
if (k & i) {
    mul_BigN(aa, a, a, carry);
    mul_BigN(bb, b, b, carry);
    add_BigN(a, aa, bb);
} else {
    lshift_BigN(b);
    sub_BigN(b, b, a);
    mul_BigN(aa, b, a, carry);
    swap_BigN(aa, a);
}


對比了前一個改進的寫法加速約 20% 。使用 `perf record` 查看結果如下

  97.43%  test     [kernel.kallsyms]      [k] mul_BigN.isra.0
   0.70%  test     [kernel.kallsyms]      [k] add_BigN.isra.0
   0.35%  test     [kernel.kallsyms]      [k] sub_BigN.isra.0
   0.25%  test     [kernel.kallsyms]      [k] init_BigN
   0.14%  test     [kernel.kallsyms]      [k] lshift_BigN

很明顯的運算的瓶頸在乘法上面，若要再改進的話要從 `mul_BigN` 下手。

### 改進 3

參考 [KYG-yaya573142 的報告](https://hackmd.io/@KYWeng/rkGdultSU) ，使用 inline asm 直接使用乘法結果的高位和低位。與前次改進比較結果如下圖，加速了約 10% 。

for (int i = 0; i < x->len; i++) {
    for (int j = 0; j < y->len; j++) {
        unsigned long long high, low;
        __asm__("mulq %3" : "=a"(low), "=d"(high) : "%0"(x->digits[i]), "rm"(y->digits[j]));
        /* add the lower of result */
        carry->digits[i + j + 1] += output->digits[i + j] > ~low;
        output->digits[i + j] += low;

        /* add the upper of result */
        if (i + j + 2 < output->max_len)
            carry->digits[i + j + 2] += output->digits[i + j + 1] > ~high;
        output->digits[i + j + 1] += high;
    }
}


### 改進 4

實作 [Karatsuba algorithm](https://en.wikipedia.org/wiki/Karatsuba_algorithm) ，假設 AB 和 CD 相乘，乘積可以看成 AC BC + AD BD 三部分，中間的部分可以改寫為 (A + B)(C + D) - AC - BD ，這樣可以少一次 O(n / 2) 的乘法，除此之外 AC BD (A + B)(C + D) 也可以用同樣的方法算乘積，以下是遞迴至某個長度的兩數相乘就改用原本乘法的比較。

- threshold = 1


- threshold = 8


- threshold = 64


- threshold = 640


由以上圖可以看出 threshold 為 8 和 64 的效果較好，太大或太小的效果較差。