[我的實作 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
```
結果顯示 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 影響的執行時間
由圖片結果，相較於沒有保留 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
```
根據 `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)是否有效改善運算效能
:::
Sign in with Wallet
Connect another wallet