2020q1 Homework4 (kcalc)

contributed by < KYG-yaya573142 >

H07: kcalc

預期目標

• 學習數值系統：浮點運算和定點數
• 整合 fibdrv 和 quiz4 開發成果
• 學習 Livepatch 運作機制
• 撰寫可適用於使用者和核心層級的程式
• 自動測試機制

Linux 核心中使用定點數的案例

我的電腦 linux kernel 版本是 5.3，以下資料都是參照此版本進行說明

$ uname -r
5.3.0-45-generic

定點數簡介

並不是所有硬體都有 FPU，而且 kernel 內也不建議進行浮點數計算 (可參閱 Linux 核心的浮點數運算)，因此須改用 定點數計算 (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 有一個用來計算定點數次方的函式 fixed_power_int，很適合用來當作理解定點數計算的範例

/**
 * 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) 可以知道 /proc/loadavg 記錄著 load average 相關的資訊

$ cat /proc/loadavg

觀察 /fs/proc/loadavg.c 會發現 /proc/loadavg 是一個虛擬檔案系統，還可以看到實際輸出資料的函式 loadavg_proc_show

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 了解 LOAD_INT 和 LOAD_FRAC 的實作

#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 可以看到負責更新avnrun 的函式

/*
 * 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 可以看到其函式實作

/*
 * 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 來計算 load average
• newload += FIXED_1-1 目的是讓所有的小數位數都進一位，但目前想不出來和 fixed_power_int 中的進位方式的差異是甚麼
• 可以看到 exp 是函式參數，也就是說會使用事先算好的對應數值，來避免進行複雜的浮點數計算

/*
 * 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 計算分析
• 註解有提到，如果會頻繁的計算 load average，須提高 fraction 的大小，否則會造成 rounding error 的累積
  • 可參閱你所不知道的 C 語言: 浮點數運算

MathEx

MathEx 如何解析給定字串

從 test-simple.c (也就是 README 中的示範程式碼) 可以得知 expr_create 是負責將字串解析為表達式 (expression) 的函式，再觀察可以發現 expr_create 中會使用 expr_next_token 尋找字串中 token 的位置，接著 expr_create 其餘的程式碼會對 token 進行分析

expr_next_token

• 負責回傳 token 的長度 (單位為 char)
• 會忽略 '#' 整行的內容
• 跳過空白與換行符號 (具體來說是跳過 isspace 為真的 char)
• 使用 flag 來建立防錯檢查機制，避免 invalid 的輸入，例如
  • EXPR_TOP 代表可接受運算子
  • EXPR_TOPEN 代表可接受左括號
  • EXPR_TCLOSE 代表可接受右括號

expr_create

token 的解析主要是在 for (;;) 這個大迴圈中進行，由於程式碼很長，請至 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

vec_expr_t es = vec_init();
vec_str_t os = vec_init();
vec_arg_t as = vec_init();

從 expression.h 查看相關的定義

#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 的最大容量

接著觀察相關的定義

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 目前看不懂用途

typedef vec(struct expr_string) vec_str_t;
struct expr_string {
    const char *s;
    int n;
};

• 主要儲存使用者自訂函數名稱的字串與長度、以及兩種括號 '(' 與 '{'

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 為例

static float add(struct expr_func *f, vec_expr_t args, void *c) {
...
}

接著定義一個 expr_func 資料結構，內有使用者自訂函數的名稱 "add" 與函數指標 add，並將此資料結構作為參數帶入 expr_create

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 的新函式

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 資料結構，連結原本的函式名稱及 patch 中定義的函式名稱

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，包含目標核心模組的名稱與對應的 klp_func 物件名稱

static struct klp_object objs[] = {
    {
        .name = "calc",
        .funcs = funcs,
    },
    {},
};

• 最後使用 klp_patch 資料結構來定義此 patch 的名稱與內容

static struct klp_patch patch = {
    .mod = THIS_MODULE,
    .objs = objs,
};

準備完 patch 的內容後，從 Livepatch - Loading 可以知道需要在 module_init 時呼叫 klp_enable_patch()

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 掛載

$ sudo insmod livepatch-calc.ko

用 dmesg 觀察是否掛載成功

$ 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 仍會導致編譯錯誤

ERROR: "expr_eval" [/home/kyg/Documents/kcalc/livepatch-calc.ko] undefined!

觀察 Makefile 會發現錯誤的原因是因為 livepatch-calc 的相依性少了 expression.o，參考 kbuild Makefile 將其加入

 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 的成果整合入 kcalc

kalc 原本已將 MathEX 整合，且在 main.c 中可以看到使用 MathEX API 的範例 user_func_nop 和 user_func_nop_cleanup，接下來仿造範例的用法，將 fibdrv 中計算費氏數列的函數以 MathEX 使用者自訂函數的形式進行改寫，並以 livepatch 的方式 patch 成正確的函式

首先在 kcalc 中先定義錯誤的使用者自訂函數 user_func_fib

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

static struct expr_func user_funcs[] = {
    {"nop", user_func_nop, user_func_nop_cleanup, 0},
+   {"fib", user_func_fib, 0},
    {NULL, NULL, NULL, 0},
};

簡單測試結果如預期

$ 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

/* 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 取用，並根據 kcalc 使用的定點數格式進行修改 (MSB 28 bits 為 mantissa)

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 中的定點數格式如下

MSB   31                               4    3   2    1   0  LSB
     +----------------------------------+------+----------+
     |     mantissa                     | sign | exponent |
     +----------------------------------+------+----------+

kcalc 中使用定點數的方式和定點數簡介中說明的方式不同，以類似 floating-point 的方式定義

\(\text{value} = \text{mantissa} \times 10^{exponent}\)

• 以 10 進位為基礎
• 不同於 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 中的 union 用法，嘗試應用在 kcalc

typedef union {
    int data;
    struct {
        int frac : 4,
        mantissa : 28;
    };
} fxp;

• 優點是 frac 的部分不需使用 GET_FRAC 就可直接取用，以及可以節省一些 bit shift 的操作
• 缺點是需特別注意 Endianness 對 bitfield 的影響，因為順序是 implementation-defined (C99 6.7.2.1-10 僅提到 bitfield 會相鄰，但沒規範順序)
  • 這同時會降低 kcalc 的 portability
  • 可參考 CMU 的範例

將 quiz4 使用到的技巧予以應用，提出對空間效率更好的實作 (待補)

從 expression.h 可以看到 MathEX 使用記憶體的方式和 quiz4 中的例子相似，最明顯的不同處為 quiz4 會先使用 stack，超過 capacity 才動態配置記憶體，而 MathEX 則是直接使用動態記憶體

/* 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 沒有做處理

除法

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
  • 但不適合用在 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 的判斷式，實作的程式碼如下

#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 中的 S32_MAX

加法

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

減法

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;
    }
...
}

• 跟加法的處理方式一樣

右移

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 的運算有誤，例如

$ echo -ne "10**2\0" > /dev/calc
$ dmesg
[1465823.655749] calc: Received 6 -> 10**2
[1465823.655763] calc: Result: 19
$ fromfix 19
1000

觀察後發現原本的寫法會讓次方的計數出錯，修改後的程式碼如下

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\)

$ 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 的值並調整整數的部分)，實作如下

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);
}

測試後，結果如下

$ 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

目前仍有多項缺點待改善

• 除了 power 之外，mult 等函式都有一樣的問題
• 這麼做會造成誤差在數值的末端開始累積，例如 \(3.14^4 = 97.211712\)，但是目前算出來是 \(97.211574\)
• frac 的部分如果超過定點數可用的範圍，最終會造成小數點位置的誤差，須對 FP2INT 以及 expr_parse_number 進行修正

reference

系統負荷是怎麼計算的?
Load average explained
Built-in Functions to Perform Arithmetic with Overflow Checking

tags: linux 2020