---
# System prepended metadata

title: fibdrv 開發

---

[我的實作 github 連結](https://github.com/backink/fibdrv)
[問題描述及原始程式碼 — 2023 年 Linux 核心設計/實作課程作業 fibdrv](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-d)


## 費氏數列

費氏數列的原始定義為 
$$ 
F(n) = F(n - 1) + F(n - 2) \\
and\\
F(0) = 0, \ F(1) = 1 
$$
接著談及 [Fast Doubling](https://chunminchang.github.io/blog/post/calculating-fibonacci-numbers-by-fast-doubling) 手法，可以得到 
$$
F(2k) = F(k)[2F(k + 1) - F(k)] \\ 
F(2k + 1) = F(k + 1)^2 + F(k)^2
$$

運用 fast doubling 開發遞迴版本 C 語言程式碼如下
```c
static inline uint64_t fast_doubling(uint32_t target)
{
    if (target <= 2)
        return !!target;

    // fib(2n) = fib(n) * (2 * fib(n+1) − fib(n))
    // fib(2n+1) = fib(n) * fib(n) + fib(n+1) * fib(n+1)
    uint64_t n = fast_doubling(target >> 1);
    uint64_t n1 = fast_doubling((target >> 1) + 1);

    // check 2n or 2n+1
    if (target & 1)
        return n * n + n1 * n1;
    return n * ((n1 << 1) - n);
}
```

* iterative 方法的時間複雜度為 $O(n)$
* fast doubling 的時間複雜度降為 $O(\log n)$

雖然減少計算量，但遞迴方法仍然有重複計算的部份

### Bottom-up Fast Doubling

觀察數值的 2 進制表示，可發現該數是如何產生，以 $87_{10}$ 為例：
```text
                87 =  1010111            (87 >> i+1)
        i = 0 : 43 = (1010111 - 1) >> 1 = 101011
        i = 1 : 21 = ( 101011 - 1) >> 1 = 10101
        i = 2 : 10 = (  10101 - 1) >> 1 = 1010 
        i = 3 :  5 = (   1010 - 0) >> 1 = 101 
        i = 4 :  2 = (    101 - 1) >> 1 = 10 
        i = 5 :  1 = (     10 - 0) >> 1 = 1
        i = 6 :  0 = (      1 - 1) >> 1 = 0
                      	     ^
                         87 的第 i 個位元
```

若是進行移項並反過來看的話會變成：
```text
                   (87 >> i+1)
        i = 6 :  1 =        0 << 1 + 1 = 1
        i = 5 :  2 =        1 << 1 + 0 = 10
        i = 4 :  5 =       10 << 1 + 1 = 101
        i = 3 : 10 =      101 << 1 + 0 = 1010
        i = 2 : 21 =     1010 << 1 + 1 = 10101
        i = 1 : 43 =    10101 << 1 + 1 = 101011
        i = 0 : 87 =   101011 << 1 + 1 = 1010111
                                     ^
                              87 的第 i 個位元
```

從 $n=0$ 開始看，可發現每次位移後，只要檢查目標數對應的位元，即可知曉下次應以 $fib(2n)$ 還是 $fib(2n+1)$ 為基礎進行右移。

整理前述觀察，可知：
1. 從最高位元的 1 開始，此時 $n=1$，而：
    - $fib(n)=1$
    - $fib(n+1)=1$
2. 若下一個位元不存在的話跳到第 3 步，否則（假設目前為 $n=k$）：
    - 透過 $fib(k)$ 以及 $fib(k+1)$ 計算 $fib(2k)$ 以及 $fib(2k+1)$
    - 檢查下一個位元：
        - 0：$n = 2k$
        - 1：$n = 2k + 1$，此時需要 $fib(n+1)$ 讓下一迭代能夠計算 $fib(2n)$ 以及 $fib(2n+1)$
    - 回到第 2 步
3. 此時 $n$ 為目標數，回傳 $fib(n)$


對應的程式碼：

這裡運用 gcc builtin function 減少運算時間
```c
static inline uint64_t fast_doubling_iter(uint64_t target)
{
    if (target <= 2)
        return !!target;

    // find first 1
    uint8_t count = 63 - __builtin_clzll(target);
    uint64_t fib_n0 = 1, fib_n1 = 1;
    
    for (uint64_t i = count, fib_2n0, fib_2n1; i-- > 0;) {
        fib_2n0 = fib_n0 * ((fib_n1 << 1) - fib_n0);
        fib_2n1 = fib_n0 * fib_n0 + fib_n1 * fib_n1;

        if (target & (1UL << i)) {
            fib_n0 = fib_2n1;
            fib_n1 = fib_2n0 + fib_2n1;

        } else {
            fib_n0 = fib_2n0;
            fib_n1 = fib_2n1;
        }
    }
    return fib_n0;
}
```

### $F(92)$ 以後的數值會出現錯誤

在原始實做中，我們使用 `long long` 資料長度做運算，即64位元有號整數，其有效範圍是 $2^{64 - 1}-1$ 至 $-2^{64}$ 之間，根據運算結果在 $F(93)$ 會出現 overflow。
```
F(91)    = 4660046610375530309 
F(92)    = 7540113804746346429

2^63 - 1 = 9223372036854775808

F(93)    = 12200160415121876738
```
實際運算結果則是
```
F(92)    = 7540113804746346429
F(93)    = -6246583658587674878
F(94)    = 1293530146158671551
```

因此若要運算 $F(93)$ 以後的數值，我們必需設計大數運算的機制。


## 解決方法

### 定義 `bn` 結構體

為了計算 $F(93)$ 以後的 Fibonacci seqeunce，我們運用動態配置記憶體來獲得我們所需的資料長度。

```c
typedef struct {
    unsigned long long *number;
    unsigned int size;
    unsigned int capacity;
} bn;
```
* number : 藉由動態配置特定長度的`long long`陣列來紀錄其數值
* size : 目前數值使用到的`long long`數量
* capacity : 目前結構中配置的`long long`數量


#### bn 的加法實作

```c
void bn_add(bn *out, bn *a, bn *b)
{
    int short_size, long_size;
    bn *larger;
    if (b->size > a->size) {
        short_size = a->size;
        long_size = b->size;
        larger = b;
    } else {
        short_size = b->size;
        long_size = a->size;
        larger = a;
    }

    bn tmp = BN_INIT;
    tmp.number =
        kmalloc(larger->capacity * sizeof(unsigned long long), GFP_KERNEL);
    memset(tmp.number, 0, larger->capacity * sizeof(unsigned long long));
    tmp.capacity = larger->capacity;

    int carry = 0;
    for (int i = 0; i < short_size; i++) {
        unsigned long long x = a->number[i];
        unsigned long long y = b->number[i];
        tmp.number[i] = x + y + carry;
        if (y > (~0UL - x - carry)) {
            carry = 1;
        } else {
            carry = 0;
        }
    }
    for (int i = short_size; i < long_size; i++) {
        unsigned long long x = larger->number[i];
        if (x == (~0UL)) {
            tmp.number[i] = 0;
            carry = 1;
        } else {
            tmp.number[i] = x + carry;
            carry = 0;
        }
    }
    tmp.size = larger->size;
    if (carry) {
        if (tmp.size == tmp.capacity)
            bn_extend(&tmp);
        tmp.number[larger->size] = 1;
        tmp.size++;
    }
    bn_set_bn(out, &tmp);
}
```

#### bn 的減法實作
```c
void bn_sub(bn *out, bn *a, bn *b)
{
    int short_len = b->size;
    bn *larger = a, *smaller = b;

    bn tmp = BN_INIT;
    bn_copy(&tmp, larger);

    for (int i = 0; i < short_len; i++) {
        if (tmp.number[i] < smaller->number[i] && i + 1 < tmp.capacity) {
            int j = i + 1;
            while (j < tmp.capacity && tmp.number[j] == 0) {
                tmp.number[j++]--;
            }
            if (j < tmp.capacity)
                tmp.number[j]--;
        }
        tmp.number[i] -= smaller->number[i];
    }
    bn_set_bn(out, &tmp);
}
```

#### bn 的左移實作
```c
void bn_lshift(bn *out, bn *num)
{
    if (!num)
        return;

    bn tmp = BN_INIT;
    bn_copy(&tmp, num);

    if (tmp.number[tmp.capacity - 1] & (1UL << 63)) {
        bn_extend(&tmp);
    }

    int tail = 0, head;
    for (int i = 0; i < tmp.size; i++) {
        if (tmp.number[i] & (1UL << 63))
            head = 1;
        else
            head = 0;
        tmp.number[i] <<= 1;
        tmp.number[i] += tail;
        tail = head;
    }
    if (tail) {
        tmp.number[num->size] = 1;
        tmp.size++;
    }
    bn_set_bn(out, &tmp);
}
```

#### bn 的乘法實作
```c 
void bn_mul(bn *out, bn *a, bn *b)
{
    bn tmp = BN_INIT;
    bn tmp2 = BN_INIT;
    bn ret = BN_INIT;
    ret.number = kmalloc(2 * sizeof(unsigned long long), GFP_KERNEL);
    memset(ret.number, 0, 2 * sizeof(unsigned long long));
    ret.size = 1;
    ret.capacity = 2;

    bn_copy(&tmp, a);
    bn_copy(&tmp2, b);

    for (int i = 0; i < tmp.size; i++) {
        for (int j = 0; j < 64; j++) {
            if (tmp.number[i] & (1UL << j)) {
                bn_add(&ret, &ret, &tmp2);
            }
            bn_lshift(&tmp2, &tmp2);
        }
    }
    kfree(tmp.number);
    kfree(tmp2.number);
    bn_set_bn(out, &ret);
}

```

#### 將 bn 轉換成字串
```c 
char *bn_to_string(bn *num)
{
    int len = (8 * sizeof(unsigned long long) * num->size) / 3 + 2;

    char *s = kmalloc(len, GFP_KERNEL);
    char *p = s;
    memset(s, '0', len - 1);
    s[len - 1] = '\0';

    for (int i = num->size - 1; i >= 0; i--) {
        for (unsigned long long j = 1UL << 63; j; j >>= 1) {
            int carry = !!(j & num->number[i]);
            for (int k = len - 2; k >= 0; k--) {
                s[k] += s[k] - '0' + carry;
                carry = s[k] > '9';
                if (carry)
                    s[k] -= 10;
            }
        }
    }

    while (*p == '0' && *(p + 1) != '\0')
        p++;
    memmove(s, p, strlen(p) + 1);

    return s;
}
```
### 正確計算 $F(93)$ 以後的數值

#### Iterative approach (Dynamic programming)

利用我們實作的 bn 結構及運算，以及 DP 演算法計算 Fibonacci sequence

參考 [你所不知道的 C 語言: goto 和流程控制篇](https://hackmd.io/@sysprog/c-control-flow)，加入 goto 的使用讓程式碼更加簡潔
```c 
void fib_bn_dp(unsigned long long target, char *buf)
{
    char *ret;
    bn *a;
    bn_init_u64(&a, 1);
    bn *b;
    bn_init_u64(&b, 1);
    bn *c;
    bn_init_u64(&c, 0);
    if (target <= 2)
        goto end;

    for (int i = 3; i <= target; i++) {
        bn_add(c, a, b);
        bn_swap(a, c);
        bn_swap(b, c);
    }

end:
    ret = bn_to_string(c);
    copy_to_user(buf, ret, strlen(ret) + 1);
    kfree(ret);
    bn_free(a);
    bn_free(b);
    bn_free(c);
}
```


#### Fast doubling approach

利用我們實作的 bn 結構及運算，搭配 fast doubling 演算法計算 Fibonacci sequence

```c 
void bn_fast_doubling(unsigned long long target, char *buf)
{
    char *ret;
    bn *fib_n;
    bn_init_u64(&fib_n, 1);
    bn *fib_n1;
    bn_init_u64(&fib_n1, 1);
    bn *fib_2n;
    bn_init_u64(&fib_2n, 0);
    bn *fib_2n1;
    bn_init_u64(&fib_2n1, 0);

    if (target <= 2)
        goto end;

    int count = 63 - __builtin_clzll(target);

    for (unsigned long long i = count; i-- > 0;) {
        bn_lshift(fib_2n, fib_n1);
        bn_sub(fib_2n, fib_2n, fib_n);
        bn_mul(fib_2n, fib_2n, fib_n);

        bn_mul(fib_2n1, fib_n, fib_n);
        bn_mul(fib_n1, fib_n1, fib_n1);
        bn_add(fib_2n1, fib_2n1, fib_n1);
        if (target & (1UL << i)) {
            bn_copy(fib_n, fib_2n1);
            bn_add(fib_n1, fib_2n, fib_2n1);
        } else {
            bn_copy(fib_n, fib_2n);
            bn_copy(fib_n1, fib_2n1);
        }
    }
end:
    ret = bn_to_string(fib_n);
    copy_to_user(buf, ret, strlen(ret) + 1);
    kfree(ret);
    bn_free(fib_n);
    bn_free(fib_n1);
    bn_free(fib_2n);
    bn_free(fib_2n1);
}
```

經過以上兩個演算法輸出結果，至 $F(100000)$ 都可以得到正確結果。

#### 利用 write 設定要使用的 Fibonacci sequence 演算法

在這個 kernel module 跟 user space 程式的互動中，user space client 並沒有寫資料到kerenl module 的需求，因此我們先設定 offset 再利用 write 讓 kernel module 讀到設定的 offset，藉由這個 offset 設定 kernel module 中所使用的 Fibonacci sequence 演算法。

```c 
static void (*fib_func)(unsigned long long, char *);

...
    
static ssize_t fib_write(struct file *file,
                         const char *buf,
                         size_t size,
                         loff_t *offset)
{
    switch (*offset) {
    case 0:
        fib_func = fib_bn_dp;
        break;
    case 1:
        fib_func = bn_fast_doubling;
        break;
    default:
        return 0;
    }
    return 1;
}
```

## 效能分析

計算 kenel module 執行時間我們使用 ktime_t 這個資料結構，以及基於 ktime_t 的 API。

對 ktime_t 做運算：
```cpp
ktime_t ktime_add(const ktime_t add1, const ktime_t add2);
ktime_t ktime_sub(const ktime_t lhs, const ktime_t rhs);
u64 ktime_to_ns(ktime_t kt);
```

在實際計算的 function 前後，我們用 `ktime_get()` 得到當前的時間，將兩者時間相減之後再將 ktime 轉換成相對應的 nanosecond 數值，並將測量到的時間作為回傳值 return 到 user space client。

```c
static ktime_t start_time, end_time;
static s64 elapsed_ns;

static ssize_t fib_read(struct file *file,
                        char *buf,
                        size_t size,
                        loff_t *offset)
{
    return fib_time_proxy(*offset, buf);
}

static long long fib_time_proxy(unsigned long long k, char *buf)
{
    start_time = ktime_get();
    fib_func(k, buf);
    end_time = ktime_get();
    elapsed_ns = ktime_to_ns(ktime_sub(end_time, start_time));

    return (long long) elapsed_ns;
}
```

接著我們撰寫 shell script 將 client 輸出的數值寫到對應的 file 中。在此我們蒐集 0 到 10000 (每間距為 2 收集一次) 的 Fibonacci 數值的運算時間，每個 Fibonacci number 收集 100 次資料，再做後續的統計分析。
```
#!/bin/bash

# script for recording kernel execution time
# first arg is the fib function number, which 0 is dp and 1 is fast doubling.
# second arg is the offset

cd fd_data
for ((i = 0; i <= 1000; i+=2))
do
    filename="${i}.txt"
    touch "$filename"

    sudo ~/linux2023/fibdrv/client 1 "$i" >> "$filename"
done
```

收集到的資料使用 python script做後續分析。先用 [z-score](https://www.investopedia.com/terms/z/zscore.asp)去除離群值，剩下的資料取平均後再作圖展示結果。

```python 
def remove_outliers(data, threshold = 2):
  z_scores = np.abs(stats.zscore(data))
  filtered_data = data[(z_scores < threshold)]
  return filtered_data
```

![](https://hackmd.io/_uploads/rkrH5InL3.png)

結果顯示 fast doubling 的耗費時間比 dynamic programming 的方法還要多，跟演算法時間複雜度預期的結果不符合。
暫時推測原因：
* 相較於 DP 演算法只有使用加法， fast doubling 運算過程使用大量乘法，然而我們目前實作的傳統乘法過程中又含有許多位移和加法操作，導致整體時間比 DP 實作還要慢
* 沒有考慮系統上其他程式的執行，導致執行時間受到影響

我們若要驗證以上推測以及優化執行時間，需要完成以下事項：
- [ ]   尋找比傳統乘法更好的乘法演算法
- [ ] 實作平方運算演算法，最佳化目前問題所需要的乘法情境
- [x] 嘗試將程式限縮在特定CPU上執行
- [ ] 參考其他 big number 實作，例如 [sysprog21/bignum](https://github.com/sysprog21/bignum)以及[The GNU
Multiple Precision Arithmetic Library](https://gmplib.org/)尋求改進空間
- [ ] 降低 copy_to_user 資料量，例如實作資料壓縮

### 保留特定 CPU 來執行 Process

taskset 可以查看指定 Process 的 CPU affinity
```shell 
taskset -p PID
```

也可以用 coremask 或是 corelist 將 precess 固定在特定的 CPU 上執行

```shell 
taskset -p COREMASK PID
taskset -cp CORELIST PID
```

例如
```shell 
taskset -p 0x11 1234
taskset -cp 0,4,5,6 1234
```

亦可指定 CPU 執行特定執行檔
```shell 
taskset COREMASK EXECUTABLE
```
雖然可以將特定 process 保留在特定 CPU 上，但我們沒辦法保證其他 process 不會在這些 CPU  上執行，因此需要在 GRUB 設定檔中加上 `isolcpus=cpu_id` 參數

```shell
sudo vim /etc/default/grub
GRUB_CMDLINE_LINUX_DEFAULT="quiet splash isolcpus=0,1"

sudo update-grub
sudo reboot
```

重新開機後，我們在原本的 shell script 中導入 taskset
```shell 
cd fd_data
for ((i = 0; i <= 1000; i+=2))
do
    filename="${i}.txt"
    touch "$filename"

    sudo taskset -c 0,1 ../client 1 "$i" >> "$filename"
done
```
來收集較無受到其他 process 影響的執行時間
![](https://hackmd.io/_uploads/Sk6aSsaL2.png)

由圖片結果，相較於沒有保留 CPU 來計算執行時間我們看到顯著的改善，執行時間的趨勢更加明顯

### 利用 perf 找出程式執行瓶頸

使用 [perf](https://perf.wiki.kernel.org/index.php/Main_Page) 工具找出程式執行瓶頸，作為後續改善之基準

參考資料：
* [Perf Wiki by kernel.org](https://perf.wiki.kernel.org/index.php/Main_Page)
* [Linux 效能分析工具: Perf](http://wiki.csie.ncku.edu.tw/embedded/perf-tutorial)
```shell 
sudo perf record -g ./client
sudo perf report --stdio -g graph,0.5
```

![](https://hackmd.io/_uploads/r1yWs6RL3.png)

根據 `perf` 紀錄的結果，我們可以發現有超過 60% 的時間都花在把二進位數字轉換成十進位上，跟原本預期的結果不同，可以作為未來改善效能的依據



:::info
TO DO
:::



## Kernel Module 除錯

參考
* Linux Foundation 演講 [Linux Kernel Debugging: Going Beyond Printk Messages - Sergio Prado, Embedded Labworks](https://www.youtube.com/watch?v=NDXYpR_m1CU&t=238s&ab_channel=TheLinuxFoundation)
* [建構 User-Mode Linux 的實驗環境](https://hackmd.io/@sysprog/user-mode-linux-env)

addr2line

cppcheck - static analysis

valgrind - dynamic analysis

:::info
TO DO
:::


## 多執行緒環境下的問題及設計

client 在執行 `read` 取得 Fibonacci number 結果前會先使用 `lseek` 設定 offset，在多執行緒環境下執行 `read` 前若有其他執行緒呼叫 `lseek`， offset 會被重新設定導致得到錯誤的 Fibonacci number 結果。因此，我們必需對 kernel module 提供的服務做資源保護。

參考以下內容：
[並行和多執行緒程式設計](https://hackmd.io/@sysprog/concurrency/https%3A%2F%2Fhackmd.io%2F%40sysprog%2FS1AMIFt0D)
[Linux 核心設計: 淺談同步機制](https://hackmd.io/@sysprog/linux-sync)
[Linux 核心設計: 多核處理器和 spinlock](https://hackmd.io/@sysprog/multicore-locks)

我們使用 mutex lock 機制在 kernel module 被 `open` 的時候上鎖，在 `close` 的時候解鎖，確保 kernel module 的 function 在執行中只有一個執行緒可以取得其服務。

```c 
static int fib_open(struct inode *inode, struct file *file)
{
    if (!mutex_trylock(&fib_mutex)) {
        printk(KERN_ALERT "fibdrv is in use");
        return -EBUSY;
    }
    return 0;
}

...
    
static int fib_release(struct inode *inode, struct file *file)
{
    mutex_unlock(&fib_mutex);
    return 0;
}
```

在 client 端，我們使用 pthread API 創造執行緒以驗證 kernel module 在多執行緒執行環境中得到正確的結果。

```c 
typedef struct {
    int func;
    int offset;
} thread_arg;

...

int main(int argc, char *argv[])
{
    pthread_t threads[THREAD_NUM];
    thread_arg args[THREAD_NUM];

    for (int i = 0; i < THREAD_NUM; i++) {
        args[i].func = atoi(argv[1]);
        args[i].offset = atoi(argv[0]) - i;
    }

    for (int i = 0; i < THREAD_NUM; i++) {
        pthread_create(&threads[i], NULL, thread_func, &args[i]);
    }

    for (int i = 0; i < THREAD_NUM; i++) {
        pthread_join(threads[i], NULL);
    }

    return 0;
}
```

:::info
TO DO
- 將 lock 做更細緻的處理(fine grained)是否有效改善運算效能
:::


