# 2023q1 Homework3 (fibdrv)
contributed by < [`ericlai1021`](https://github.com/ericlai1021/fibdrv) >
## 開發環境
```shell
$ 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: 39 bits physical, 48 bits virtual
Byte Order: Little Endian
CPU(s): 8
On-line CPU(s) list: 0-7
Vendor ID: GenuineIntel
Model name: Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz
CPU family: 6
Model: 142
Thread(s) per core: 2
Core(s) per socket: 4
Socket(s): 1
Stepping: 10
CPU max MHz: 3400.0000
CPU min MHz: 400.0000
BogoMIPS: 3600.00
Virtualization features:
Virtualization: VT-x
Caches (sum of all):
L1d: 128 KiB (4 instances)
L1i: 128 KiB (4 instances)
L2: 1 MiB (4 instances)
L3: 6 MiB (1 instance)
NUMA:
NUMA node(s): 1
NUMA node0 CPU(s): 0-7
$ uname -r
5.19.0-35-generic
```
## 前期準備
* 取得[原始程式碼](https://github.com/ericlai1021/fibdrv)後，進行編譯及測試
* 觀察產生的 `fibdrv.ko` 核心模組
```shell
$ modinfo fibdrv.ko
```
得到以下輸出
```shell
filename: /home/eric/linux2023/fibdrv/fibdrv.ko
version: 0.1
description: Fibonacci engine driver
author: National Cheng Kung University, Taiwan
license: Dual MIT/GPL
srcversion: 9A01E3671A116ADA9F2BB0A
depends:
retpoline: Y
name: fibdrv
vermagic: 5.19.0-35-generic SMP preempt mod_unload modversions
```
* 觀察 `fibdrv.ko` 核心模組在 Linux 核心掛載後的行為
模組載入核心：
```shell
$ sudo insmod fibdrv.ko
```
* 查看 `/dev` 目錄下可以發現的確有一個名為 `fibonacci` 的 character device
```shell
$ ls -l /dev/fibonacci
crw------- 1 root root 509, 0 三 16 13:44 /dev/fibonacci
```
繼續觀察
```shell
$ cat /sys/class/fibonacci/fibonacci/dev
509:0
```
可以注意到 `509` 這個數字，不同電腦設備會產生不同數字，試著對照 [`fibdrv.c`](https://github.com/ericlai1021/fibdrv/blob/master/fibdrv.c) ，發現關鍵在 alloc_chrdev_region() ， 此函式會向 linux kernel 申請設備號。
```shell
$ cat /sys/module/fibdrv/version
0.1
```
這和 [`fibdrv.c`](https://github.com/ericlai1021/fibdrv/blob/master/fibdrv.c) 透過 `MODULE_VERSION` 所指定的版本號碼相同。
```shell
$ lsmod | grep fibdrv
fibdrv 16384 0
$ cat /sys/module/fibdrv/refcnt
0
```
這兩道命令的輸出都是 0，意味著目前的 [reference counting](https://en.wikipedia.org/wiki/Reference_counting)。
* 掛載 `fibdrv.ko` 核心模組後，執行 `client` 並觀察輸出
```shell
$ sudo ./cilent
Reading from /dev/fibonacci at offset 0, returned the sequence 0.
Reading from /dev/fibonacci at offset 1, returned the sequence 1.
Reading from /dev/fibonacci at offset 2, returned the sequence 1.
...
Reading from /dev/fibonacci at offset 92, returned the sequence 7540113804746346429.
Reading from /dev/fibonacci at offset 93, returned the sequence 7540113804746346429.
Reading from /dev/fibonacci at offset 94, returned the sequence 7540113804746346429.
```
可以發現在 offset 大於 92 之後都會得到同樣的值，猜測是因為 overflow ，於是先去查看 [`client.c`](https://github.com/ericlai1021/fibdrv/blob/master/client.c) ，基本上該檔案只是對 `/dev/fibonacci` 此裝置做開檔、讀檔、寫檔的動作，所以 fibonacci sequence 的計算都在 [`fibdrv.c`](https://github.com/ericlai1021/fibdrv/blob/master/fibdrv.c) 檔案內，馬上去查看，發現以下程式碼
```c
/* MAX_LENGTH is set to 92 because
* ssize_t can't fit the number > 92
*/
#define MAX_LENGTH 92
```
再往下看到 `fib_device_lseek()` ，此函式中會判斷欲找尋的位置是否大於 `MAX_LENGTH` ，若大於則回傳 `MAX_LENGTH` ，因此要讀取大於 F(92) 的數值都會讀到 F(92)。
## Linux 效能分析
參考 [作業說明 （三）](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-c) 當中的設定
### 限定 CPU 給特定的程式使用
使用 `isolcpus` 這個 Linux 核心起始參數，修改 `/etc/default/grub` 設定檔，在 `GRUB_CMDLINE_LINUX_DEFAULT` 欄位加入 `isolcpus=0` 使開機時保留第 0 個 CPU 核給那些被 [taskset](https://man7.org/linux/man-pages/man1/taskset.1.html) 特別指定的行程使用。
```shell
GRUB_CMDLINE_LINUX_DEFAULT="quiet splash isolcpus=0"
```
修改完輸入 `sudo update-grub` 以更新設定，重新開機後確認設定是否生效。
```shell
$ cat /sys/devices/system/cpu/isolated
0
$ taskset -cp 1
pid 1's current affinity list: 1-7
```
可以發現第 0 個 CPU 確實被保留下來，再使用 `taskset` 查看 PID=1 的行程的處理器親和性，此時已經不包含第 0 個 CPU。
### 時間量測
使用 [ktimer 相關的 API](https://www.kernel.org/doc/html/latest/core-api/timekeeping.html) 來量測時間。 改寫 `fibdrv.c` 當中的 `fib_read()` 和 `fib_write()` 兩個函式。
首先宣告一個 `ktime_t` 的變數:
```c
static ktime_t kt;
```
針對 `fib_read()` 的修改如下：
```c
static long long fib_time_proxy(long long k)
{
kt = ktime_get();
long long result = fib_sequence(k);
kt = ktime_sub(ktime_get(), kt);
return result;
}
/* calculate the fibonacci number at given offset */
static ssize_t fib_read(struct file *file,
char *buf,
size_t size,
loff_t *offset)
{
return (ssize_t) fib_time_proxy(*offset);
}
```
針對 `fib_write()` 的修改如下：
```c
/* write operation is skipped */
static ssize_t fib_write(struct file *file,
const char *buf,
size_t size,
loff_t *offset)
{
return ktime_to_ns(kt);
}
```
我們可發現到 write 在此核心模組暫時沒作用，所以可以用來輸出上一次執行 fib_sequence() 的時間。
在 `client.c` 中先用 lseek(fd, i, SEEK_SET) 將檔案的讀寫位置指向特定位置 （即第 i 個位置），接著再用 read() 取得 Fib(i) 的數值，最後再使用 write() 取得計算 Fib(i) 的時間，相對應的程式碼如下：
```c
for (int i = 0; i <= offset; i++) {
lseek(fd, i, SEEK_SET);
sz = read(fd, buf, 1);
sz = write(fd, write_buf, strlen(write_buf));
}
```
## 改進費氏數列
### Fast Doubling
參考作業說明 [(一)](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-a#-%E8%B2%BB%E6%B0%8F%E6%95%B8%E5%88%97) 、 [(四)](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-d) 以迭代方式實作 bottom-up Fast Doubling。
使用 Fast Doubling 可以得出以下公式
$$
\begin{split}
F(2k) &= F(k)[2F(k+1) - F(k)] \\
F(2k+1) &= F(k+1)^2+F(k)^2
\end{split}
$$
對應的虛擬碼如下：
```python
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;
```
搭配 bitwise 操作求得 `number of binary digit in n` 以及判斷 `current binary digit` 是否為 1。
* 首先討論如何判斷當前位元是否為 1，想法很簡單，就是將 1 往左位移 i 個位元再跟 n 做 AND 運算即可求得 n 的第 i 位元是否為 1。
對應的程式碼如下：
```c
...
if (n & (1 << i)) {
...
}
...
```
* 接著來討論如何求得 n 的二進制位數，這邊可以選擇是否搭配 [clz / ctz](https://en.wikipedia.org/wiki/Find_first_set) 一類的指令來實作。
* 使用 clz 指令取得 n 的 leading 0-bits 數量
因為 n 的型態為 `long long` ，所以首先將 `sizeof(long long)` 向左位移 3 個位元即可求得 `long long` 所需的位元數，接著再減去 n 的 leading 0-bits 數量就可以求得 n 的二進制位數。
對應的程式碼如下：
```c
/* calculate the number of binary digits in k using clz */
int nbits = (sizeof(n) << 3) - __builtin_clzll(n);
```
* 相反的，如果不使用 clz 指令的話，我的作法是簡單的透過一個迴圈產生跟 clz 指令同樣的結果
想法是從 n 的 MSB 往下逐一判斷並搭配變數 count 表示 n 的 leading 0-bits 數，若當前位元為 1 則跳出迴圈，否則讓 `count++` ， 最後 `(sizeof(long long) << 3) - count` 即為 n 的 leading 0-bits 數量。
對應的程式碼如下：
```c
/* calculate the number of binary digits in k without using clz */
int nbits = 0;
int llsize = sizeof(n) << 3;
for (int i = llsize - 1; i >= 0; i--) {
if (n >> i)
break;
nbits++;
}
nbits = llsize - nbits;
```
### 實驗結果 （尚未排除干擾效能分析的因素）
* 比較使用遞迴方法與 Fast Doubling （搭配 clz 指令）計算費氏數列的時間
可以發現 Fast Doubling 在計算 Fib(0) 時會花費大量時間，猜測是跟使用 clz 指令有關，目前正在做實驗研究，若我們將 y 軸範圍縮小到 2500 ns 會得到以下結果
即便尚未排除干擾效能的因素，仍可以看到 Fast Doubling 明顯優於遞迴方式
* 目前對於 Fast Doubling 計算 Fib(0) 花費大量時間的解法為在 `fib_sequence()` 的一開始加上提前中止的條件處理
```c
if (k <= 0)
return k;
```
### 排除干擾效能分析的因素
* 抑制 [address space layout randomization](https://en.wikipedia.org/wiki/Address_space_layout_randomization) (ASLR)
```shell
$ sudo sh -c "echo 0 > /proc/sys/kernel/randomize_va_space"
```
* 設定 scaling_governor 為 performance。準備以下 shell script，檔名為 `performance.sh`:
```shell
for i in /sys/devices/system/cpu/cpu*/cpufreq/scaling_governor
do
echo performance > ${i}
done
```
之後再用 `$ sudo sh performance.sh` 執行
* 針對 Intel 處理器，關閉 turbo mode
```shell
$ sudo sh -c "echo 1 > /sys/devices/system/cpu/intel_pstate/no_turbo"
```
* 關閉 [SMT](https://en.wikipedia.org/wiki/Simultaneous_multithreading) (Intel 稱它為 [Hyper-Threading](https://en.wikipedia.org/wiki/Hyper-threading))
```shell
$ sudo sh -c "echo off > /sys/devices/system/cpu/smt/control"
```
### 實驗結果 （排除干擾效能分析的因素）
* 比較使用遞迴方法與 Fast Doubling （搭配 clz 指令）計算費氏數列的時間
可以看出來兩種方式的曲線都明顯平滑許多
* 比較在 Fast Doubling 上有無使用 clz 指令的差別
可以看到使用 clz 指令對計算費氏數列有顯著的改善。
## 計算 Fib(92) 之後的 Fibonacci 數
> 參考作業說明 [(四)](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-d#%E5%88%9D%E6%AD%A5%E6%94%AF%E6%8F%B4%E5%A4%A7%E6%95%B8%E9%81%8B%E7%AE%97) 初步支援大數運算
### 支援大數運算結構體
為了能正確計算 92 項以後的費氏數列，引入長度可變動的數值表示法，藉由動態配置不同大小記憶體來呈現更大範圍的整數，結構體的定義如下
值得注意的地方是，我透過觀察發現即便是使用 [fast doubling](https://hackmd.io/@ericlai1021/linux2023q1-fibdrv#Fast-Doubling) 在計算費氏數列的過程中，所有會用到自定義結構體 `bn` 宣告的變數 (e.g., `a`, `b`, `t1`, `t2` ) 皆為非負整數，所以我就沒有定義 `sign` 來表示數值的正負。
```c
/*
* number[size - 1] = msb, number[0] = lsb
* sign = 1 for negative number
*/
typedef struct _bn {
unsigned int *number;
unsigned int size;
} bn;
```
* `number` : 指向無號整數陣列的開頭
* `size` : 配置給 number array 的大小，單位為 `sizeof(unsigned int)`
### iterative 實作方法
初步嘗試使用大數運算來實作迭代作法的費氏數列計算，首先會需要加入以下功能
1. `bn_to_string` ，將 `bn` 型態變數轉換成字串
2. `bn_add` ，將兩個 `bn` 型態變數相加
3. `bn_cpy` ，將一個 `bn` 型態變數複製到另一個 `bn` 型態變數（即做 assign 的動作）
在講解以上功能之前要先介紹支援 `bn` 型態的記憶體配置方法以及釋放方法，這邊為了能更加彈性的使用 `bn` 結構體，因此皆採用動態配置記憶體的方式來使用 `bn` 結構體
### bn_alloc
先配置空間給 `bn` 結構體，再配置 `sizeof(unsigned int) * size` 的大小給 `number` 指標變數，最後將所有數值初始化為 0
參考教材 [kmalloc vs. vmalloc](https://hackmd.io/@sysprog/linux-memory#kmalloc-vs-vmalloc) 及網路材料 [difference between kmalloc and vmalloc](https://www.emblogic.com/blog/10/difference-between-kmalloc-and-vmalloc/) 後選擇使用 kmalloc 來配置記憶體是因為 kmalloc 會配置實體連續的記憶體空間，接著就可以簡單的使用 `memset` 來初始化此空間，但當沒有足夠的實體記憶體空間時 kmalloc 就會回傳錯誤，通常是回傳 `NULL`
```c
/*
* allocate a bn structure with the given size
* the value is initialized to 0
*/
bn *bn_alloc(size_t size)
{
bn *new = (bn *) kmalloc(sizeof(bn), GFP_KERNEL);
new->number =
(unsigned int *) kmalloc(sizeof(unsigned int) * size, GFP_KERNEL);
// QUESTION? physical or virtual contiguous memory allocation
memset(new->number, 0, sizeof(unsigned int) * size);
new->size = size;
return new;
}
```
### bn_free
依序將 `number` 以及 `bn` 結構體釋放
```c
/*
* free entire bn data structure
* return 0 on success, -1 on error
*/
int bn_free(bn *src)
{
if (src == NULL)
return -1;
kfree(src->number);
kfree(src);
return 0;
}
```
### bn_to_string
由於大數沒辦法直接以數值的形式列出，這裡改用字串來呈現，轉換的概念很簡單，就是將大數除以 10 取餘數取得個位數，接著再將大數更新為商數，不斷重複此步驟直到大數為 0
```c
memset(s, '0', len - 1);
s[len - 1] = '\0';
int cur = len - 2;
while (!bn_is_zero(src)) {
unsigned long long remainder = 0U;
// src->number[size - 1] is msb
for (int i = src->size - 1; i >= 0; i--) {
unsigned int quotient = 0U;
// walk through every bit of bn
for (unsigned int d = 1U << 31; d > 0; d >>= 1) {
if (src->number[i] & d) {
remainder = (remainder << 1) | 1U;
} else {
remainder = (remainder << 1) | 0U;
}
quotient = (quotient << 1) | (remainder / 10U);
remainder %= 10;
}
src->number[i] = quotient;
}
s[cur] += remainder;
cur--;
}
```
### bn_add
實作加法的概念很直觀，就是從 lsb 開始相加，將進位的部份加到下一輪的相加結果
```c
/* c = a + b */
void bn_add(const bn *a, const bn *b, bn *c)
{
bn_resize(c, MAX(a->size, b->size) + 1);
unsigned long long carry = 0;
for (int i = 0; i < c->size; i++) {
unsigned int tmp1 = (i < a->size) ? a->number[i] : 0;
unsigned int tmp2 = (i < b->size) ? b->number[i] : 0;
carry += (unsigned long long) tmp1 + tmp2;
c->number[i] = carry & 0x00000000FFFFFFFF;
carry >>= 32;
}
if (!c->number[c->size - 1] && c->size > 1)
bn_resize(c, c->size - 1);
}
```
### bn_cpy
使用 `memcpy` 將一個 `bn` 結構體的數值部份複製到另一個 `bn` 結構體的數值部份
```c
/*
* copy the value from src to dest
* return 0 on success, -1 on error
*/
int bn_cpy(bn *dest, const bn *src)
{
if (bn_resize(dest, src->size) < 0)
return -1;
memcpy(dest->number, src->number, src->size * sizeof(unsigned int));
return 0;
}
```
### fast doubling 實作方法
使用大數運算實作 fast doubling 會需要額外實作以下功能
1. `bn_sub` ，將兩個 `bn` 型態變數相減
2. `bn_lshift` ，將 `bn` 型態變數向左位移
3. `bn_mult` ，將兩個 `bn` 型態變數相乘
### bn_sub
實作減法的概念類似加法，一樣從 lsb 開始相減，若相減結果為負數，則向下一位借位
```c
/* c = a - b */
void bn_sub(const bn *a, const bn *b, bn *c)
{
bn_resize(c, MAX(a->size, b->size));
long long borrow = 0;
for (int i = 0; i < c->size; i++) {
unsigned int tmp1 = (i < a->size) ? a->number[i] : 0;
unsigned int tmp2 = (i < b->size) ? b->number[i] : 0;
borrow = (long long) tmp1 - tmp2 - borrow;
if (borrow < 0) {
c->number[i] = borrow + (1LL << 32);
borrow = 1;
} else {
c->number[i] = borrow;
borrow = 0;
}
}
int clz_ints = bn_clz(c) / 32;
if (clz_ints == c->size)
--clz_ints;
bn_resize(c, c->size - clz_ints);
}
```
### bn_shift
對 `bn` 結構體左移 `shift` 個 bit ，因為左移超過 32 bits 的處理相對複雜，因此就先做出比較簡單的版本
```c
/* dest = src << shift */
void bn_lshift(bn *dest, bn *src, size_t shift)
{
size_t z = bn_clz(src);
shift %= 32; // only handle shift within 32 bits
if (!shift)
return;
if (shift > z)
bn_resize(src, src->size + 1);
bn_cpy(dest, src);
/* bit shift */
for (int i = src->size - 1; i > 0; i--)
dest->number[i] =
src->number[i] << shift | src->number[i - 1] >> (32 - shift);
dest->number[0] <<= shift;
}
```
### bn_mult
實作大數乘法的概念類似硬體乘法器 （如下圖），從 `multiplier` 的 lsb 開始判斷，若當前 bit 為 1 ，則將 `multiplicand` 加至 `product` ，再將 `product` 往左位移 1 bit ，重複上述直到 `multiplier` 的判斷到達 msb 為止
```c
/* c = a * b */
bn *tmp1 = bn_alloc(1);
bn_cpy(tmp1, a);
if (b->number[0] & 1U)
bn_add(c, tmp1, c);
for (int i = 1; i < bn_digit(b); i++) {
bn_lshift(tmp1, tmp1, 1);
if (b->number[i / 32] & (1U << (i % 32)))
bn_add(c, tmp1, c);
}
bn_free(tmp1);
int clz_ints = bn_clz(c) / 32;
if (clz_ints == c->size)
--clz_ints;
bn_resize(c, c->size - clz_ints);
```
### 正確計算 $F(92)$ 以後的數值
使用以上實作出來的大數 API 來正確計算 92 項以後的費氏數列，程式邏輯皆參照原始的實作，首先展示迭代演算法
```c
void bn_fib(bn *dest, unsigned int n)
{
bn_resize(dest, 1);
if (n < 2) { // Fib(0) = 0, Fib(1) = 1
dest->number[0] = n;
return;
}
bn *a = bn_alloc(1);
bn *b = bn_alloc(1);
dest->number[0] = 1;
for (unsigned int i = 1; i < n; i++) {
bn_cpy(b, dest);
bn_add(dest, a, dest);
bn_cpy(a, b);
}
bn_free(a);
bn_free(b);
}
```
接著是 fast doubling 演算法的實作
```c
void bn_fib_fast(bn *dest, unsigned int n)
{
bn_resize(dest, 1);
if (n < 2) { // Fib(0) = 0, Fib(1) = 1
dest->number[0] = n;
return;
}
int nbits = (sizeof(n) << 3) - __builtin_clz(n);
bn *a = bn_alloc(1);
a->number[0] = 0U;
bn *b = bn_alloc(1);
b->number[0] = 1U;
for (int i = nbits - 1; i >= 0; i--) {
bn *t1 = bn_alloc(1);
bn_lshift(t1, b, 1);
bn_sub(t1, a, t1);
bn_mult(a, t1, t1);
bn *t2 = bn_alloc(1);
bn_mult(b, b, b);
bn_mult(a, a, a);
bn_add(a, b, t2);
bn_cpy(a, t1);
bn_cpy(b, t2);
if (n & (1U << i)) {
bn_add(a, b, t1);
bn_cpy(a, b);
bn_cpy(b, t1);
}
bn_free(t1);
bn_free(t2);
}
bn_cpy(dest, a);
bn_free(a);
bn_free(b);
}
```
使用 [/scripts/verify.py](https://github.com/ericlai1021/fibdrv/blob/master/scripts/verify.py) 來驗證，至少能夠正確計算至第 3000 項
## 改善大數運算
後續改進皆更新在 [2023期末專題](https://hackmd.io/@sysprog/HyyMNz7S3)
Sign in with Wallet
Connect another wallet