# 2022q1 Homework3 (fibdrv)
contributed by < `jim12312321` >
> [作業說明](https://hackmd.io/@sysprog/linux2022-fibdrv)
## 實驗環境
```shell
$ gcc --version
gcc (Ubuntu 11.2.0-7ubuntu2) 11.2.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): 20
On-line CPU(s) list: 0-19
Thread(s) per core: 1
Core(s) per socket: 14
Socket(s): 1
NUMA node(s): 1
Vendor ID: GenuineIntel
CPU family: 6
Model: 154
Model name: 12th Gen Intel(R) Core(TM) i7-12700H
Stepping: 3
CPU MHz: 2700.000
CPU max MHz: 6000.0000
CPU min MHz: 400.0000
BogoMIPS: 5376.00
Virtualization: VT-x
L1d cache: 336 KiB
L1i cache: 224 KiB
L2 cache: 8.8 MiB
NUMA node0 CPU(s): 0-19
```
## 自我檢查清單
- [x] 研讀上述 ==Linux 效能分析的提示== 描述，在自己的實體電腦運作 GNU/Linux，做好必要的設定和準備工作 $\to$ 從中也該理解為何不希望在虛擬機器中進行實驗;
- [x] 研讀上述費氏數列相關材料 (包含論文)，摘錄關鍵手法，並思考 [clz / ctz](https://en.wikipedia.org/wiki/Find_first_set) 一類的指令對 Fibonacci 數運算的幫助。請列出關鍵程式碼並解說
- [ ] 複習 C 語言 [數值系統](https://hackmd.io/@sysprog/c-numerics) 和 [bitwise operation](https://hackmd.io/@sysprog/c-bitwise)，思考 Fibonacci 數快速計算演算法的實作中如何減少乘法運算的成本;
- [ ] 研讀 [KYG-yaya573142 的報告](https://hackmd.io/@KYWeng/rkGdultSU)，指出針對大數運算，有哪些加速運算和縮減記憶體操作成本的舉措？
- [ ] `lsmod` 的輸出結果有一欄名為 `Used by`，這是 "each module's use count and a list of referring modules"，但如何實作出來呢？模組間的相依性和實際使用次數 ([reference counting](https://en.wikipedia.org/wiki/Reference_counting)) 在 Linux 核心如何追蹤呢？
> 搭配閱讀 [The Linux driver implementer’s API guide » Driver Basics](https://www.kernel.org/doc/html/latest/driver-api/basics.html)
- [x] 注意到 `fibdrv.c` 存在著 `DEFINE_MUTEX`, `mutex_trylock`, `mutex_init`, `mutex_unlock`, `mutex_destroy` 等字樣，什麼場景中會需要呢？撰寫多執行緒的 userspace 程式來測試，觀察 Linux 核心模組若沒用到 mutex，到底會發生什麼問題。嘗試撰寫使用 [POSIX Thread](https://en.wikipedia.org/wiki/POSIX_Threads) 的程式碼來確認。
---
## 前置準備工作
:::spoiler 排除干擾效能的因素
- 孤立特定的 cpu
- 在 /etc/default/grub 中用 sudo 權限編輯，並在 `GRUB_CMDLINE_LINUX` 後加入 `isolcpus=0` ，以空出第 0 個 cpu
```
GRUB_CMDLINE_LINUX="isolcpus=0"
```
- 更新 grub (也要使用 sudo 權限)
```
＄sudo update-grub
```
- 重開機
- 大功告成！
```
pid 1's current affinity list: 1-19
```
- 抑制 address space layout randomization
```
$ sudo sh -c "echo 0 > /proc/sys/kernel/randomize_va_space"
```
- 設定 scaling_governor 為 performance
- 寫腳本，檔名為 `performance.sh`
```
for i in /sys/devices/system/cpu/cpu*/cpufreq/scaling_governor
do
echo performance > ${i}
done
```
- 執行腳本
```
$ sudo sh performance.sh
```
- 關閉 turbo mode
```
$ sudo sh -c "echo 1 > /sys/devices/system/cpu/intel_pstate/no_turbo"
```
- 把上述設定全部整合在 `performance.sh` 中，便於後續重新執行。
:::
---
## 初探 fibdrv
- 將 [K04: fibdrv](https://hackmd.io/@sysprog/linux2022-fibdrv) 中提到可以測量時間的方法加進 fibdrv.c
- 需注意不須使用 `DEFINE_KTIME` 進行初始化，於這個 [PATCH](https://github.com/spotify/linux/commit/272705c5979c114e63dbfcd28ea15093038a4c42) 時，`DEFINE_KTIME` 已被移除。
- 修改 Makefile 會造成錯誤的指令刪除
```
@diff -u out scripts/expected.txt
```
- 這道指令雖然能讓我們看到輸出檔案跟 `scripts/expected.txt` 的差別，但是 `scripts/expected.txt` 中的資料也只是完全沒修改時會印出的內容罷了，留著這到指令會造成`make check` 時錯誤(我不知道為什麼，不使用 `make check` 單獨使用這道指令時又沒問題！)且會讓版面很亂，故刪除。
- 將 `./client` 指定在特定的 CPU 運行並紀錄結果
```
$ sudo taskset -c 0 ./client >./output/output.txt
```
- 在 `client.c` 中加入對 userspace time 的測量
- 耗費時間散佈圖
---
## 摘錄 Fibonacci 數關鍵手法，並思考 clz / ctz 一類的指令對 Fibonacci 數運算的幫助
### 費式數列 - 基本定義
$F_n=\left\{
\begin{aligned}
0, n=0 \\
1, n=1 \\
F_{n-2}+F_{n-1},n>1
\end{aligned}
\right.$
因此由費式數列的基本定義可以得知，除了前兩項，費式數列的算法為 $F_n = F_{n-2}+F_{n-1}$ 。
### fast Doubling
更詳細的推導在 [K04: fibdrv](https://hackmd.io/@sysprog/linux2022-fibdrv#-%E8%B2%BB%E6%B0%8F%E6%95%B8%E5%88%97) 中有說明，這邊就不重複，僅列出結論。
$F(2k)=F(k)[2F(k+1)-F(k)]$
$F(2k+1)=F(k+1)^2+F(k)^2$
這表示我們可以用 $F(k)$ 和 $F(k+1)$ 求出 $F(2k)$ 和 $F(2k+1)$ ，舉例來說：
>用 fast doubling 求出 F(10):
>此時我們可知 $k=5$ ，表示 F(10) 可透過 F(5) 和 F(6) 求出，而 F(5) 和 F(6) 又分別可透過 F(2),F(3) 和 F(3),F(4) 求出...，以此類推，最後我們可以得知要求出 F(10) ，可透過 F(0),F(1),F(2),F(3),F(4),F(5),F(6) 求出答案。
>比起原先定義中從 F(0) 一路算到 F(8),F(9) 才能求出 F(10) 明顯更快。
所以費式數列的定義可進一步變成：
$F_n=\left\{
\begin{aligned}
0, n=0 \\
1, n=1 \\
F(k)[2F(k+1)-F(k)],n\ is\ even\ and\ k = \dfrac{n}{2} \\
F(k+1)^2+F(k)^2,n\ is\ odd\ and\ k = \lfloor\dfrac{n}{2}\rfloor
\end{aligned}
\right.$
### 考量到硬體加速的 $F(n)$
>1. 省略 $F(0)$，直接從 $F(1)$ 開始;
>2. clz/[ffs](https://www.kernel.org/doc/htmldocs/kernel-api/API-ffs.html): 先去除數字 MSB 起算的開頭 0 位元，因為真正的數字不包含 leading 0s，所以計算它也沒用，因此需要 [clz](https://en.wikipedia.org/wiki/Find_first_set#CLZ) 計算有多少 leading 0s
* 遇到 `0` $\rightarrow$ 進行 fast doubling，也就是求 $F(2n)$ 和 $F(2n+1)$
* 遇到 `1` $\rightarrow$ 進行 fast doubling，也就是先求 $F(2n)$ 和 $F(2n+1)$，再求 $F(2n+2)$
關於 clz / ctz 一類的指令如何在運算費式數列中加速，如同上面 [K04: fibdrv](https://hackmd.io/@sysprog/linux2022-fibdrv#-%E8%B2%BB%E6%B0%8F%E6%95%B8%E5%88%97) 中的說明，可以先去除 leading 0s ，直接從 first set bit 開始算起。
- 耗費時間圖
- n $\rightarrow$ normal, f $\rightarrow$ fast doubling
- k $\rightarrow$ kernel space, u $\rightarrow$ user space
從前 92 個的計算結果看來，有無進行 fast doubling 的差異不大。
:::spoiler normal fib seq 程式碼
```c=
static long long fib_sequence(long long k)
{
/* a before b */
long long a,b;
a = 0;
b = 1;
for (int i = 1; i <= k; i++) {
a += b;
_swap(a,b);
}
return a;
}
```
:::
:::spoiler fast doubling fib seq 程式碼
```c=
static long long fast_fib_seq(long long k)
{
if(k==0)
return 0;
/*a before b*/
long long a,b,tmp1,tmp2;
a = 0;
b = 1;
int idx = __builtin_clzll(k)+1;
for(;idx<=MAX_LL;idx++){
tmp1 = a*(2*b-a);
tmp2 = b*b + a*a;
a = tmp1;
b = tmp2;
if(_get_bit(k,(MAX_LL-idx))){ // current bit == 1
tmp1 = a+b;
a = b;
b = tmp1;
}
}
return a;
}
```
:::
---
## 減少乘法運算的成本
以下是 fast doubling 的虛擬碼 - from [K04: fibdrv](https://hackmd.io/@sysprog/linux2022-fibdrv)
```c
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;
```
### 嘗試改寫平方運算
可以發現其中大部分的乘法運用都是在做平方運算，所以如果能針對平方運算做優化，則理論上能減少整體運算成本。
```c=
#define MAX_INT 32
#define SET_BIT(x,n) (x | 1 << n)
#define CHECK_BIT(x,n) (x >> n & 1)
/* return a*a
* main idea is using 多項式平方.
* (a+b+c)^2 = a^2 + b^2 + c^2 + 2ab + 2ac + 2bc.
*/
int square_opt(int a){
int ret = 0;
// idx == int index of first set
int fs_idx = MAX_INT -__builtin_clz(a)-1;
int cnt = fs_idx;
do{
if(CHECK_BIT(a,cnt)){
ret += SET_BIT(0,cnt<<cnt);
for(int cur = cnt+1;cur <= fs_idx;cur++){ // handle 2ab
if(CHECK_BIT(a,cur))
ret += SET_BIT(0,cnt+cur+1); // 2ab
}
}
}while(--cnt >= 0);
return ret;
}
```
簡單測試 (使用 [Online C Compiler](https://www.onlinegdb.com/online_c_compiler) 和 `<time.h>` 進行測量)
```
using operator "*" :
5656 * 5656 = 31990336
normal square : 100 ns.
using square_opt:
5656 * 5656 = 31990336
opt square : 472 ns.
```
目前的結果看來是沒有優化成功的，反而成果越糟 :cry: 繼續思考有沒有改進空間。
於是我換了一個思路，利用將十進位換成二進位 (e.g. $5=2^2+2^0$)，然後再利用乘法分配律相加。
```c=
/*
* a * a
* = a * (2^n + ... + 2 ^k)
* = a << n + ... + a << k
*/
int square_opt_v2(int a){
int ret = 0;
int cnt = MAX_INT-__builtin_clz(a)-1;
do{
if(CHECK_BIT(a,cnt)){
ret += a << cnt;
}
}while(--cnt >= 0);
return ret;
}
```
而經過和上面一樣的測試環境和方法，以下是得到的結果。
```
using operator "*" :
5656 * 5656 = 31990336
cost time :108 ns
using square_opt_v2 :
5656 * 5656 = 31990336
cost time :167 ns
```
可以看到成果更接近單純使用 `*` 進行次方運算的速度了，但整體來說還是略慢。
我覺得 v2 的效能損耗主要在第 11 行中的 `if(CHECK_BIT(a,cnt))` ，如果能改成 branchless 應該有機會進一步增進效能，因此基於 v2 的 v3 (branchless 的 v2) 如下：
```c=
#define MASK(x,n) (-CHECK_BIT(x,n))
int square_opt_v3(int a){
int ret = 0;
int cnt = MAX_INT-__builtin_clz(a)-1;
do{
ret += (a << cnt) & MASK(a,cnt);
}while(--cnt >= 0);
return ret;
}
```
利用 `CHECK_BIT(x,n)` 的輸出只有 0 或 1 的特性，製作 MASK 達成 branchless 。
而以下是當前的測試成果：
```c=
using operator "*" :
5656 * 5656 = 31990336
cost time :166 ns
using square_opt_v3 :
5656 * 5656 = 31990336
cost time :153 ns
```
終於有機會比用 `*` 還快了 :tada:
接著需要經過測試才可以知道是不是穩定的表現更好。
- `*` 與 `square_opt_v3` 差異圖
其中 diff 為調整前減去調整後，意即值大於 0 表示有改進，反之小於 0 表示比原先表現更差。
> 計算 123456 ^ 2 1000 次並紀錄每次的 diff
可以從目前的測試結果看到， v3 依舊無法比原先的乘法算的更快，也沒有任何測試中能得到比原先表現更好的數據。
### 嘗試改寫乘法運算
看著程式碼發呆一下後，突然發現其實我在實作的 `square_opt` 轉個彎後其實就是在實作二進位的乘法！因此我稍微改一下，附上程式跟測試結果圖。
- 程式碼
```c=
/* return a * b */
long long mul_opt(long long a,long long b){
long long ret = 0;
long long cnt = MAX_LL-__builtin_clzll(a)-1;
do{
ret += (b << cnt) & MASK(a,cnt);
}while(--cnt >= 0);
return ret;
}
```
- `*` 與 `mul_opt` 差異圖
> 計算 123456 * 654321 1000 次並紀錄每次的 diff
還是無法增進效能，甚至也還無法逼近原先的效能。
---
## 擴增費式數列的可運算範圍 <br>- 加法
### 實作成果
- 可運算最大值
目前到 500-th 都是正確的
```
sudo taskset -c 0 ./client >./output/output.txt
...
Reading from /dev/fibonacci at offset 499, returned the sequence 86168291600238450732788312165664788095941068326060883324529903470149056115823592713458328176574447204501. cost time in kernel: 1682468 ns,cost time in userspace: 1682894 ns.
Reading from /dev/fibonacci at offset 500, returned the sequence 139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125. cost time in kernel: 1691187 ns,cost time in userspace: 1691595 ns.
...
```
- 耗費時間散佈圖 （log）
### 發想與細節
實作程式碼 ： [My fibdrv](https://github.com/jim12312321/fibdrv)
在看到[作業說明](https://hackmd.io/@sysprog/linux2022-fibdrv#%E8%A8%88%E7%AE%97-F93-%E5%8C%85%E5%90%AB-%E4%B9%8B%E5%BE%8C%E7%9A%84-Fibonacci-%E6%95%B8---%E4%BD%BF%E7%94%A8%E6%95%B8%E5%AD%97%E5%AD%97%E4%B8%B2%E4%B8%A6%E5%A5%97%E7%94%A8-quiz2-SSO-Small-String-Optimization)後，覺得[大數實作 （by AdrianHuang）](https://github.com/AdrianHuang/fibdrv/tree/big_number) 的作法很有意思因此想自己嘗試做做看利用數字字串相加的功能。
- 結構體
```cpp=
typedef struct fib_node {
char data[256];
} fib_node;
```
就是個單純拿來放數字的結構（現在想想好像用 char 的二維陣列就能實作了）。
當初看到[大數實作 （by AdrianHuang）](https://github.com/AdrianHuang/fibdrv/tree/big_number) 中的結構體 `xs` 很複雜，想簡化成簡單明瞭的形式。本結構體一樣能達到擴增到費氏數列第 500 個的功能，因此我不太清楚當初 `xs` 為什麼要這樣設計。
- 字串加法
```cpp=
#include <linux/string.h>
/*
* add string a and b mathematically.
* add the char in a and b backward,and record carry if necessary.
* reverse the answer in the end.
*/
static void string_add(char *a, char *b, char *out)
{
short tmp = 0, carry = 0;
long len_a = strlen(a), len_b = strlen(b);
long len_max = _max(len_a, len_b);
for (int i = 0; i < len_max; i++) {
// cppcheck-suppress shiftTooManyBitsSigned
tmp = ((len_a - i - 1) >> 63) ? 0 : (a[len_a - i - 1] - '0');
// cppcheck-suppress shiftTooManyBitsSigned
tmp += ((len_b - i - 1) >> 63) ? 0 : (b[len_b - i - 1] - '0');
if (carry)
tmp += carry;
if (tmp >= 10) {
tmp -= 10;
carry = 1;
} else {
carry = 0;
}
out[i] = tmp + '0';
}
if (carry) {
out[len_max] = carry + '0';
len_max++;
}
out[len_max] = '\0';
_reverse(out);
}
```
發想︰
看到[大數實作 （by AdrianHuang）](https://github.com/AdrianHuang/fibdrv/tree/big_number) 的實作中要 reverse 很多次且要確保 a 恆大於 b ，我想要把這些步驟省下來，因此直接從 a 和 b 的尾端加回來，最後再將輸出 reverse 即可，且透過 `_max` 讓 a 不一定要恆大於 b 也能運算，其他邏輯跟[大數實作 （by AdrianHuang）](https://github.com/AdrianHuang/fibdrv/tree/big_number) 相異不大。
實作細節︰
1. 將 a 和 b 從後往前逐字元相加
2. 將得出字串反轉得到答案
---
## 擴增費式數列的可運算範圍 <br>- fast doubling
### 實作成果
- 可運算最大值
目前到 500-th 都是正確的
```
sudo taskset -c 0 ./client >./output/output.txt
...
Reading from /dev/fibonacci at offset 499, returned the sequence 86168291600238450732788312165664788095941068326060883324529903470149056115823592713458328176574447204501. cost time in kernel: 619861 ns,cost time in userspace: 620154 ns.
Reading from /dev/fibonacci at offset 500, returned the sequence 139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125. cost time in kernel: 609854 ns,cost time in userspace: 610124 ns.
...
```
- 耗費時間散佈圖 （log）
- 與用加法實作費氏數列的差異 （log）
可以看到利用 fast doubling 後越是後面的費氏數列運算的越快。
### 發想與細節
實作程式碼 ： [My fast fibdrv](https://github.com/jim12312321/fibdrv/blob/master/fast_fibdrv.c)
在 fast doubling 實作中，實作邏輯就按照[作業說明](https://hackmd.io/@sysprog/linux2022-fibdrv#-%E8%B2%BB%E6%B0%8F%E6%95%B8%E5%88%97)中提供的虛擬碼。
- 十進位轉二進位
在 fast doubling 中需要將 `k` 轉成 2 進位後從高位逐位判斷，在本實作中用 stack 儲存十進位轉二位的位數。
```c=180
/*
* convert decimal to binary.
* store the results in stack.
*/
static void d2b(long long a)
{
push(a & 1);
a >>= 1;
if (a > 0)
d2b(a);
}
```
- 字串減法
```c=115
/*
* out = a - b
* Make sure that a is always greater than b.
*/
static void string_sub(char *a, char *b, char *out)
{
// cppcheck-suppress variableScope
short tmp = 0, borrow = 0;
for (int i = 0; i < strlen(a); i++) {
// cppcheck-suppress shiftTooManyBitsSigned
tmp = ((strlen(a) - i - 1) >> 63) ? 0 : (a[strlen(a) - i - 1] - '0');
if (borrow < 0)
tmp -= 1;
// cppcheck-suppress shiftTooManyBitsSigned
tmp -= ((strlen(b) - i - 1) >> 63) ? 0 : (b[strlen(b) - i - 1] - '0');
if (tmp < 0) {
tmp += 10;
borrow = -1;
} else {
borrow = 0;
}
if ((i + 1 == strlen(a)) && tmp == 0)
break;
out[i] = tmp + '0';
}
if (tmp != 0) {
out[strlen(a)] = '\0';
} else {
out[strlen(a) - 1] = '\0';
}
_reverse(out);
}
```
實作方式基本跟字串加法相同，與加法不同的是在加法中要紀錄進位 `carry` ，而在減法中要紀錄借位 `borrow` 。
- 字串乘法
```c=148
/*
* Make sure out have been initialized with 0 before.
*/
static void string_mul(char *a, char *b, char *out)
{
short tmp = 0, carry = 0;
for (int ib = 0; ib < strlen(b); ib++) {
carry = 0;
for (int ia = 0; ia < strlen(a); ia++) {
tmp = b[strlen(b) - ib - 1] - '0';
tmp *= a[strlen(a) - ia - 1] - '0';
tmp += out[ia + ib] - '0';
if (carry)
tmp += carry;
if (tmp >= 10) {
carry = tmp / 10;
tmp %= 10;
} else {
carry = 0;
}
out[ia + ib] = tmp + '0';
}
out[ib + strlen(a)] = carry + '0';
}
/*len(a)-1+len(b)-1+1*/
tmp = strlen(a) + strlen(b) - 1;
if (carry)
tmp++;
out[tmp] = '\0';
_reverse(out);
}
```
整體作法就是直式乘法的思維，須注意 `out` 要先將裡面的值初始為全為 0 的字串，這樣在抓取先前預算過的值才不會出錯。 （`tmp += out[ia + ib] - '0';`）
---
## 大數運算的加速策略與運算邏輯
### [sysprog21/bignum](https://github.com/sysprog21/bignum)
在老師實作的大數運算中主要由 apm 和 bn 所組成，其中 apm 負責處理高精度運算而 bn 則是利用 apm 中提供的各項基礎運算，包裝成更高階的 API 以利使用。
#### Arbitrary-Precision arithmetic (apm)
高精度運算是利用數字陣列來模擬大數運算，舉例來說：
> uint8_t 的 MAX 是 255 ，那如果要表達 300 怎麼辦？
> 1. 令一個 `uint8_t *digit` 的數字陣列來儲存(同時也可以指定此陣列的 size，本例中假設 size = 3)
> 2. 已知 $300 = 0\times256^2+1\times256^1+44\times256^0$
> 3. 則 300 可以用 digit[2] = 0,digit[1] = 1,digit[0] = 44 表示
>
> 視覺化表示就是
>\begin{array}{rrrr}
> 300: & 00000000 & 00000001 & 00101100 \\
> digit: & digit[2] & digit[1] & digit[0] \\
> value: & 0 & 1 & 44 \\
>\end{array}
>最後要輸出給使用者看的時候要再經過格式轉換，也就是 [format.c](https://github.com/sysprog21/bignum/blob/master/format.c) 在做的事。
- 加法/減法
在加減法中，運算方法其實與正常數值運算雷同，根據指定的 size 從低位元組運算至高位元組並記錄 carry/borrow 。
以加法為例：
> size = 1
> 240 + 50 = 290
>
> \begin{array}{rrr}
> 240 + 50: & & \\
> & 1111 & 0000 \\
> + & 0011 & 0010 \\ \hline
> 290: & 0010 & 0010 \\
> \end{array}
>
> carry = 1
- 乘法
演算法採用 Karatsuba algorithm ，要運算 $U*V$ 時，首先會先將 U 和 V 視為 $U = U_1*2^N+U_0$, $V = V_1*2^N+V_0$，再進行乘法運算。
$U*V = (2^{2N}+2^N)U_1V_1+2^N(U_1-U_0)(V_0-V_1)+(2^N+1)U_0V_0$
而在 APM 中數值的表達方式，本身就符合 $U = U_1*2^N+U_0$ 的形式。
>舉例來說：
>假設 4 個 bits 為一組，則
>\begin{array}{rrr}
> 28 & & \\
> = & \ \ 0001 &\ \ 1100 \\
> = & 1 * 2^4 & +12 \\
>\end{array}
在乘法中首先會設定一個閾值 (`KARATSUBA_MUL_THRESHOLD`)，若 $\lfloor\dfrac{size}{2}\rfloor$ 小於閾值則執行一般直式乘法把 $U_1V_1$, $U_0V_0$ 和 $(U_1-U_0)(V_0-V_1)$ 分別算出來並放在對應的位置上進行加減法運算，即可得到兩數乘積，否則繼續依據閾值拆分。
> 舉例來說：
> 假設 4 個 bits 為一組且 `KARATSUBA_MUL_THRESHOLD` 為 1
> 要運算 $7*28$，則
> $size = 2$,$\ \ 7=0*2^4+7$,$\ \ 28=1*2^4+12$
>$\rightarrow$ $U_1V_1 = 0$,$\ \ U_0V_0 = 84$, $\ \ (U_1-U_0)(V_0-V_1) = -77$
>(為求簡潔，以下皆用 10 進位取代 2 進位，$0010\Rightarrow/2$)
>\begin{array}{rrrr|r}
> 7 * 28 & & & & \\
> & & /5 & /4 & 84 \\
> + & /5 & /4 & & (2^4+1) * 84 \\
> - & /4 & /13 & & - 2^4 * 77 \\ \hline
> & /0 & /12 & /4 & 196 \\
>\end{array}
>\begin{array}{rr}
> 196 & & \\
> = & \ \ 1100 &\ \ 0010 \\
> = & 12 * 2^4 & +4 \\
>\end{array}
上述例子中 28 和 7 分別被拆成兩項，若是進行正常的大數乘法(逐 digit 相乘相加)，則需要 $2^2$，考慮兩個大數且都被需拆 $N$ 次，則需要 $(N+1)^2$ 次乘法運算。
而使用 Karatsuba algorithm 則會需要 $3*N$ 次乘法運算。
因此使用 Karatsuba algorithm 對於大數運算能減少許多乘法運算帶來的成本。
#### bn
在 bn 中進一步組合 APM 中寫好的高精度運算函式成為面向使用者的 API 。
除此之外，也將在 apm 運算中常用的 digit, size 等封裝成 bn 結構體。
```c=
typedef struct {
apm_digit *digits; /* Digits of number. */
apm_size size; /* Length of number. */
apm_size alloc; /* Size of allocation. */
unsigned sign : 1; /* Sign bit. */
} bn, bn_t[1];
```
### [KYG-yaya573142 的報告](https://hackmd.io/@KYWeng/rkGdultSU)
#### 加速運算
- 改寫 fast donbling 的運算式
將原先老師提供的 fast doubling (簡化中間過程)，以下簡稱為 ver1：
:::info
$\begin{split}
\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}\\ \\ &=
\begin{bmatrix}
F(n+1)^2 + F(n)^2\\
F(n)F(n+1) + F(n-1)F(n)
\end{bmatrix}
\end{split}$
$\Downarrow$
\begin{split}
F(2k) & =F(k)[2F(k+1)-F(k)] \\
F(2k+1) & =F(k+1)^2+F(k)^2 \\
\end{split}
:::
使用不同的 Q-Matrix 改寫成 (簡化中間過程)，以下簡稱為 ver2：
:::info
$\begin{split}
\begin{bmatrix}
F(2n-1) \\
F(2n)
\end{bmatrix} &=
\begin{bmatrix}
0 & 1 \\
1 & 1
\end{bmatrix}^{2n}
\begin{bmatrix}
F(0) \\
F(1)
\end{bmatrix}\\ \\ &=
\begin{bmatrix}
F(n)^2 + F(n-1)^2\\
F(n)F(n) + F(n)F(n-1)
\end{bmatrix}
\end{split}$
$\Downarrow$
\begin{split}
F(2k-1) & =F(k)^2 + F(k-1)^2 \\
F(2k) & =F(k)[2F(k-1) + F(k)] \\
\end{split}
:::
在 [KYG-yaya573142 的報告](https://hackmd.io/@KYWeng/rkGdultSU)聲稱 `ver2` 可以少掉一次迴圈及避免使用減法，避免使用減法這點很好理解，`ver2` 的確沒有減法，又因為在 [bn_sub](https://github.com/KYG-yaya573142/fibdrv/blob/optimize-bn/bn_kernel.c#L296) 的實現方法為反轉 sign bit 達成減法效果，因此減少減法，的確會達到優化的目的。
而針對少掉一次迴圈的部分我存有一些疑惑，首先為了在同一個基準上比較，首先應該將 ver1 原先的實現方法也加上 `__bulitin_clz` ，以下為簡易測試函式。
```c=
unsigned int fib_fdoubling_ver1(unsigned int n)
{
if (n < 2) { //Fib(0) = 0, Fib(1) = 1
return n;
}
unsigned int f1 = 0; //F(k)
unsigned int f2 = 1; //F(k+1)
unsigned int k1;
unsigned int k2;
for (unsigned int i = 1U << (31 - __builtin_clz(n)); i; i >>= 1) {
/* F(2k) = F(k) * [ 2 * F(k+1) – F(k) ] */
k1 = f1 * (2 * f2 - f1);
/* F(2k+1) = F(k)^2 + F(k+1)^2 */
k2 = f1*f1 + f2*f2;
if (n & i) {
f1 = k2;
f2 = k1;
f2 += k2;
} else {
f1 = k1;
f2 = k2;
}
}
return f1;
}
unsigned int fib_fdoubling_ver2(unsigned int n)
{
if (n < 2) { //Fib(0) = 0, Fib(1) = 1
return n;
}
unsigned int f1 = 0; // F(k-1)
unsigned int f2 = 1; // F(k)
unsigned int k1;
unsigned int k2;
for (unsigned int i = 1U << (30 - __builtin_clz(n)); i; i >>= 1) {
/* F(2k-1) = F(k)^2 + F(k-1)^2 */
k1 = f2*f2 + f1*f1;
/* F(2k) = F(k)[2 * F(k-1) + F(k)] */
k2 = f2 * (2 * f1 + f2);
if (n & i) {
f1 = k2;
f2 = k1;
f2 += k2;
} else {
f1 = k1;
f2 = k2;
}
}
return f2;
}
```
其中關鍵的原因在於 `ver1` 的 $F(2k)$ 更新方式需要透過 $F(k)*value$ ，因此必須從 n 的 last set bit 開始迴圈，才可運算正確的結果，若一樣用 `1U << (30 - __builtin_clz(n))` ，則會因為一開始 $F(k)$ 為 0 進而無法正確更迭費式數列的值。
若要使用與 `ver2` 相同的迴圈起點， `1U << (30 - __builtin_clz(n))` (換句話說就是從 n 的 last set bit 的下一個 bit 開始)，則就必須滿足以下任一條件：
1. 由於跳過 last set bit ，因此必須更改 `f1` 與 `f2` 的值
2. 推導其他運算 fast doubling 的算式
在 `ver2` 中利用改變 Q-Matrix 使得運算 $F(2k)$ 時不會受到 $F(k)$ 為 0 的影響(條件 (2) 的解法)。
但從條件 (1) 著手一樣可以正確處理費式數列的運算，已知由於 $F(k)$ 為 0 會造成運算錯誤，因此只要將其改成 1 後續即可正確運算 (如同以下的 `ver1_1`)。
```diff=
+ unsigned int fib_fdoubling_ver1_1(unsigned int n)
- unsigned int fib_fdoubling_ver1(unsigned int n)
{
...
+ unsigned int f1 = 1;
- unsigned int f1 = 0;
...
+ for (unsigned int i = 1U << (30 - __builtin_clz(n)); i; i >>= 1) {
- for (unsigned int i = 1U << (31 - __builtin_clz(n)); i; i >>= 1) {
...
}
```
另外 [bn_fib_fdoubling](https://github.com/KYG-yaya573142/fibdrv/blob/optimize-bn/bn_kernel.c#L379) 中本來就有對 n == 1 or 0 的情況做特殊處理，所以更改 `f1` 的值也不用擔心 n == 0 時會有錯誤的回傳值。
```c=
void bn_fib_fdoubling(bn *dest, unsigned int n)
{
...
if (n < 2) { // Fib(0) = 0, Fib(1) = 1
dest->number[0] = n;
return;
}
...
}
```
最後，我認為在 [KYG-yaya573142 的報告](https://hackmd.io/@KYWeng/rkGdultSU)造成時間小幅度減少 3% 的原因，主要是來自引入 `ver2` 的算法後，必須搭配使用 `clz` 避免運算結果出錯，而引入 `clz` 便會在 n 值越小的時候，減少越多次迴圈運算，達成整體效能的優化，而非因為減少一次迴圈及取代減法。
- 改善 `bn_do_add` 的實作
透過將原先一次將 `c->number[0 ... c->size -1]` 全部運算的迴圈改成兩個迴圈，變成一個負責 `c->number[0 ... bsize -1]` ，若 `asize > bsize` 則再利用一個迴圈將 `c->number[bsize ... asize-1]` 算完。
不過我認為改善後的實作仍有改進空間，目前想到的有兩點：
1. $asize \geq bsize$
在改善後的實作中須確保 $asize \geq bsize$ 才可以後續運算，因此在 `bn_do_add` 中使用以下指令確保 $asize \geq bsize$ 恆成立。
```c=
if (a->size < b->size) // a->size >= b->size
SWAP(a, b);
```
但有沒有可能基於費式數列的特性，本來就可以確保 $asize \geq bsize$ 恆成立？
先來看看 `bn_add` 會使用在什麼地方，除了 `bn_sub` 中用來實現減法外，其餘都使用在計算費式數列中，讓我們列出 `bn_fib_fdoubling` 中跟 `bn_add` 有關的操作。
>1. bn_add(k1, f2, k1); // k1 = 2 * F(k-1) + F(k)
>2. bn_add(k2, k1, f1); // f1 = k1 + k2 = F(2k-1) now
>3. bn_add(f1, f2, f2); // f2 = F(2k+2)
接著將三個會用到 `bn_add` 的時機分開討論，先講 2 和 3 ，他們很直覺，在第 2 點中此時的 k1 和 k2 分別為 $F(k)^2$ 和 $F(k-1)^2$ ，由費式數列的單調遞增特性可知 $k1 \geq k2$ 恆成立，第三點中也可以很輕易的判斷執行 `bn_add` 時 $f1 \geq f2$ 恆成立。
而在第 1 點中， k1 和 f2 分別為 $2*F(k-1)$ 和 $F(k)$，接著可以經過一些簡單的推導得知，$k1 \geq f2$ 恆成立。
\begin{split}
2 * F(k-1) &\ \ \geq &\ \ F(k) \\
F(k-1) + F(k-1) &\ \ \geq &\ \ F(k-1)+F(k-2) \\
F(k-1) &\ \ \geq &\ \ F(k-2) \\
\end{split}
在確保所有使用 `bn_do_add` 的操作都可以確保 $a \geq b$ 恆成立後，便可以改寫 `bn_do_add`。
```diff=
/* |c| = |a| + |b|
+ * make sure asize always bigger than or equal to bsize
*/
static void bn_do_add(const bn *a, const bn *b, bn *c)
{
...
- if (a->size < b->size) // a->size >= b->size
- SWAP(a, b);
...
}
```
然後改寫 `bn_fib_fdoubling` 中的 `bn_add` (原本的寫法就已經符合 $a \geq b$ 因此僅加上註解)。
```diff=
void bn_fib_fdoubling(bn *dest, unsigned int n)
{
...
+ // k1 size always bigger than f2 size
bn_add(k1, f2, k1); // k1 = 2 * F(k-1) + F(k)
...
+ // k2 size always bigger than k1 size
bn_add(k2, k1, f1); // f1 = k1 + k2 = F(2k-1) now
if (n & i) {
...
+ // f1 size always bigger than f2 size
bn_add(f1, f2, f2); // f2 = F(2k+2)
}
}
...
}
```
2. 修改第二個迴圈
第二個迴圈主要的功能為將 a 中剩下的數值複製到 c 中並且加上 carry 。但考慮到進位的次數與 `asize - bsize` 的差值會隨著費式數列越大而越大，因此可以再拆解第二個迴圈，讓第二個迴圈負責處理 carry 而第三個負責處理剩下的 a->number 複製到 c->number 。
```diff=
/* |c| = |a| + |b| */
static void bn_do_add(const bn *a, const bn *b, bn *c)
{
...
bn_data carry;
carry = _add_partial(a->number, b->number, bsize, c->number);
if (asize != bsize) { // deal with the remaining part if asize > bsize
+ int i = bsize;
+ for (; i < asize & !!carry; i++) { // !!carry == 0 if carry == 0, == 1 if others
- for (int i = bsize; i < asize; i++) {
bn_data tmp1 = a->number[i];
carry = (tmp1 += carry) < carry;
c->number[i] = tmp1;
}
+ for (; i < asize; i++) {
+ c->number[i] = a->number[i];
+ }
}
if (carry) {
c->number[asize] = carry;
++(c->size);
}
}
```
#### 縮減記憶體操作成本
- 降低 malloc 與 memcpy 的使用
原先的實作(節錄片段)
```c=
void bn_fib_fdoubling(bn *dest, unsigned int n)
{
...
/* F(2k) = F(k) * [ 2 * F(k+1) – F(k) ] */
bn_cpy(k1, f2); // k1 = F(k+1)
bn_lshift(k1, 1); // k1 = 2 * F(k+1)
bn_sub(k1, f1, k1); // k1 = 2 * F(k+1) - F(k)
bn_mult(k1, f1, k1); // k1 = F(k) * [ 2 * F(k+1) - F(k) ]
...
bn_cpy(f1,k1); // f1 = k1
...
}
```
可以發現為了將運算結果都儲存在 `k1` 中，因此每次運算都會遇到以下問題：
1. 資料複製
2. 乘法運算時 `c == a or c == b` 動態配置暫存記憶體
而改進的做法便是重構算式，將 `k2` 也拿來儲存計算結果，從而避免問題 (2) ，再利用 `bn_swap` 省去 `memcpy` 的運用解決問題 (1) 。
```c=
void bn_fib_fdoubling(bn *dest, unsigned int n)
{
...
/* F(2k) = F(k) * [ 2 * F(k+1) – F(k) ] */
/* F(2k+1) = F(k)^2 + F(k+1)^2 */
bn_lshift(f2, 1, k1);// k1 = 2 * F(k+1)
bn_sub(k1, f1, k1); // k1 = 2 * F(k+1) – F(k)
bn_mult(k1, f1, k2); // k2 = k1 * f1 = F(2k)
bn_mult(f1, f1, k1); // k1 = F(k)^2
bn_swap(f1, k2); // f1 <-> k2, f1 = F(2k) now
...
}
```
- 引入 memory pool 的概念
在 `bn` 中加入 `capacity` 儲存實際配置的長度，而 `size` 變成當前使用的大小，進而減少頻繁呼叫 `realloc` 所產生的效能損耗。
另外在 [KYG-yaya573142 的報告](https://hackmd.io/@KYWeng/rkGdultSU)中提到
>至少能正確計算至第 100000 項
是因為 unsigned int 能表達的數值上界為 4294967295($2^{32}-1$，在 LP64 或 LLP64 系統)，因此能表達的費式數列上界為 $2^{32*(2^{32}-1)}$ ，考慮可能不需要計算到這麼後面的費式數列值，還可以引入 bit field 的技巧將 size,capacity 和 sign 整合。
```diff
typedef struct _bn {
unsigned int *number; /* ptr to number */
- unsigned int size; /* length of number */
- unsigned int capacity; /* total allocated length, size <= capacity */
- int sign;
+ unsigned int size: SIZE;
+ unsigned int capacity: CAPACITY;
+ unsigned int sign: SIGN;
} bn;
```
更進一步的節省 `bn` 的空間，不過這樣改寫的話原先很多函式都要一並改寫，整體利弊尚未可知。
- 因 word size 制宜
根據當前系統不同的 word size 調整 `bn` 中各項屬性的資料型態，不僅能使 `bn` 能夠在不同系統運行，同時也確保當使用 64-bit 的 CPU 時能夠增加計算效能與儲存空間。
:::warning
TODO:
針對 [KYG-yaya573142 的報告](https://hackmd.io/@KYWeng/rkGdultSU)減少大數運算的成本中的以下章節撰寫報告
- [ ] 改善 bn_mult 的效能
- [ ] 引入 bn_sqr
- [ ] 實作 Karatsuba algorithm
:::
---
## mutex 在 fibdrv 中的應用場景與移除時會發生的事情
參考了 [KYG-yaya573142 的報告](https://hackmd.io/@KYWeng/rkGdultSU)撰寫一個多執行緒程式來測試。
```c=
void *worker(void *arg){
int fd = open("/dev/fibonacci",O_RDWR);
if (fd < 0){
perror("Failed to open character device");
goto end;
}
int idx = *((int *)arg);
char buf[100];
for (int offset = 80; offset <=92; offset++){
lseek(fd,offset,SEEK_SET);
long long result = read(fd,buf,0);
printf("thread %d F(%d): %lld\n",idx,offset,result);
}
end:
close(fd);
}
int main(){
pthread_t test[2];
int idx[2] = {1,2};
pthread_create(&test[0],NULL,worker,&idx[0]);
pthread_create(&test[1],NULL,worker,&idx[1]);
for(int i = 0;i < 2; i++)
pthread_join(test[i],NULL);
return 0;
}
```
- 未移除 mutex
在原先的程式中，利用 `mutex_trylock` 和 `mutex_unlock` 來避免程式進入 critical section ，從而保證每個執行緒的內容正確。而其中使用的 `mutex_trylock` 當發現該區段已經被 lock 時，會直接返回，屬於非阻斷式 I/O 。
而這樣的特點也可以經由實驗發現，當反覆嘗試原本的 fibdrv 後，會出現兩種結果，一種是出現錯誤訊息的，另一種是兩個執行緒都成功的完成任務。
```
$ sudo ./thread
thread 1 F(80): 23416728348467685
Failed to open character device: Device or resource busy
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
thread 1 F(83): 99194853094755497
thread 1 F(84): 160500643816367088
```
```
$ sudo ./thread
thread 1 F(80): 23416728348467685
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
...
thread 2 F(80): 23416728348467685
thread 2 F(81): 37889062373143906
thread 2 F(82): 61305790721611591
...
```
第一種情況發生的原因是因為兩個執行緒針對 fibdrv 的讀寫順序是相互交雜的，因此當 fibdrv 已經被 lock 時，另一個執行緒嘗試 trylock 失敗，因此印出錯誤訊並結束該執行緒。
而第二種狀況則是運氣好，兩個執行緒的讀寫沒有交雜，因此兩個執行緒都完成各自的任務。
- 將 trylock 改成 lock
在改成使用 `mutex_lock` 後，當執行緒程式存取 mutex 時，若該區段已被 lock ，則會持續等待直到 unlock ，屬於阻斷式 I/O 。
```
$ sudo ./thread
thread 1 F(80): 23416728348467685
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
...
thread 2 F(80): 23416728348467685
thread 2 F(81): 37889062373143906
thread 2 F(82): 61305790721611591
...
```
- 去除 mutex
去除 mutex 後， thread 將可以依照原先的順序讀取，不用等待其他 thread 釋出程式驅動。
```
$ sudo ./thread
thread 2 F(80): 23416728348467685
thread 2 F(81): 37889062373143906
thread 2 F(82): 61305790721611591
thread 2 F(83): 99194853094755497
thread 1 F(80): 23416728348467685
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
thread 2 F(84): 160500643816367088
thread 2 F(85): 259695496911122585
thread 1 F(83): 99194853094755497
thread 1 F(84): 160500643816367088
```
可以看到 thread 1 和 2 交互讀取 fibdrv ，不過有趣的是可以發現其中每個費氏數列的值並沒有出錯，這是因為在原先撰寫的測試程式中， file descriptor 並不是共享資源， thread 間沒有共用同一個 fd 自然也不會有衝突。
在把 fd 改為全域變數之後，可以發現 thread 所得到的值出現錯誤，意味著進入了 critical section 。
```
$ sudo ./thread
thread 1 F(80): 23416728348467685
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
...
thread 1 F(89): 1779979416004714189
thread 2 F(80): 1779979416004714189
thread 2 F(81): 37889062373143906
thread 2 F(82): 61305790721611591
```
不過猜測由於測試程式僅有兩個執行緒，因此在進行 lseek 和 read 的操作時，要出現如 thread 1 lseek 但還沒 read 時， offset 已經被 tread 2 改變造成結果不正確的 race condition 反而機率不大，如果採用像 [ KYG-yaya573142 的報告](https://hackmd.io/@KYWeng/rkGdultSU)中用 sleep 的方法，就可以很明顯的看到效果。
```
$ sudo ./thread
thread 1 F(80): 23416728348467685
thread 2 F(80): 37889062373143906
thread 1 F(81): 37889062373143906
thread 1 F(82): 61305790721611591
thread 2 F(81): 99194853094755497
thread 1 F(83): 61305790721611591
thread 1 F(84): 160500643816367088
thread 2 F(82): 259695496911122585
```
或是單純的增加 thread 也可以較容易觀察到錯誤（增加到 20 個 thread）
```
sudo ./thread
...
thread 18 F(83): 99194853094755497
thread 6 F(84): 61305790721611591
thread 3 F(84): 160500643816367088
thread 19 F(80): 23416728348467685
...
```
---
## 額外議題
### xs 的設計與記憶體開銷
在實作大數加法的時候，我直接利用 char* 來儲存數字，企圖"簡化" xs 的架構來達到同樣可以大數運算的效果，但就在多看幾次 [SSO (Small/Short String Optimization)](http://wdv4758h.github.io/notes/cpp/string-optimization.html) 和聽老師上課後才稍微理解為何 xs 要這樣設計。
- xs 的架構
```c=
typedef union {
/* allow strings up to 15 bytes to stay on the stack
* use the last byte as a null terminator and to store flags
* much like fbstring:
* https://github.com/facebook/folly/blob/master/folly/docs/FBString.md
*/
char data[16];
struct {
uint8_t filler[15],
/* how many free bytes in this stack allocated string
* same idea as fbstring
*/
space_left : 4,
/* if it is on heap, set to 1 */
is_ptr : 1, is_large_string : 1, flag2 : 1, flag3 : 1;
};
/* heap allocated */
struct {
char *ptr;
/* supports strings up to 2^MAX_STR_LEN_BITS - 1 bytes */
/* MAX_STR_LEN_BITS is 54 in this case,
* so xs supports strings up to 2^54 -1 bytes */
size_t size : MAX_STR_LEN_BITS,
/* capacity is always a power of 2 (unsigned)-1 */
capacity : 6;
/* the last 4 bits are important flags */
};
} xs;
```
可以看到 xs 在佔用 16 bytes 且小字串的情況下，可以利用 15 bytes 儲存字串資料，利用 space_left 來計算剩餘可用的 bytes ，且 15 bytes 可用 4 個 bits 就可以儲存資料剩餘大小 ($2^4-1 = 15$) 。最後剩下 4 個 bits 的 `is_ptr` 和 `is_large_string` 可以分別對應到 [FBString](https://github.com/facebook/folly/blob/main/folly/docs/FBString.md) 中提到的中字串與大字串，而 `flag2` 與 `flag3` 我認為只是把剩下的 bits 定義完，實際上並不會用到。
至於在處理中字串與大字串時則利用 `heap allocated` 的 struct 來表達。
- xs 記憶體開銷
xs 這樣的設計使得在小字串時可以直接用 stack 中的記憶體，避免動態配置新記憶體，空間利用上也更完善，直接在原先配置的空間中紀錄了 size, capaity 和盡可能的儲存字串資料。
---
Sign in with Wallet
Connect another wallet