# 2023q1 Homework3 (fibdrv)
contributed by < [fewletter](https://github.com/fewletter/fibdrv) >
## 實驗環境
```
$ uname -r
5.15.0-67-generic
$ gcc --version
gcc (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4.0
$ lscpu
Architecture: x86_64
CPU op-mode(s): 32-bit, 64-bit
Byte Order: Little Endian
Address sizes: 39 bits physical, 48 bits virtual
CPU(s): 6
On-line CPU(s) list: 0-5
Thread(s) per core: 1
Core(s) per socket: 6
Socket(s): 1
NUMA node(s): 1
Vendor ID: GenuineIntel
CPU family: 6
Model: 158
Model name: Intel(R) Core(TM) i5-9500 CPU @ 3.00GHz
Stepping: 10
CPU MHz: 3000.000
CPU max MHz: 4400.0000
CPU min MHz: 800.0000
BogoMIPS: 6000.00
Virtualization: VT-x
L1d cache: 192 KiB
L1i cache: 192 KiB
L2 cache: 1.5 MiB
L3 cache: 9 MiB
NUMA node0 CPU(s): 0-5
```
> 作業要求
:::info
- [ ] 研讀上述 ==Linux 效能分析的提示== 描述，在自己的實體電腦運作 GNU/Linux，做好必要的設定和準備工作 $\to$ 從中也該理解為何不希望在虛擬機器中進行實驗;
- [x] 研讀上述費氏數列相關材料 (包含論文)，摘錄關鍵手法，並思考 [clz / ctz](https://en.wikipedia.org/wiki/Find_first_set) 一類的指令對 Fibonacci 數運算的幫助。請列出關鍵程式碼並解說
- [x] 複習 C 語言 [數值系統](https://hackmd.io/@sysprog/c-numerics) 和 [bitwise operation](https://hackmd.io/@sysprog/c-bitwise)，思考 Fibonacci 數快速計算演算法的實作中如何減少乘法運算的成本;
- [ ] 學習指出針對大數運算的加速運算和縮減記憶體操作成本的舉措，紀錄你的認知和疑惑
- [ ] 注意到 `fibdrv.c` 存在著 `DEFINE_MUTEX`, `mutex_trylock`, `mutex_init`, `mutex_unlock`, `mutex_destroy` 等字樣，什麼場景中會需要呢？撰寫多執行緒的 userspace 程式來測試，觀察 Linux 核心模組若沒用到 mutex，到底會發生什麼問題。嘗試撰寫使用 [POSIX Thread](https://en.wikipedia.org/wiki/POSIX_Threads) 的程式碼來確認。
:::
## 開發紀錄
### Linux 時間測量和效能分析
參考[作業說明：Linux 時間測量和效能分析](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-c#-Linux-%E6%95%88%E8%83%BD%E5%88%86%E6%9E%90%E7%9A%84%E6%8F%90%E7%A4%BA)和 [yanjiew](https://hackmd.io/@yanjiew/linux2023q1-fibdrv) 所撰寫的時間測試檔案。
為了要達到時間測量的目的，依照[作業說明：Linux 時間測量和效能分析](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-c#-Linux-%E6%95%88%E8%83%BD%E5%88%86%E6%9E%90%E7%9A%84%E6%8F%90%E7%A4%BA)中的步驟實作
* 抑制 address space layout randomization (ASLR)
```cpp
$ sudo sh -c "echo 0 > /proc/sys/kernel/randomize_va_space"
```
* 設定 scaling_governor 為 performance。準備以下 shell script，檔名為 performance.sh:
```cpp
for i in /sys/devices/system/cpu/cpu*/cpufreq/scaling_governor
do
echo performance > ${i}
done
```
* 在測試時間時，利用 `$ sudo sh performance.sh` 執行，可藉由
```cpp
cat /sys/devices/system/cpu/cpu*/cpufreq/scaling_governor
```
* 檢查 cpu 是否已經進入 performance 模式，測試結果：
```cpp
performance
performance
performance
performance
performance
performance
```
* 針對 Intel 處理器，關閉 turbo mode
```cpp
$ sudo sh -c "echo 1 > /sys/devices/system/cpu/intel_pstate/no_turbo"
```
#### 時間測量
為了能夠在測量程式在 kernel 中運行的時間，將沒有用到的函式 `fib_write` 來作為輸出時間的函式。
首先先了解 `ktime_t` 此 api 如何運作，從 [`ktime_t`](https://docs.kernel.org/core-api/timekeeping.html#c.ktime_get) 中我們可以知道其中的 `ktime_get()` 是用來標記時間，所以只要在程式開始和結束時都使用`ktime_get()` 就可以將程式執行時間擷取出來，程式碼如下：
```cpp
static ktime_t kt;
```
在 `fib_read()` 中計時
```cpp
static ssize_t fib_read(struct file *file,
char *buf,
size_t size,
loff_t *offset)
{
ktime_t start = ktime_get();
long long result = fib_sequence(*offset);
kt = ktime_sub(ktime_get(), start);
return (ssize_t) result;
}
```
從 `fib_write()` 中回傳時間
```cpp
static ssize_t fib_write(struct file *file,
const char *buf,
size_t size,
loff_t *offset)
{
return (ssize_t) kt;
}
```
參考 [yanjiew](https://hackmd.io/@yanjiew/linux2023q1-fibdrv) 所撰寫的時間測試檔案，並將其修改成可以輸出 `data.txt` 來紀錄時間，測試結果
```cpp
data.txt
Time on this sequence:44
Time on this sequence:24
Time on this sequence:32
Time on this sequence:35
Time on this sequence:37
Time on this sequence:39
Time on this sequence:42
Time on this sequence:42
...
```
在 `scripts` 中寫一個簡單的程式可以將上述資料畫圖
```python
data = np.loadtxt('data.txt')
if __name__ == '__main__':
fig = plt.figure()
plt.title("python scripts plot")
plt.xlabel("nth fabonacci")
plt.ylabel("time(ns)")
plt.plot(data,label = 'fast doubling')
plt.legend()
fig.savefig('data.png')
```
針對原始費氏數列的遞迴程式碼進行時間測量就會得到下面的圖
```cpp
static long long fib_sequence(long long k)
{
/* FIXME: C99 variable-length array (VLA) is not allowed in Linux kernel. */
long long f[k + 2];
f[0] = 0;
f[1] = 1;
for (int i = 2; i <= k; i++) {
f[i] = f[i - 1] + f[i - 2];
}
return f[k];
}
```
既然可以從輸出 kernel 的執行時間，那接著便是測試 user space 的時間，利用 `time.h` 中的 [`clock_gettime`](https://man7.org/linux/man-pages/man2/clock_getres.2.html) 進行時間的測量，從其中文件可以看到 `clockid_t` 和 `timespec` 兩個引數型態
```cpp
int clock_gettime(clockid_t clockid, struct timespec *tp)
```
往下搜尋可以看到 `timespec` 這個結構體，記錄著以奈秒(ns)為時間的單位
```cpp
struct timespec {
__kernel_time_t tv_sec; /* seconds */
long tv_nsec; /* nanoseconds */
};
```
要利用此結構體只需要將全部的成員換成 ns 的時間單位就好
```cpp
clock_gettime(CLOCK_MONOTONIC, &t1);
long long time = t1.tv_sec * 1e9 + t1.tv_nsec;
```
`clockid_t` 則是選擇 `CLOCK_MONOTONIC`，主要是因為只要形成 process 開始和結束時候的時間標記，程式碼實作如下:
```cpp
clock_gettime(CLOCK_MONOTONIC, &t1);
sz = read(fd, buf, 1);
clock_gettime(CLOCK_MONOTONIC, &t2);
long long ut = (long long)(t2.tv_sec * 1e9 + t2.tv_nsec)
- (t1.tv_sec * 1e9 + t1.tv_nsec);
```
#### 利用統計方法消除極端值
參考[作業參考:用統計手法去除極端值](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-c#%E7%94%A8%E7%B5%B1%E8%A8%88%E6%89%8B%E6%B3%95%E5%8E%BB%E9%99%A4%E6%A5%B5%E7%AB%AF%E5%80%BC)和 [4ce43c4](https://github.com/colinyoyo26/fibdrv/commit/4ce43c4e7ee0b80c4ce9e6dcd995bbfdf239206c) 來實作 python scripts，實作此 scripts 的目的為降低極端值對圖的影響，如果像上面那張圖只憑計算一次的時間來畫圖，很可能會有極端值出現降低圖的可信度
上圖為原始算法所算的測試時間，在 `kernel` 中的測試時間隨著費氏數列的項數而增加，但是在 `user` 中的時間卻變動不大，這部分還在研究中。
#### 修改 Makefile
為了將上述命令不需要在每次測試程式時都重打一次，將上述命令寫進 Makefile 中並且命名為 test
```cpp
test: all
sudo sh -c "echo 0 > /proc/sys/kernel/randomize_va_space"
sudo sh performance.sh
sudo sh -c "echo 1 > /sys/devices/system/cpu/intel_pstate/no_turbo"
$(MAKE) unload
$(MAKE) load
sudo taskset 0x1 ./test_time
$(MAKE) unload
python3 scripts/plot.py
```
### 研讀費氏數列相關材料
參考[作業說明一](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-a)和 [Fast Fibonacci algorithms](https://www.nayuki.io/page/fast-fibonacci-algorithms) ，可以知道費式數列有多種計算方法，第一種為公式解
$$F(n) = \frac{1}{\sqrt{5}}(\frac{1+\sqrt{5}}{2})^n-\frac{1}{\sqrt{5}}(\frac{1-\sqrt{5}}{2})^n$$
此方法最大的缺點在於 $\sqrt{5}$ 此數字在電腦中計算時無法取得確切數值，所以當計算$(\frac{1+\sqrt{5}}{2})^n$ 時，誤差會因為 n 次方而增加。
第二種方法為直觀的加法，利用費氏數列原本的定義來計算
$$F(k) = F(k-1) + F(k-2)$$
此種方法的缺點則在會重複計算某些不必要的數值，只為了形成 $F(k-1)$ 和 $F(k-2)$，並且每次計算時必須要每項都重新計算一次。
第三種則是 Fast Doubling，為了能夠減少計算量觀察其中規律
$$ \begin{bmatrix}
F(2) \\
F(1)
\end{bmatrix}
= \begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}
\begin{bmatrix}
F(1) \\
F(0)
\end{bmatrix} $$
$$\to\ \begin{bmatrix}
F(n+1) \\
F(n)
\end{bmatrix}
= \begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}^{n}
\begin{bmatrix}
F(1) \\
F(0)
\end{bmatrix} $$
$$\to\ \begin{bmatrix}
F(2n+1) \\
F(2n)
\end{bmatrix}
= \begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}^{2n}
\begin{bmatrix}
F(1) \\
F(0)
\end{bmatrix} $$
$$\to\ \begin{bmatrix}
F(2n+1) \\
F(2n)
\end{bmatrix}
= \begin{bmatrix}
F(n+1) & F(n) \\
F(n) & F(n-1)
\end{bmatrix}
\begin{bmatrix}
F(n+1) & F(n) \\
F(n) & F(n-1)
\end{bmatrix}
\begin{bmatrix}
1 \\
0
\end{bmatrix} $$
$$\to\ \begin{bmatrix}
F(2n+1) \\
F(2n)
\end{bmatrix}
= \begin{bmatrix}
F(n+1)^2 + F(n)^2 \\
F(n)F(n+1) + F(n)F(n-1)
\end{bmatrix} $$
$$\to\ \begin{bmatrix}
F(2n+1) \\
F(2n)
\end{bmatrix}
= \begin{bmatrix}
F(n+1)^2 + F(n)^2 \\
F(n)[2F(n+1) - F(n)]
\end{bmatrix} $$
從以上算式可以注意到費氏數列的規則，$F(2n+1)=F(n+1)^2 + F(n)^2$ 和 $F(2n) = F(n)[2F(n+1) - F(n)]$，接著仔細觀察所有的 $2n+1$ 和 $2n$ 的費氏數列都可以由 $n$，$n+1$ 兩種數值所組成，那我們每次只需要迭代目標一半的數值，一直迭代到目標就是我們的答案。
那目前為止從上面的資訊可以推斷程式碼可以是使用迭代或是遞迴，那規則呢 ? 我們何時需要使用 $F(2n+1)=F(n+1)^2 + F(n)^2$ 或 $F(2n) = F(n)[2F(n+1) - F(n)]$，這可以從目標的二進位值來進行分析，假設我們想求 $F(8)$，那 $8$ 的二進位如下:
$$8_{10} = 1000_{2}$$
要如何從 $0000$ 變成 $1000$，首先要 $+1$ 然後左移三位
```cpp
(0000 + 1) << 3 = 1000
```
對比到費氏數列中，就可以理解到為什麼在[作業說明](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-a#-%E8%B2%BB%E6%B0%8F%E6%95%B8%E5%88%97)中，一開始 $F(0) F(1)$ 變成 $F(1) F(2)$ 一定是 $+1$ ，因為二進位的最高位要是 $1$，才會有大於 $0$ 的值。
| | $+1$ | $1\times2$ | $2\times2$ | $4\times2$ |
|:------:| ------ |:----------:|:----------:|:----------:|
| $F(0)$ | $F(1)$ | $F(2)$ | $F(4)$ | $F(8)$ |
| $F(1)$ | $F(2)$ | $F(3)$ | $F(5)$ | $F(9)$ |
以上是項數最低位為 $0$ 的情況，下面舉例項數最低位為 $1$ 的情況，以 $F(13)$ 為例
$$13_{10} = 1101_{2}$$
| | $+1$ | $1\times2+1$ | $3\times2$ | $6\times2+1$ |
|:------:| ------ |:------------:|:----------:|:------------:|
| $F(0)$ | $F(1)$ | $F(3)$ | $F(6)$ | $F(13)$ |
| $F(1)$ | $F(2)$ | $F(4)$ | $F(7)$ | $F(14)$ |
從以上情況可以歸納出兩個規則
* 遇到最低位元為 $0$ 的時候，可以直接利用上面所推導的算式直接算出 $F(2n)$ 和 $F(2n+1)$
$$F(2n) = F(n)[2F(n+1) - F(n)] \\
F(2n+1)=F(n+1)^2 + F(n)^2$$
* 遇到最低位元為 $1$ 的時候，除了要計算 $F(2n)$ 和 $F(2n+1)$ 外，還要利用費氏數列的規則去計算 $F(2n+2)$
$$F(2n) = F(n)[2F(n+1) - F(n)] \\
F(2n+1)=F(n+1)^2 + F(n)^2 \\
F(2n+2)=F(2n+1)+F(2n)$$
### 改進費氏數列
Fast doubling 程式實作:
以下程式碼將上面的想法用最直接的方式實作一次
```cpp
for (uint64_t i = 1LL << 63; i; i >>= 1) {
long long n0 = f[0] * (2 * f[1] - f[0]);
long long n1 = f[0] * f[0] + f[1] * f[1];
if (k & i) {
f[0] = n1;
f[1] = n0 + n1;
}
else {
f[0] = n0;
f[1] = n1;
}
}
```
從上圖可以看到光是只有將算法從一般的 `iterative` 換成 `fast doubling` 運行時間就降低了不少，但是秉持著任何程式碼都有可以進步的空間，所以繼續優化程式碼，上述程式碼看似沒什麼問題，但是我很想把分支拿掉。
#### 去除分支
為了把分支拿掉，思考良久，寫出適合此情況的 bitmask
```cpp
for (uint64_t i = 1LL << 63; i; i >>= 1) {
long long n0 = f[0] * (f[1] * 2 - f[0]);
long long n1 = f[0] * f[0] + f[1] * f[1];
uint64_t mask = !(k & i)-1;
f[0] = (n0 & ~mask) + (n1 & mask);
f[1] = (n0 & mask) + n1;
}
```
首先先考慮 `f[0]` 和 `f[1]` 的組合，可以從有分支的情況下看到
```cpp
if (k & i) {
f[0] = n1;
f[1] = n0 + n1;
}
else {
f[0] = n0;
f[1] = n1;
}
```
代表 `f[0]` 和 `f[1]` 是由 `n0` 和 `n1` 所組成，而 `f[0]` 當中 `n0` 和 `n1` 是交替出現，在 `f[1]` 中，`n0` 是只有 `k & i ` 不等於 0 時出現，所以可以很清楚的看出其中的 bitmask 的關係
```cpp
f[0] = (n0 & ~mask) + (n1 & mask);
f[1] = (n0 & mask) + n1;
```
而什麼樣的 bitmask 會依照 `&mask` 和 `&~mask` 去保留和棄置此項，答案就是 `1111` 和 `0000`，既然目標是利用 bitwise operation 將 `k & i` 形成 `1111` 和 `0000`，所以就開始湊可以形成這樣的算式
假設 $k = 7_{10}=0111_{2}$
```cpp
i = 1000
k & i = 0000
!(k & i) - 1 = 0000
f[0] = n0; f[1] = n1;
i >= 1 = 1100
k & i = 0100
!(k & i) - 1 = 1111
f[0] = n1; f[1] = n0 + n1;
```
這樣湊出所需要的 bitmask
從上圖看來兩個方法速度上並沒有明顯差別，可能是測試的項數不夠大又或者是有分支的方法速度本身就和 bitwise operation 差不多快，這些都還需要進一步實驗才可得知。
#### 利用 `clz` 改寫
從 gcc 中的 [`__builtin_`](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html) 提供 `__builtin_clz()` 可以用來計算出現 1 的最高位元在左邊數過來第幾位，利用此函式修改的程式碼如下
```cpp
for (uint32_t i = 1UL << (31 - __builtin_clz(k)); i; i >>= 1) {
long long n0 = f[0] * ((f[1] << 1) - f[0]);
long long n1 = f[0] * f[0] + f[1] * f[1];
uint64_t mask = !(k & i)-1;
f[0] = (n0 & ~mask) + (n1 & mask);
f[1] = (n0 & mask) + n1;
}
```
測試的結果顯示利用 `clz` 的程式碼比較慢，並且有隨著項數增加而時間越長的趨勢，所以去找此函式的程式碼了解為何會出現此情況
```cpp
int clz (int x){
if (x == 0) return 32;
int n = 0;
if ((x & 0xFFFF0000) == 0) {n += 16; x <<= 16;}
if ((x & 0xFF000000) == 0) {n += 8; x <<= 8;}
if ((x & 0xF0000000) == 0) {n += 4; x <<= 4;}
if ((x & 0xC0000000) == 0) {n += 2; x <<= 2;}
if ((x & 0x80000000) == 0) n += 1;
return n;
}
```
從程式碼中可以看到 `__builtin_clz()` 在函式內已經使用了 bitwise operation 來計算數值的 leading zeros，所以在函式內移過一次， for 迴圈開始時又移了一次，時間上會比只有在 for 迴圈移一次還多，至於隨著項數增加而時間越長的趨勢則需要去實作大數運算才能確認。
### 支援大數運算
#### 創建 `bn` 結構體
為了能夠計算超過 `long long` 位元組所限制的費氏數列，所以利用 `bn` 結構體來實作費氏數列。
**`bn` 結構體**
```cpp
typedef struct _bn {
unsigned int *number;
unsigned int size;
} bn;
```
從結構體中可以將大數視為 $2^{32}$ 進位的數字，`size` 則是此大數的大小，單位則為 32 位元，以費氏數列 $F(92) = 7540113804746346429$ 做舉例，此數字超過 32 位元必須要用 `long long` 才能儲存，以 `bn` 結構體表示為
> number[0] = 7540113804746346429 & 4294967295
> number[1] = (7540113804746346429 >> 32) & 4294967295
> size = 2
而問題來了，當有一個數值超過了 `long long` 所能表示的範圍的時候，我們要如何去表示，所以在這裡參考了 [作業說明：預先計算儲存 $F(n)$ 所需的空間](https://hackmd.io/@sysprog/linux2023-fibdrv/%2F%40sysprog%2Flinux2023-fibdrv-b#%E9%A0%90%E5%85%88%E8%A8%88%E7%AE%97%E5%84%B2%E5%AD%98-Fn-%E6%89%80%E9%9C%80%E7%9A%84%E7%A9%BA%E9%96%93) ，可以利用公式解先算出此數值的位元數
$$F(n) bit = \frac{-1160964 + 694242*n}{10^6} + 1$$
```cpp
int fibn_per_32bit(int k)
{
return k < 2 ? 1 : ((long) k * 694242 - 1160964) / (1000000) + 1;
}
```
接著是如何創建一個 `bn` 結構體，以下是參考自[作業說明](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)
```cpp
bn *bn_alloc(size_t size)
{
bn *new = kmalloc(sizeof(bn), GFP_KERNEL);
new->number = kmalloc(sizeof(int) * size, GFP_KERNEL);
memset(new->number, 0, sizeof(int) * size);
new->size = size;
return new;
}
```
有 `bn_alloc` 就要有 `bn_free`
```cpp
int bn_free(bn *src)
{
if (src == NULL)
return -1;
kfree(src->number);
kfree(src);
return 0;
}
```
#### `bn` 運算
然後觀察費氏數列需要什麼運算符號，從[最原始的算法](https://github.com/sysprog21/fibdrv/blob/master/fibdrv.c#L27)中可以看到需要加法和互換兩種運算式
**`bn_add`**
此程式在做 $2^{32}$ 進位加法，其中的運算與 $10$ 進位加法完全相同，先算低位元的加法，用 `unsigned long long int` 去儲存結果，若有進位的數值自然會在右移 32 位元之後出現。
```cpp
void bn_add(bn *a, bn *b, bn *c)
{
unsigned long long int 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 int) tmp1 + tmp2;
c->number[i] = carry;
carry >>= 32;
}
}
```
**`bn_swap`**
互換兩個 `bn` 的位置，在迭代時常會用到
```cpp
void bn_swap(bn *a, bn *b)
{
bn tmp = *a;
*a = *b;
*b = tmp;
}
```
**`bn_to_string`**
透過 ASCII code 轉換，將數字和 `'0'` 的差以 char 陣列儲存，並且從高位逐一檢查是否有進位，若有進位則處理成沒有進位的數字
```cpp
char *bn_to_string(bn src)
{
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 int d = 1U << 31; d; d >>= 1) {
/* binary -> decimal string */
int carry = !!(d & src.number[i]);
for (int j = len - 2; j >= 0; j--) {
s[j] += s[j] - '0' + carry; // double it
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;
}
```
**`bn_multSSA`**
利用 [Schönhage–Strassen Algorithm](https://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm) 來實作乘法，Schönhage–Strassen Algorithm 的主要觀念是講兩個要相乘的數字分成若干的小數字相乘後相加，而在大數中，我們已經有分好的小數字也就是 `bn` 結構體中的 `number`，所以應用在大數中的 Schönhage–Strassen Algorithm 會變成以下型式
```
[]為此大數的第幾個number
a[2] a[1] a[0]
x b[1] b[0]
------------------------------------------------------------------------
a[2]*b[0] a[1]*b[0] a[0]*b[0]
a[2]*b[1] a[1]*b[1] a[0]*b[1]
------------------------------------------------------------------------
a[2]*b[1] a[2]*b[0]+a[1]*b[1] a[1]*b[0]+a[0]*b[1] a[0]*b[0]
```
在大數 $a$ 和大數 $b$ 相乘之中，$(a_2b_1,a_2b_0+a_1b_1,a_1b_0+a_0b_1,a_0b_0)$ 為 $a$ 和 $b$ 的線性卷積，而若大數 $c = ab$，那 $c$ 當中的 $number$ 就可以對應到 $a$ 和 $b$ 的線性卷積
```
a[2]*b[1] a[2]*b[0]+a[1]*b[1] a[1]*b[0]+a[0]*b[1] a[0]*b[0]
| | | |
| | | |
c[3] c[2] c[1] c[0]
```
從上面 $c$ 與 $ab$ 的對應關係可以看出 $c[x] = \sum\limits_{x = 0 ,y = 0}^{(a \to size) - 1}a[y]b[x-y]$，其對應的程式碼則如下
```cpp
unsigned long long int carry = 0;
for (int i = 0; i < c->size; i++) {
int j = a->size;
while (j >= 0) {
unsigned int tmp1 = ((j <= i) && j < a->size) ? a->number[j] : 0;
unsigned int tmp2 = ((i - j) < b->size) ? b->number[i - j] : 0;
carry += (unsigned long long int) tmp1 * tmp2;
j--;
}
c->number[i] = carry;
carry >>= 32;
}
```
而為了不讓大數 $c$ 的大小在 fast doubling 時一直膨脹，所以在乘法開始執行之前，使用下面程式碼來限制 $c$ 的大小
```cpp
int d = bn_msb(a) + bn_msb(b);
d = DIV_ROUNDUP(d, 32) + !d; // round up, min size = 1
bn_resize(c, d);
```
在測試時發現 $F(293)$ 和 [The first 500 Fibonacci numbers](https://blog.abelotech.com/posts/first-500-fibonacci-numbers/) 當中 $F(293)$ 數字並不符合
> 293 7654090467756936378415884538348976340768064993978954512095813
```cpp
F(293) = 7654090467756936378415203973615134463841138244764090975672901
```
檢查過後發現在進行加法時，$c[2]$ 的總和超過了 unsigned long long int 的上限，所以將 unsigned long long int 改成 __uint128_t ，由此可計算超過 1000 項費氏數列
```diff
-unsigned long long int carry = 0;
+__uint128_t carry = 0;
for (int i = 0; i < c->size; i++) {
int j = a->size;
while (j >= 0) {
unsigned int tmp1 = ((j <= i) && j < a->size) ? a->number[j] : 0;
unsigned int tmp2 = ((i - j) < b->size) ? b->number[i - j] : 0;
- carry += (unsigned long long int) tmp1 * tmp2;
+ carry += (__uint128_t) tmp1 * tmp2;
j--;
}
c->number[i] = carry;
carry >>= 32;
}
```
#### `bn` 支援費氏數列運算
使用 `bn` 重新實作一次原始費氏數列運算
```cpp
void bn_fib(bn *ret, unsigned int n)
{
int nsize = fibn_per_32bit(n) / 32 + 1;
bn_resize(ret, nsize);
if (n < 2) {
ret->number[0] = n;
return;
}
bn *f1 = bn_alloc(nsize);
bn *tmp = bn_alloc(nsize);
ret->number[0] = 0;
f1->number[0] = 1;
for (unsigned int i = 1; i < n + 1; i++) {
bn_add(ret, f1, tmp);
bn_swap(f1, ret);
bn_swap(f1, tmp);
}
bn_free(f1);
bn_free(tmp);
}
```
修改 `fib_read`，使得 `fib_read` 可以回傳字串到 user space 中
```cpp
bn *fib = bn_alloc(1);
bn_fib(fib, *offset);
char *p = bn_to_string(*fib);
copy_to_user(buf, p, len);
bn_free(fib);
```
調整測試上限，便可獲得超過原本 $F(92)$ 的數值，測試結果如下
```cpp
Reading from /dev/fibonacci at offset 493, returned the sequence 4801994309516818408452995190383973646671835537213025106017041525333156522800379742496995285687643305513.
Reading from /dev/fibonacci at offset 494, returned the sequence 7769790006581794836315417026718114982822736329999469724305586980435409187727711750121668968517028834617.
Reading from /dev/fibonacci at offset 495, returned the sequence 12571784316098613244768412217102088629494571867212494830322628505768565710528091492618664254204672140130.
Reading from /dev/fibonacci at offset 496, returned the sequence 20341574322680408081083829243820203612317308197211964554628215486203974898255803242740333222721700974747.
Reading from /dev/fibonacci at offset 497, returned the sequence 32913358638779021325852241460922292241811880064424459384950843991972540608783894735358997476926373114877.
Reading from /dev/fibonacci at offset 498, returned the sequence 53254932961459429406936070704742495854129188261636423939579059478176515507039697978099330699648074089624.
Reading from /dev/fibonacci at offset 499, returned the sequence 86168291600238450732788312165664788095941068326060883324529903470149056115823592713458328176574447204501.
Reading from /dev/fibonacci at offset 500, returned the sequence 139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125
```
甚至可以測到 $F(1000)$ 的費氏數列
```cpp
Reading from /dev/fibonacci at offset 993, returned the sequence 1497068822810020534399232227533005158498637007856137201747455252741318759630046838167766787649593980126381385456182553867019240362815567640483762190938317523078461741297449202133947902177195835489685370439138.
Reading from /dev/fibonacci at offset 994, returned the sequence 2422308238804407089268820758059737271890428034963283291423237466738923487762205937312441964529423332944990116110066323372235058978516564015277004931960761593992730308301490464065917906637454573375206452747367.
Reading from /dev/fibonacci at offset 995, returned the sequence 3919377061614427623668052985592742430389065042819420493170692719480242247392252775480208752179017313071371501566248877239254299341332131655760767122899079117071192049598939666199865808814650408864891823186505.
Reading from /dev/fibonacci at offset 996, returned the sequence 6341685300418834712936873743652479702279493077782703784593930186219165735154458712792650716708440646016361617676315200611489358319848695671037772054859840711063922357900430130265783715452104982240098275933872.
Reading from /dev/fibonacci at offset 997, returned the sequence 10261062362033262336604926729245222132668558120602124277764622905699407982546711488272859468887457959087733119242564077850743657661180827326798539177758919828135114407499369796465649524266755391104990099120377.
Reading from /dev/fibonacci at offset 998, returned the sequence 16602747662452097049541800472897701834948051198384828062358553091918573717701170201065510185595898605104094736918879278462233015981029522997836311232618760539199036765399799926731433239718860373345088375054249.
Reading from /dev/fibonacci at offset 999, returned the sequence 26863810024485359386146727202142923967616609318986952340123175997617981700247881689338369654483356564191827856161443356312976673642210350324634850410377680367334151172899169723197082763985615764450078474174626.
Reading from /dev/fibonacci at offset 1000, returned the sequence 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875.
```
將上面的算法跑 50 次後算平均可以發現不管是在 user space 還是 kernal to user 所花的時間都大到誇張，表示在程式在處理 `copy_to_user` 時花了大量的時間，而上圖中 kernel 看似沒有太大的變化，但當我們單獨拉出 kernel 運行時間的圖
就可以發現，kernel 運行時間也是隨著費氏數列的增加而變多
**實作大數 fast doubling**
```cpp
void bn_fastdoub_fib(bn *ret, unsigned int n)
{
int nsize = fibn_per_32bit(n) / 32 + 1;
bn_resize(ret, nsize);
if (n < 2) {
ret->number[0] = n;
return;
}
bn_reset(ret);
bn *f1 = bn_alloc(nsize);
f1->number[0] = 1;
bn *n0 = bn_alloc(nsize);
bn *n1 = bn_alloc(nsize);
for (uint32_t i = 1UL << (31 - __builtin_clz(n)); i; i >>= 1) {
bn *tmp1 = bn_alloc(1);
bn *tmp2 = bn_alloc(1);
bn *tmp3 = bn_alloc(1);
/*n0 = f[0] * (f[1] * 2 - f[0]);*/
bn_add(f1, f1, n0);
bn_sub(n0, ret, n0);
bn_multSSA(n0, ret, tmp1);
/*n1 = f[0] * f[0] + f[1] * f[1];*/
bn_multSSA(ret, ret, tmp2);
bn_multSSA(f1, f1, tmp3);
bn_add(tmp2, tmp3, n1);
if (n & i) {
bn_cpy(ret, n1);
bn_cpy(f1, tmp1);
bn_add(f1, n1, f1);
} else {
bn_cpy(ret, tmp1);
bn_cpy(f1, n1);
}
}
bn_free(n0);
bn_free(n1);
bn_free(f1);
}
```
從上面兩張圖可以看到相比於原始費氏數列運算速度快上許多，但是在 kernel space 的運行速度和 user space 的運行速度相比還是相差太多，所以開啟 perf 檢查到底是什麼造成 user space 的時間增加
```shell
$ sudo perf record -g --call-graph dwarf ./test_time
$ sudo perf report --stdio -g graph,0.5,caller
96.86% 0.00% test_time libc-2.31.so [.] __GI___read (inlined)
vfs_read
|
--96.49%--fib_read
|
|--92.59%--bn_to_string
|
--3.46%--bn_fastdoub_fib
|
|--1.51%--bn_multSSA
| |
| --0.51%--bn_resize
|
--0.88%--bn_alloc
```
Sign in with Wallet
Connect another wallet