---
tags: linux2022
---

# 2022q1 Homework3 (fibdrv)
contributed by < `ibat10clw` >
> [fibdrv](https://github.com/ibat10clw/fibdrv)

## 費氏數列
### 版本一 (DP)
為理解 client 和 fibdrv 之間的互動，我的第一個版本只是將 `f[i] = f[i-1] + f[i-2]` 以字元陣列儲存，並使用字元加法做改寫，同時根據作業說明中 clinet.c 的提示使用 `copy_to_user` 將 fibdrv 計算結果的字元陣列取出。

static void add(char *a, char *b, char *out)
{
    int len_a = strlen(a);
    int len_b = strlen(b);
    int i = 0;
    int carry = 0;
    int diff = len_a - len_b;
    for (i = len_b - 1; i >= 0; i--) {
        int sum = (a[i + diff] - '0') + (b[i] - '0') + carry;
        out[i + diff] = '0' + sum % 10;
        carry = sum / 10;
    }
    for (i = diff - 1; i >= 0; i--) {
        int sum = (a[i] - '0') + carry;
        out[i] = '0' + sum % 10;
        carry = sum / 10;
    }
    if (carry) {
        if (diff == 0)
            len_a++;
        for (i = len_a; i >= 0; i--) {
            out[i + 1] = out[i];
        }
        out[len_a + 1] = 0;
        out[0] = '0' + carry;
    }
    out[len_a] = 0;
}

加法部份是將傳入的字元陣列以較短的為準對齊，依序從最右邊相加到比較短陣列結束後，同時需要紀錄當前的相加是否需要進位，若有進位需求的話，則需要把結果累計到下一次的計算，再把剩下的數字填入 buf。

static long my_fibo(int k, char *buf)
{
    char f3[SIZE];
    char f2[SIZE] = "1";
    char f1[SIZE] = "0";
    for (int i = 2; i <= k; ++i) {
        add(f2, f1, f3);
        strncpy(f1, f2, strlen(f2) + 1);
        strncpy(f2, f3, strlen(f3) + 1);
    }
    if (k == 0) {
        copy_to_user(buf, f1, strlen(f1) + 1);
        return 0;
    }
    if (k == 1) {
        copy_to_user(buf, f2, strlen(f2) + 1);
        return 1;
    }
    copy_to_user(buf, f3, strlen(f3) + 1);
    return k;
}

關於費氏數列的計算是以 DP 的方式依序將 `f(n-1)` 與 `f(n-2)` 相加後把結果存到 `f3` ，之後分別把 `f2` 和 `f3` 複製給 `f1` 和 `f2`以便進行下一輪的計算，而 buf 是由 `bif_read` 中的參數傳入的，因此最後只要將結果放入 `copy_to_user` 即可在 client 中將其印出。
### 版本二 (Fast Doubling)

若要使用 fast doubling 的方法計算費氏數列還需要引入乘法與減法的操作

void sub(char *a, char *b, char *out)
{
    int len_a = strlen(a);
    int len_b = strlen(b);
    int i = 0;
    int borrow = 0;
    int idx = len_a - len_b;
    for (i = len_b - 1; i >= 0; i--) {
        int diff = (a[i + idx] - '0') - (b[i] - '0') - borrow;
        if (diff < 0) {
            diff += 10;
            borrow = 1;
        } else
            borrow = 0;
        out[i + idx] = '0' + diff;
    }
    for (i = idx - 1; i >= 0; i--) {
        int diff = (a[i] - '0') - borrow;
        if (diff < 0) {
            diff += 10;
            borrow = 1;
        } else
            borrow = 0;
        out[i] = '0' + diff;
    }
    out[len_a] = 0;
    i = 0;
    while (out[i] == '0')
        i++;
    if (i != 0) {
        for (int j = 0; j < len_a; ++j)
            out[j] = out[j + i];
    }
}

減法與加法的想法是類似的，不過從需要檢查進位變成需要檢查借位，同樣的以較短的輸入為準對齊後從右邊開始計算回來，每個回合將對應的位數相減後判斷是否需要借位，借助著直式減法的思路，當需要借位時可以將當前的差 +10 後把前一個位數 -1 ，透過這種方法便可以達成以字元鎮列為儲存單位的減法

最後的 while 迴圈是判斷是否有 leading zero ，若有的話就將有效位數全部向前挪動

void mul(char *a, char *b, char *out)
{
    int len_a = strlen(a);
    int len_b = strlen(b);
    int N = len_a + len_b;
    int i;
    memset(out, '0', N);
    out[N] = 0;
    for (i = len_b - 1; i >= 0; i--) {
        for (int j = len_a - 1; j >= 0; j--) {
            int digit1 = a[j] - '0';
            int digit2 = b[i] - '0';
            int pos = 1 + i + j;
            int carry = out[pos] - '0';
            int multiplication = digit1 * digit2 + carry;
            out[pos] = (multiplication % 10) + '0';
            out[pos - 1] += (multiplication / 10);
        }
    }

    i = 0;
    while (out[i] == '0')
        i++;
    if (i != 0) {
        for (int j = 0; j < N; ++j)
            out[j] = out[j + i];
    }
}

字元的乘法與前面的流程稍有不同，但一樣是模擬直式的計算流程取得結果
首先把兩字元陣列的長度相加，並初始化為 0 以儲存計算結果

接著兩層的 for 迴圈分別看過所有存在陣列內的數字，兩兩計算相乘的結果，同時每次計算時還得先找出要將結果填在 out 的哪個位置，若有進位的話則在將其填在 pos - 1 的地方，當所有的數字都看過後即可取得結果

最底下的 while 迴圈同樣是去除存在陣列中的 leading zero

static long fast_fib(int k, char *buf)
{
    if (k == 0) {
        copy_to_user(buf, "0", 2);
        return k;
    }

    if (k == 1) {
        copy_to_user(buf, "1", 2);
        return k;
    }
    int stack[100];
    int top = 0;
    while (k > 1) {
        stack[top++] = k;
        k /= 2;
    }
    char f1[SIZE] = "1";
    char f2[SIZE] = "1";
    char tmp1[SIZE];
    char tmp2[SIZE];
    char res1[SIZE];
    while (top--) {
        int n = stack[top];
        // f1 = f1 * (2 * f2 - f1)
        mul("2", f2, tmp1);
        sub(tmp1, f1, tmp2);
        mul(tmp2, f1, res1);
        // f2 = f1 ** 2 + f2 ** 2
        mul(f1, f1, tmp1);
        mul(f2, f2, tmp2);
        add(tmp2, tmp1, f2);
        if (n % 2 == 1) {
            add(f2, res1, tmp1);
            strncpy(f1, f2, strlen(f2) + 1);
            strncpy(f2, tmp1, strlen(tmp1) + 1);
        } else {
            strncpy(f1, res1, strlen(res1) + 1);
        }
    }
    copy_to_user(buf, f1, strlen(f1) + 1);
    return k;
}

在使用 fast doubling 的方法時需要根據當前的 fibonacci number 是偶數項或奇數項來決定下次的計算方法，因此這邊使用了一個 stack 來紀錄項數的奇偶，以便決定後續的計算方式。

當項數為奇數的時候，就用原始定義的 `fib(n) = fib(n-1) + fib(n-2)` 把數列往前推進一項，等待全部計算完畢後同樣使用 copy_to_user 將結果傳出。

### 效能評測
根據作業說明的提示，可以使用 `ktime` 相關的 API 在 kernel 中紀錄計算 fibonacci 的時間，再搭配還沒有被使用上的 `fib_write` 從 client 中獲得紀錄的時間

static long long fib_time_proxy(long long k, char *buf)
{
    kt = ktime_get();
    long long result = fast_fib(k, buf);
    kt = ktime_sub(ktime_get(), kt);

    return result;
}

在呼叫 fast_fib 的前後先計算時間差，之後在從 client 呼叫 write 得到結果

for (int i = 0; i <= offset; i++) {
    lseek(fd, i, SEEK_SET);
    clock_gettime(CLOCK_REALTIME, &t1);
    sz = read(fd, buf, sizeof(buf));
    clock_gettime(CLOCK_REALTIME, &t2);
    sz = write(fd, write_buf, strlen(write_buf));
    printf("%d %lld %ld %lld\n", i, sz, t2.tv_nsec - t1.tv_nsec, (t2.tv_nsec - t1.tv_nsec) - sz);//k, u, k2u
}

此外在 client 中同樣在 read 前後各紀錄一次時間，就可以計算執行時間
而 `copy_to_user` 的開銷則可以透過總執行時間扣掉在 kernel 的時間，空缺的那一塊就是 `copy_to_user` 的開銷

## 大數運算
上面的費氏數列實作雖然能計算出正確的答案，但由於使用字串的方式模擬 10 進位的運算，每次運算相當於只能處理 1 個位數，以現在 CPU 的規格來說並沒有完全發揮他的效能，因此先全面改寫運算的規則後再重新實作費氏數列。
### 用於儲存大數的資料結構

typedef struct _bn {
    unsigned int number[BNSIZE];
    unsigned int size;
} bn;

`unsigned int` 共可以儲存 32 bit 的資料，使用一個陣列儲存運算的結果，並且根據 size 來決定目前有效資料的長度，LSB 存於 index 0 的位置
示意圖如下:
若要解釋以上資料則從 index 0 開使將數字以 2 進位的方式連接，可得以下結果:
`00000000 00000000 000000001 11000011`
如此便可以將 2 進位的格式在輸出時轉換為 10 進位
### 將大數以十進位的方式輸出
這個部份參考了 [KYG-yaya573142](https://hackmd.io/@KYWeng/rkGdultSU#%E5%A4%A7%E6%95%B8%E9%81%8B%E7%AE%97) 的實作

static char *bn_to_string(const bn *src)
{
    // log10(x) = log2(x) / log2(10) ~= log2(x) / 3.322
    size_t len = (8 * sizeof(int) * src->size) / 3 + 2;
    char *s = kmalloc(len, GFP_KERNEL);
    char *p = s;
    memset(s, '0', len - 1);
    s[len - 1] = '\0';

    for (int i = src->size - 1; i >= 0; i--) {
        for (unsigned long d = 1U << 31; d; d >>= 1) {
            // double logical not can make the result to 1
            int carry = !!(d & src->number[i]);
            for (int j = len - 2; j >= 0; j--) {
                s[j] += s[j] - '0' + carry;
                carry = (s[j] > '9');
                if (carry)
                    s[j] -= 10;
            }
        }
    }
    // skip leading zero
    while (p[0] == '0' && p[1] != '\0') {
        p++;
    }
    memmove(s, p, strlen(p) + 1);
    return s;
}

一般將 2 進位的數字轉換為 10 進位時，可以借助下列公式
$\displaystyle\sum_{b=1}^{N} {2^N}\times b$ `for b = lsb to msb`

範例如下: `00000000 00000000 000000001 11000011`
根據 bit=1 所在的位置可以得出該 2 進位數的 10 進位表示為
$2^8+2^7+2^6+2^1+2^0$
由於輸入的 bn 可能非常的大，因此需要直接將其轉為字元陣列的形式保存，演算法描述如下
1. 從 bn 的 MSB 開始讀過所有的有效位元(根據size取捨)
* 對應程式碼第 10 行的 for 迴圈
3. 由於 bn 的每個欄位是 32 bit，因此使用 32 bit 的無號數對每個 bit 做 bitwise and 判斷是否為1，若出現的是 1 就將字元陣列 +1
* 對應程式碼 11 行的 for 迴圈
5. 從初始化字元陣列的尾巴開始把每個位數扣掉 ASCII 的偏移後與自己相加，相當於乘以 2 的效果，由於是從 MSB 開始做轉換，當第一個 bit 1 出現到轉換結束執行的次數恰好能將其恢復成 power of 2，同時也需要判斷字元陣列表示的十進位數是否需要進位
* 對應程式碼 14 行的 for 迴圈

如此便能夠在輸出時將 bn 轉化為十進位表示的字串，並且透過 `copy_to_user` 將結果傳到 userspace
### 大數加法
bn 的儲存格式在選擇時使用了 unsigned integer，unsigned integer 在 overflow 發生的時候會將結果 `mod` 該型態的最大值後再做保存，根據這個特性只需要判斷是否發生 overflow，若發生的話就將 bn 的下一個欄位 +1 來做進位，如此便能完成大數加法

#define MAXSLOT 0xffffffffUL+1UL
static void add_bn(bn *a, bn *b, bn *c)
{
    unsigned long carry = 0;
    int i;
    c->size = a->size;
    for (i = 0; i < b->size; ++i) {
        unsigned int sum = a->number[i] + b->number[i] + carry;
        carry =
            ((unsigned long) a->number[i] + b->number[i] + carry) >= MAXSLOT;
        c->number[i] = sum;
    }
    for (; i < a->size; ++i) {
        unsigned int sum = a->number[i] + carry;
        carry = ((unsigned long) a->number[i] + carry) >= MAXSLOT;
        c->number[i] = sum;
    }
    if (carry) {
        c->size++;
        c->number[i]++;
    }
}

其中 MAXSLOT 定義為 `UINT_MAX+1` 是為了 `branchless` 的實作，只要當前要相加的 bn 和 carry 加在一起後判斷是否會超過 `MAXSLOT` 若會的話則代表 overflow 發生了，這時就要將 carry 設為 1，這部份的操作則是透過 logical operator 來達成