# 2020q1 Homework4 (kcalc)
contributed by < [`KYG-yaya573142`](https://github.com/KYG-yaya573142/kcalc) >
> [H07: kcalc](https://hackmd.io/@sysprog/linux2020-kcalc)
## 預期目標
* 學習數值系統：浮點運算和定點數
* 整合 fibdrv 和 quiz4 開發成果
* 學習 [Livepatch](https://www.kernel.org/doc/Documentation/livepatch/livepatch.txt) 運作機制
* 撰寫可適用於使用者和核心層級的程式
* 自動測試機制
## Linux 核心中使用定點數的案例
我的電腦 linux kernel 版本是 5.3，以下資料都是參照此版本進行說明
```shell
$ uname -r
5.3.0-45-generic
```
### 定點數簡介
並不是所有硬體都有 [FPU](https://en.wikipedia.org/wiki/Floating-point_unit)，而且 kernel 內也不建議進行浮點數計算 (可參閱 [Linux 核心的浮點數運算](https://hackmd.io/@sysprog/linux2020-kcalc#-Linux-%E6%A0%B8%E5%BF%83%E7%9A%84%E6%B5%AE%E9%BB%9E%E6%95%B8%E9%81%8B%E7%AE%97))，因此須改用 [定點數計算 (Fixed-point arithmetic)](https://en.wikipedia.org/wiki/Fixed-point_arithmetic)
* 小數點的位置固定
* 以 10 進位為例，精確度到小數點後 3 位，$number_{10} \times 10^3 = 定點數_{10}$
* 以 2 進位為例，精確度到小數點後 4 位，$number_{2} \times 2^4 = 定點數_{2}$
定點數的運算有以下法則
* 定點數加法與減法 - 可直接進行
* 定點數乘法 - 結果須再**除** $b^n$ (b 為進制，n 為小數位數)
* 定點數除法 - 結果須再**乘** $b^n$ (b 為進制，n 為小數位數)
### `fixed_power_int`
在研究 Load Average 的時候發現 [linux/kernel/sched/loadavg.c](https://elixir.bootlin.com/linux/v5.3/source/kernel/sched/loadavg.c#L110) 有一個用來計算定點數次方的函式 `fixed_power_int`，很適合用來當作理解定點數計算的範例
```cpp
/**
* fixed_power_int - compute: x^n, in O(log n) time
*
* @x: base of the power
* @frac_bits: fractional bits of @x
* @n: power to raise @x to.
*
* By exploiting the relation between the definition of the natural power
* function: x^n := x*x*...*x (x multiplied by itself for n times), and
* the binary encoding of numbers used by computers: n := \Sum n_i * 2^i,
* (where: n_i \elem {0, 1}, the binary vector representing n),
* we find: x^n := x^(\Sum n_i * 2^i) := \Prod x^(n_i * 2^i), which is
* of course trivially computable in O(log_2 n), the length of our binary
* vector.
*/
static unsigned long
fixed_power_int(unsigned long x, unsigned int frac_bits, unsigned int n)
{
unsigned long result = 1UL << frac_bits;
if (n) {
for (;;) {
if (n & 1) {
result *= x;
result += 1UL << (frac_bits - 1);
result >>= frac_bits;
}
n >>= 1;
if (!n)
break;
x *= x;
x += 1UL << (frac_bits - 1);
x >>= frac_bits;
}
}
return result;
}
```
* 以 2 進位來思考次方項 $x^n = x^{\sum n_i \cdot 2^i} ,\ n_i\in{0, 1}$
* 例如 $x^{11} = x^{1 \cdot 2^0} \times x^{1 \cdot 2^1} \times x^{0 \cdot 2^2} \times x^{1 \cdot 2^3}$
* `if (n & 1)` 是否為 true 即對應到 $n_i$ 是否為 1
* `1UL << frac_bits` 代表定點數 1 (2 進位，小數點後有 `frac_bits` 位)
* `x += 1UL << (frac_bits - 1)` 會讓 `x` 進行無條件進位
* 因為是進行定點數乘法，需再除 $2^{frac\_bits}$，也就是 `x >>= frac_bits`
### `loadavg_proc_show`
從 [man proc(5)](http://man7.org/linux/man-pages/man5/proc.5.html) 可以知道 `/proc/loadavg` 記錄著 load average 相關的資訊
```shell
$ cat /proc/loadavg
```
觀察 [/fs/proc/loadavg.c](https://elixir.bootlin.com/linux/v5.3/source/fs/proc/loadavg.c#L33) 會發現 `/proc/loadavg` 是一個虛擬檔案系統，還可以看到實際輸出資料的函式 `loadavg_proc_show`
```cpp
static int loadavg_proc_show(struct seq_file *m, void *v)
{
unsigned long avnrun[3];
get_avenrun(avnrun, FIXED_1/200, 0);
seq_printf(m, "%lu.%02lu %lu.%02lu %lu.%02lu %ld/%d %d\n",
LOAD_INT(avnrun[0]), LOAD_FRAC(avnrun[0]),
LOAD_INT(avnrun[1]), LOAD_FRAC(avnrun[1]),
LOAD_INT(avnrun[2]), LOAD_FRAC(avnrun[2]),
nr_running(), nr_threads,
idr_get_cursor(&task_active_pid_ns(current)->idr) - 1);
return 0;
}
```
觀察 `seq_printf` 使用的格式 `%lu.%02lu`，然後從 [include/linux/sched/loadavg.h](https://elixir.bootlin.com/linux/v5.3/source/include/linux/sched/loadavg.h#L43) 了解 `LOAD_INT` 和 `LOAD_FRAC` 的實作
```cpp
#define LOAD_INT(x) ((x) >> FSHIFT)
#define LOAD_FRAC(x) LOAD_INT(((x) & (FIXED_1-1)) * 100)
```
* `%lu.%02lu` 分別使用 `LOAD_INT` 與 `LOAD_FRAC` 來組成完整的 10 進位數值
* `LOAD_INT` 就是直接取用 interger
* `LOAD_FRAC` 會取用 fraction 的部分
* `FIXED_1-1` 可視為 fraction 的 mask
* 乘上 100 的目的是保留 10 進位下的小數點後 2 位
* `FIXED_1/200` 可理解為是 10 進位下的定點數 0.005，由於結果會保留小數點後 2 位，加上此數值目的是進行無條進件位
* load average 的資料實際上儲存在 `avnrun` 這個 global array
### Load Average
從 [linux/kernel/sched/loadavg.c](https://elixir.bootlin.com/linux/v5.3/source/kernel/sched/loadavg.c#L331) 可以看到負責更新`avnrun` 的函式
```cpp
/*
* calc_load - update the avenrun load estimates 10 ticks after the
* CPUs have updated calc_load_tasks.
*
* Called from the global timer code.
*/
void calc_global_load(unsigned long ticks)
{
unsigned long sample_window;
long active, delta;
sample_window = READ_ONCE(calc_load_update);
if (time_before(jiffies, sample_window + 10))
return;
/*
* Fold the 'old' NO_HZ-delta to include all NO_HZ CPUs.
*/
delta = calc_load_nohz_fold();
if (delta)
atomic_long_add(delta, &calc_load_tasks);
active = atomic_long_read(&calc_load_tasks);
active = active > 0 ? active * FIXED_1 : 0;
avenrun[0] = calc_load(avenrun[0], EXP_1, active);
avenrun[1] = calc_load(avenrun[1], EXP_5, active);
avenrun[2] = calc_load(avenrun[2], EXP_15, active);
WRITE_ONCE(calc_load_update, sample_window + LOAD_FREQ);
/*
* In case we went to NO_HZ for multiple LOAD_FREQ intervals
* catch up in bulk.
*/
calc_global_nohz();
}
```
`calc_load` 是實際計算 load average 的函式，在 [/include/linux/sched/loadavg.h](https://elixir.bootlin.com/linux/v5.3/source/include/linux/sched/loadavg.h#L18) 可以看到其函式實作
```cpp
/*
* a1 = a0 * e + a * (1 - e)
*/
static inline unsigned long
calc_load(unsigned long load, unsigned long exp, unsigned long active)
{
unsigned long newload;
newload = load * exp + active * (FIXED_1 - exp);
if (active >= load)
newload += FIXED_1-1;
return newload / FIXED_1;
}
```
* 根據 [exponential smoothing](https://en.wikipedia.org/wiki/Exponential_smoothing) 來計算 load average
* `newload += FIXED_1-1` 目的是讓所有的小數位數都進一位，但目前想不出來和 `fixed_power_int` 中的進位方式的差異是甚麼
* 可以看到 `exp` 是函式參數，也就是說會使用事先算好的對應數值，來避免進行複雜的浮點數計算
```cpp
/*
* These are the constant used to fake the fixed-point load-average
* counting. Some notes:
* - 11 bit fractions expand to 22 bits by the multiplies: this gives
* a load-average precision of 10 bits integer + 11 bits fractional
* - if you want to count load-averages more often, you need more
* precision, or rounding will get you. With 2-second counting freq,
* the EXP_n values would be 1981, 2034 and 2043 if still using only
* 11 bit fractions.
*/
extern unsigned long avenrun[]; /* Load averages */
extern void get_avenrun(unsigned long *loads, unsigned long offset, int shift);
#define FSHIFT 11 /* nr of bits of precision */
#define FIXED_1 (1<<FSHIFT) /* 1.0 as fixed-point */
#define LOAD_FREQ (5*HZ+1) /* 5 sec intervals */
#define EXP_1 1884 /* 1/exp(5sec/1min) as fixed-point */
#define EXP_5 2014 /* 1/exp(5sec/5min) */
#define EXP_15 2037 /* 1/exp(5sec/15min) */
```
* `EXP_1` 、`EXP_5` 和 `EXP_15` 是事先算好並轉換成定點數的數值
* 例如 $1/e^{(5/60)} \times \text{FIXED_1}\approx 1884_{10}$
* 使用 11 個 bits 作為 fraction，剩餘的部分是 interger
* 可參閱 [Linux Kernel Load Average 計算分析](http://brytonlee.github.io/blog/2014/05/07/linux-kernel-load-average-calc/)
* 註解有提到，如果會頻繁的計算 load average，須提高 fraction 的大小，否則會造成 rounding error 的累積
* 可參閱[你所不知道的 C 語言: 浮點數運算](https://hackmd.io/@sysprog/c-floating-point#Rounding)
## [MathEx](https://github.com/jserv/mathex)
### MathEx 如何解析給定字串
從 `test-simple.c` (也就是 README 中的示範程式碼) 可以得知 `expr_create` 是負責將字串解析為表達式 (expression) 的函式，再觀察可以發現 `expr_create` 中會使用 `expr_next_token` 尋找字串中 token 的位置，接著 `expr_create` 其餘的程式碼會對 token 進行分析
#### `expr_next_token`
* 負責回傳 token 的長度 (單位為 char)
* 會忽略 `'#'` 整行的內容
* 跳過空白與換行符號 (具體來說是跳過 [isspace](http://man7.org/linux/man-pages/man3/isspace.3p.html) 為真的 char)
* 使用 flag 來建立防錯檢查機制，避免 invalid 的輸入，例如
* `EXPR_TOP` 代表可接受運算子
* `EXPR_TOPEN` 代表可接受左括號
* `EXPR_TCLOSE` 代表可接受右括號
#### `expr_create`
token 的解析主要是在 `for (;;)` 這個大迴圈中進行，由於程式碼很長，請至 [expression.c](https://github.com/jserv/MathEX/blob/master/expression.c) 閱覽對應的原始碼。另外，`expr_create` 會使用三個 vector 物件 `es` 、 `os` 與 `as`，會在後面進行討論。
* 一開始會先略過空白與換行符號
* 684 行會將常數 token 直接存進 `es`
* 687 行會處理當 token 是 operator 的情況
* 如果是一般的 operator 就存進 `os`
* 如果是 comma operator 同時最接近 `os` 內容是 `'{'`，代表 operand 是函式的參數，會存進 `as`
* 719 行會處理當 token 是變數或函式名稱的情況
* 可以得知 `id` 與 `idn` 是變數或函式的名稱與名稱的長度
* 因此 543 行的 `if (idn > 0)` 是在處理 token 是變數與函式的情況
* 變數會存進 `es`
* 使用 `expr_func` 來判斷函式名稱是否有對應的定義，若有則將函式名稱存進 `os`
* 實際將函式指標登入至 `es` 的程式碼在 667 行
* 若 token 是函式名稱，會將 parent flag 改為 `EXPR_PAREN_EXPECTED` 來告知接下來遇到的左括號會是函式的左括號而不是一般的括號，並在 574 行改以 `'{'` 存入 `os` 作為標示
* 587 行會開始根據 `es` 、 `os` 與 `as` 的內容，將最終結果整合進 `es` (概念類似展開方程式)
* 最後在 750 行將 `es` 作為最後結果回傳
#### `es` 、 `os` 與 `as`
在 `expr_create` 開頭定義了三個 vector 物件 `es` 、 `os` 與 `as`
```cpp
vec_expr_t es = vec_init();
vec_str_t os = vec_init();
vec_arg_t as = vec_init();
```
從 [expression.h](https://github.com/jserv/MathEX/blob/master/expression.h) 查看相關的定義
```cpp
#define vec_init() { NULL, 0, 0 }
#define vec(T) \
struct { \
T *buf; \
int len; \
int cap; \
}
typedef vec(struct expr) vec_expr_t;
typedef vec(struct expr_string) vec_str_t;
typedef vec(struct expr_arg) vec_arg_t;
```
* 定義的 vector 儲存以下資訊
* `buf` 儲存的目標資料的起始位置
* `len` 為此 vector 現有的資料長度
* `cap` 為此 vector 的最大容量
接著觀察相關的定義
```cpp
typedef vec(struct expr) vec_expr_t;
struct expr {
int type;
union {
struct { float value; } num;
struct { float *value; } var;
struct { vec_expr_t args; } op;
struct {
struct expr_func *f;
vec_expr_t args;
void *context;
} func;
} param;
};
```
* `es` 可視為保存解析後的表達式的 vector
* `struct expr` 內還會有 `vec_expr_t`，代表可以是樹狀結構
* 根據 `type` 來決定用哪種方式取用 `union` 內的資料 (`type` 定義在 expression.c 的 `enum expr_type` 內)
* `num` 單純儲存數值
* `var` 儲存變數
* `op` 儲存 operands，使用的 operator 則是根據 `type` 決定
* `func` 儲存使用者自訂的函式以及對應的 operands
* `func.f` 會指向存放有對應函式的 `struct expr_func` 物件，該物件內的 `exprfn_t f` 才是真正的 function pointer
* `func.args` 存放使用者函數對應的 operands
* `func.context` 目前看不懂用途
```cpp
typedef vec(struct expr_string) vec_str_t;
struct expr_string {
const char *s;
int n;
};
```
* 主要儲存使用者自訂函數名稱的字串與長度、以及兩種括號 `'('` 與 `'{'`
```cpp
typedef vec(struct expr_arg) vec_arg_t;
struct expr_arg {
int oslen;
int eslen;
vec_expr_t args;
};
```
* 主要儲存函式的參數
* `oslen` 與 `eslen` 是用來定位相對於參照的函式或 operator 的位置
### `MathEx` 中，`nop` 一類使用者自訂的函數如何實作
為了使用自訂的函數，首先需先定義自訂函數，以自訂函數 `add` 為例
```cpp
static float add(struct expr_func *f, vec_expr_t args, void *c) {
...
}
```
接著定義一個 `expr_func` 資料結構，內有使用者自訂函數的名稱 `"add"` 與函數指標 `add`，並將此資料結構作為參數帶入 `expr_create`
```cpp
static struct expr_func user_funcs[] = {
{ "add", add, 0 },
{ NULL, NULL, 0 },
};
struct expr *e = expr_create(s, strlen(s), &vars, user_funcs);
```
在 `expr_create` 中會分析字串，跟自訂函數有關的部分主要是
* 719 行會判斷 token 可能是函式名稱
* 下一次迴圈，在 556 行使用 `expr_func` 來判斷函式名稱是否有對應的定義 (從 `user_funcs`)，若有則將函式名稱存進 `os`
* 最後在 667 行處會將函式指標登入至 `es`，至此就可以在表達式中呼叫使用者自訂的函數了
## livepatch
### 用法簡介
livepatch 模組的整體寫法與 kernel module 幾乎相同，kalc 中其實已經寫好可以運作的範本，以下簡單分析 kalc 使用 livepatch 的方式
* 首先定義欲 patch 的新函式
```cpp
void livepatch_nop_cleanup(struct expr_func *f, void *c)
{
...
}
int livepatch_nop(struct expr_func *f, vec_expr_t args, void *c)
{
...
}
```
* 建立 [`klp_func`](https://elixir.bootlin.com/linux/latest/source/include/linux/livepatch.h#L27) 資料結構，連結原本的函式名稱及 patch 中定義的函式名稱
```cpp
static struct klp_func funcs[] = {
{
.old_name = "user_func_nop",
.new_func = livepatch_nop,
},
{
.old_name = "user_func_nop_cleanup",
.new_func = livepatch_nop_cleanup,
},
{},
};
```
* 建立 patch 物件 [`klp_object`](https://elixir.bootlin.com/linux/latest/source/include/linux/livepatch.h#L106)，包含目標核心模組的名稱與對應的 `klp_func` 物件名稱
```cpp
static struct klp_object objs[] = {
{
.name = "calc",
.funcs = funcs,
},
{},
};
```
* 最後使用 [`klp_patch`](https://elixir.bootlin.com/linux/latest/source/include/linux/livepatch.h#L106) 資料結構來定義此 patch 的名稱與內容
```cpp
static struct klp_patch patch = {
.mod = THIS_MODULE,
.objs = objs,
};
```
準備完 patch 的內容後，從 [Livepatch - Loading](https://www.kernel.org/doc/html/latest/livepatch/livepatch.html#loading) 可以知道需要在 `module_init` 時呼叫 `klp_enable_patch()`
```cpp
static int livepatch_calc_init(void)
{
return klp_enable_patch(&patch);
}
static void livepatch_calc_exit(void)
{
klp_unregister_patch(&patch);
}
module_init(livepatch_calc_init);
module_exit(livepatch_calc_exit);
```
使用 kbuild Makefile 編譯後，就可將此 patch 掛載
```shell
$ sudo insmod livepatch-calc.ko
```
用 `dmesg` 觀察是否掛載成功
```shell
$ dmesg
[1144023.202333] livepatch: enabling patch 'livepatch_calc'
[1144023.216455] livepatch: 'livepatch_calc': starting patching transition
[1144024.750831] livepatch: 'livepatch_calc': patching complete
```
### 在 livepatch-calc 中使用 MathEX API
`livepatch-calc` 的目的是在不卸載核心模組的前提下，修補 kcalc 中錯誤的函式，因此在 `livepatch-calc` 中勢必會使用到 MathEX API。
然而實際使用後，會發現儘管已 `#include "expression.h"`，但使用 MathEX API 仍會導致編譯錯誤
```shell
ERROR: "expr_eval" [/home/kyg/Documents/kcalc/livepatch-calc.ko] undefined!
```
觀察 Makefile 會發現錯誤的原因是因為 `livepatch-calc` 的相依性少了 `expression.o`，參考 [kbuild Makefile](https://www.kernel.org/doc/Documentation/kbuild/makefiles.txt) 將其加入
```diff
KDIR=/lib/modules/$(shell uname -r)/build
obj-m += calc.o
obj-m += livepatch-calc.o
calc-objs += main.o expression.o
+livepatch-calc-y = livepatch-calc.o expression.o
ccflags-y := -std=gnu99 -Wno-declaration-after-statement
```
再次編譯會發現仍然出錯，使用 `make V=1` 查看錯誤訊息
```
...
make[2]: Circular /home/kyg/Documents/kcalc/livepatch-calc.o <- /home/kyg/Documents/kcalc/livepatch-calc.o dependency dropped.
...
ld -m elf_x86_64 -z max-page-size=0x200000 -r -o /home/kyg/Documents/kcalc/calc.o /home/kyg/Documents/kcalc/main.o /home/kyg/Documents/kcalc/expression.o
ld -m elf_x86_64 -z max-page-size=0x200000 -r -o /home/kyg/Documents/kcalc/livepatch-calc.o /home/kyg/Documents/kcalc/expression.o
```
發現如果 module 的名稱和 source file 相同，會導致 Makefile 誤判相依性，因此解決方法就是 module 和 source file 不要使用相同的名稱。由於 `test.sh` 等測試腳本都是使用 `livepatch-calc` 這個名稱，所以我將 `livepatch.c` 更改為 `klp-calc.c`，更改完後即可正常編譯
### 配合使用 livepatch 將 [fibdrv](https://hackmd.io/BZjwYw1FQ1O0NPY5kIS7VQ) 的成果整合入 kcalc
kalc 原本已將 [MathEX](https://github.com/jserv/mathex) 整合，且在 [main.c](https://github.com/KYG-yaya573142/kcalc/blob/master/main.c#L103) 中可以看到使用 MathEX API 的範例 `user_func_nop` 和 `user_func_nop_cleanup`，接下來仿造範例的用法，將 fibdrv 中計算費氏數列的函數以 MathEX 使用者自訂函數的形式進行改寫，並以 livepatch 的方式 patch 成正確的函式
首先在 kcalc 中先定義錯誤的使用者自訂函數 `user_func_fib`
```cpp
noinline int user_func_fib(struct expr_func *f, vec_expr_t args, void *c)
{
(void) c;
/* echo what you inputed */
return expr_eval(&vec_nth(&args, 0));
}
```
記得將此自訂函數使用 MathEX API 註冊，名稱使用 `fib`
```diff
static struct expr_func user_funcs[] = {
{"nop", user_func_nop, user_func_nop_cleanup, 0},
+ {"fib", user_func_fib, 0},
{NULL, NULL, NULL, 0},
};
```
簡單測試結果如預期
```shell
$ echo -ne "fib(666)\0" > /dev/calc
$ dmesg
[1155511.080853] calc: Device successfully closed
[1155542.870672] calc: Received 9 -> fib(666)
[1155542.870682] calc: Result: 10656
[1155542.870687] calc: Device successfully closed
$ fromfixed 10656
666
```
接下來要製作 patch，我選擇直接更新本來的 livepatch-calc.c
```diff
/* clang-format off */
static struct klp_func funcs[] = {
{
.old_name = "user_func_nop",
.new_func = livepatch_nop,
},
{
.old_name = "user_func_nop_cleanup",
.new_func = livepatch_nop_cleanup,
},
+ {
+ .old_name = "user_func_fib",
+ .new_func = livepatch_fib,
+ },
{},
};
```
`fib` 實作的部分直接從 [fibdrv](https://github.com/KYG-yaya573142/fibdrv/blob/master/fibdrv.c#L43) 取用，並根據 [kcalc 使用的定點數格式](#kcalc-使用的定點數格式)進行修改 (MSB 28 bits 為 mantissa)
```cpp
noinline int livepatch_fib(struct expr_func *f, vec_expr_t args, void *c)
{
if (args.len != 1) {
pr_err("fib(): too many or no argument\n");
return -1;
}
int n = expr_eval(&vec_nth(&args, 0)); /* fixed-point */
int frac = GET_FRAC(n); /* frac */
n >>= 4; /* mantissa */
/* fixed-point -> interger */
if (frac > 0) {
for (; frac; frac--)
n *= 10;
}
if (frac < 0) {
for (; frac; frac++)
n /= 10;
}
if (n < 2) { /* F(0) = 0, F(1) = 1 */
return (n << 4);
}
long long fib[2];
unsigned int ndigit = 32 - __builtin_clz(n);
fib[0] = 0; /* F(k) */
fib[1] = 1; /* F(k+1) */
/* fast doubling algorithm */
for (unsigned int i = 1U << (ndigit - 1); i; i >>= 1) {
long long k1 = fib[0] * (fib[1] * 2 - fib[0]);
long long k2 = fib[0] * fib[0] + fib[1] * fib[1];
if (n & i) {
fib[0] = k2;
fib[1] = k1 + k2;
} else {
fib[0] = k1;
fib[1] = k2;
}
}
return FP2INT(fib[0], 0);
}
```
* 因為 fib 是 MathEX 的使用者自訂函數，輸入與輸出的定點數格式需符合原本的格式
* 由於直接進行加法會相當複雜 (類似手動處理浮點數運算)，我選擇先將定點數轉為整數再進行運算
* `FP2INT` 和 `GET_FRAC` 是 expression.c 中的輔助函數
* kcalc 中定點數的轉換方式後文會進行討論
## kcalc
### kcalc 使用的定點數格式
kcalc 中的定點數格式如下
```shell
MSB 31 4 3 2 1 0 LSB
+----------------------------------+------+----------+
| mantissa | sign | exponent |
+----------------------------------+------+----------+
```
kcalc 中使用定點數的方式和[定點數簡介](#定點數簡介)中說明的方式不同，以類似 floating-point 的方式定義
$\text{value} = \text{mantissa} \times 10^{exponent}$
* 以 10 進位為基礎
* 不同於 [IEEE 754](https://en.wikipedia.org/wiki/IEEE_754) 中使用 bias 來調整 exponent 的正負號，kcalc 是將 LSB 第 4 bit 作為 sign bit 使用，因此可以將 exponent 視作長度為 4 bits 的 signed interger 進行思考
* 由於 `NAN_INT` 定義為 `(1 << 4) | ((1 << 4) - 1)`，所以當 exponent 是 -1 時，會將 mantissa 進位並調整 exponent 為 -2 (在 expression.c 中的函式 `FP2INT` 可以看到這種操作)
### 使用 union
參考 [quiz2](https://hackmd.io/OXUSiB-IQiqC17pIK367Yw) 中的 union 用法，嘗試應用在 kcalc
```cpp
typedef union {
int data;
struct {
int frac : 4,
mantissa : 28;
};
} fxp;
```
* 優點是 frac 的部分不需使用 `GET_FRAC` 就可直接取用，以及可以節省一些 bit shift 的操作
* 缺點是需特別注意 [Endianness](https://en.wikipedia.org/wiki/Endianness) 對 bitfield 的影響，因為順序是 implementation-defined (C99 6.7.2.1-10 僅提到 bitfield 會相鄰，但沒規範順序)
* 這同時會降低 kcalc 的 portability
* 可參考 [CMU 的範例](https://wiki.sei.cmu.edu/confluence/display/c/EXP11-C.+Do+not+make+assumptions+regarding+the+layout+of+structures+with+bit-fields)
### 將 quiz4 使用到的技巧予以應用，提出對空間效率更好的實作 (待補)
從 expression.h 可以看到 MathEX 使用記憶體的方式和 quiz4 中的例子相似，最明顯的不同處為 quiz4 會先使用 stack，超過 capacity 才動態配置記憶體，而 MathEX 則是直接使用動態記憶體
```cpp
/* Simple expandable vector implementation */
static inline int vec_expand(char **buf, int *length, int *cap, int memsz)
{
if (*length + 1 > *cap) {
void *ptr;
int n = (*cap == 0) ? 1 : *cap << 1;
#ifdef __KERNEL__
ptr = krealloc(*buf, n * memsz, GFP_KERNEL);
#else
ptr = realloc(*buf, n * memsz);
#endif
if (!ptr)
return -1; /* allocation failed */
*buf = (char *) ptr;
*cap = n;
}
return 0;
}
```
之前 quiz4 記憶體相關的分析尚未完成，待完成後會互相整合
### MathEx 對於負數的除法結果不正確、乘法和加法的 overflow 沒有做處理
#### 除法
```diff
static int divid(int a, int b)
{
int frac1 = GET_FRAC(a);
int frac2 = GET_FRAC(b);
int n1 = GET_NUM(a);
int n2 = GET_NUM(b);
if (n1 == 0 && n2 == 0)
return NAN_INT;
if (n2 == 0)
return INF_INT;
+ if (n1 < 0) {
+ n1 = -n1;
+ n2 = -n2;
+ }
while (n1 * 10 < ((1 << 27) - 1)) {
--frac1;
n1 *= 10;
printf("n1: %d, frac: %d\n", n1, frac1);
}
int n3 = n1 / n2;
int frac3 = frac1 - frac2;
return FP2INT(n3, frac3);
}
```
* 實際進行除法的 `int n3 = n1 / n2` 本來就能處理負數，導致問題的是 `while` 迴圈的判斷沒有考慮到 `n1 < 0` 的情況
* 因此只需要修正當 `n1 < 0` 的情況即可
#### 乘法
乘法時可能有 2 種產生 overflow 的情況
* 同號，可能產生 positive overflow
* 異號，可能產生 negative overflow
我能想到偵測 overflow 的方式有
* 將結果除以原本的 `n1`，若無 overflow 應該要等於 `n2`
* 但可能因數值的精確度限制而產生誤判
* 先產生 overflow 才偵測的到
* 根據與 `int` 最大 ($T{max}$) 與最小值 ($T{min}$) 的距離判斷
* 使用 gcc 的 [Built-in Functions to Perform Arithmetic with Overflow Checking](https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins.html)
* 但不適合用在 kcalc 中，因為是使用自訂的數值型態
考量 overflow 是 undefined behavior，最好的方式應該是預防 overflow 而不是發生後才偵測，因此最後決定使用第二個方法。首先根據 `n1` 與 `n2` 是否為同號來列出不產生 overflow 的條件
* 同號時
* $n_1 \times n_2 < T_{max}$
* 異號時
* $T_{min} < n_1 \times n_2$
根據上兩式，將 `n1` 或是 `n2` 移項就可以得到避免 overflow 的判斷式，實作的程式碼如下
```cpp
#define S28_MIN (-(1 << 27))
#define S28_MAX ((1 << 27) - 1)
#define FRAC_MIN (-(1 << 3))
#define FRAC_MIN ((1 << 3) - 1)
/* (2^27 - 1) ~ -2^27 */
static int mult(int a, int b)
{
if (!a || !b) {
return 0;
}
int frac1 = GET_FRAC(a);
int frac2 = GET_FRAC(b);
int n1 = GET_NUM(a);
int n2 = GET_NUM(b);
if ((n1 > 0 && (n2 > S28_MAX / n1 || n2 < S28_MIN / n1)) ||
(n2 == -1 && n1 == S28_MIN) ||
(n1 < 0 && (n1 < S28_MAX / n2 || n1 < S28_MIN / n2))) {
/* num overflow */
pr_err("mult: num overflow\n");
return -1;
}
if ((frac1 > 0 && frac1 > (FRAC_MAX - frac2)) ||
(frac1 < 0 && frac1 < (FRAC_MIN - frac2))) {
/* frac overflow */
pr_err("mult: frac overflow\n");
return -1;
}
int n3 = n1 * n2;
return FP2INT(n3, (frac1 + frac2));
}
```
* 特別注意 $-T{max}$ 是存在的數值，但是 $-T{min}$ 本身就會產生 overflow，因此需特別處理分母可能為 -1 的狀況
* `(n2 == -1 && n1 == S28_MIN)`
* 若 `n1` 或 `n2` 其中一者為 0 會直接回傳 0，所以不用考慮分母為 0 的情況
* `S28_MAX` 的名稱是參考 [linux/include/linux/limits.h](https://github.com/torvalds/linux/blob/master/include/linux/limits.h) 中的 `S32_MAX`
#### 加法
```cpp
static int plus(int a, int b)
{
...
if ((n1 > 0 && n2 > (S28_MAX - n1)) ||
(n1 < 0 && n2 < (S28_MIN - n1))) {
/* plus overflow */
pr_err("plus: overflow\n");
return -1;
}
...
}
```
* 跟乘法中 frac 的處理方式一樣
* 注意移項的問題，例如不能寫成 `n1 > 0 && n1 > (S28_MAX - n2)`，因為當 `n2` 為負數時 `S28_MAX - n2` 會發生 overflow
#### 減法
```cpp
static int minus(int a, int b)
{
...
if ((n1 < 0 && n2 > (S28_MAX + n1)) ||
(n1 > 0 && n2 < (S28_MIN + n1))) {
/* plus overflow */
pr_err("minus: overflow\n");
return -1;
}
...
}
```
* 跟加法的處理方式一樣
#### 右移
```diff
static int right_shift(int a, int b)
{
/* FIXME: should use 2-base? */
- return divid(a, mult(2 << 4, b));
+ return divid(a, power(2 << 4, b));
}
```
* 應該只是單純的寫錯成 `mult`
#### 次方
測試時發現 power 的運算有誤，例如
```shell
$ echo -ne "10**2\0" > /dev/calc
$ dmesg
[1465823.655749] calc: Received 6 -> 10**2
[1465823.655763] calc: Result: 19
$ fromfix 19
1000
```
觀察後發現原本的寫法會讓次方的計數出錯，修改後的程式碼如下
```diff
static int power(int a, int b)
{
int frac1 = GET_FRAC(a);
int frac2 = GET_FRAC(b);
int n1 = GET_NUM(a);
int n2 = GET_NUM(b);
while (frac2 != 0) {
if (frac2 < 0) {
n2 /= 10;
++frac2;
} else {
n2 *= 10;
--frac2;
}
}
int on1 = n1;
int of1 = frac1;
- n1 = 1;
- for (int i = 0; i < n2; ++i) {
+ for (int i = 1; i < n2; ++i) {
frac1 += of1;
n1 *= on1;
}
return FP2INT(n1, frac1);
}
```
#### 動態調整數值
因為 kcalc 使用定點數方式的關係，可能會產生計算結果是 kcalc 可容納的數值，但是計算過程會出現 kcalc 定點數格式無法容納的數值。
舉例來說，考量計算 $3.14^{4} \approx 97.211$
```shell
$ echo -ne "3.14**4\0" > /dev/calc
$ dmesg
[ 895.091423] calc: Received 8 -> 3.14**4
[ 895.091435] calc: Result: 919916808
$ fromfix 919916808
.57494800000000000000
```
這是因為 power 在計算時，實際上會計算 $314^4$，導致 overflow，因此我想到的方法為當超過可用範圍時，動態的移動數值的小數點位置 (也就是調整 exponent 的值並調整整數的部分)，實作如下
```diff
static int power(int a, int b)
{
int frac1 = GET_FRAC(a);
int frac2 = GET_FRAC(b);
int n1 = GET_NUM(a);
int n2 = GET_NUM(b);
while (frac2 != 0) {
if (frac2 < 0) {
n2 /= 10;
++frac2;
} else {
n2 *= 10;
--frac2;
}
}
int on1 = n1;
int of1 = frac1;
for (int i = 1; i < n2; ++i) {
+ while ((on1 > 0 && (n1 > S28_MAX / on1)) ||
+ (n1 == -1 && on1 == S28_MIN) ||
+ (on1 < 0 && (on1 < S28_MAX / n1 || on1 < S28_MIN / n1))) {
+ ++frac1;
+ n1 /= 10;
+ }
frac1 += of1;
n1 *= on1;
}
return FP2INT(n1, frac1);
}
```
測試後，結果如下
```shell
$ echo -ne "3.14**4\0" > /dev/calc
$ dmesg
[ 1256.659897] calc: Received 8 -> 3.14**4
[ 1256.659908] calc: Result: 1555385194
$ fromfixed 1555385194
97.21157400000000000000
```
:::warning
目前仍有多項缺點待改善
* 除了 `power` 之外，`mult` 等函式都有一樣的問題
* 這麼做會造成誤差在數值的末端開始累積，例如 $3.14^4 = 97.211712$，但是目前算出來是 $97.211574$
* frac 的部分如果超過定點數可用的範圍，最終會造成小數點位置的誤差，須對 `FP2INT` 以及 `expr_parse_number` 進行修正
:::
## reference
[系統負荷是怎麼計算的?](http://tinylab.org/how-to-calc-load-and-part1/)
[Load average explained](https://wiki.nix-pro.com/view/Load_average_explained)
[Built-in Functions to Perform Arithmetic with Overflow Checking](https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins.html)
###### tags: `linux 2020`
Sign in with Wallet
Connect another wallet