# 2020q3 Homework6 (kcalc) contributed by <`joe-U16`> [Github](https://github.com/joe-U16/kcalc-fixed) ###### tags: `sysprog2020` ## :penguin: 作業要求 * 在 GitHub 上 fork [kcalc-fixed](https://github.com/sysprog21/kcalc-fixed)，目標是擴充數值運算器的功能，例如新增 `sqrt`, `Sigma` ($\sum$), `Pi` ($\Pi$) 等等，可從 [math-expression-evaluator](https://github.com/bugwheels94/math-expression-evaluator) 的 "Supported symbols" 選出幾項來實作。過程中應一併完成以下: - [x] 擴充數值運算器的時候，應一併提供測試案例，且該儘量反映出有效的數值範圍 - [ ] 留意乘法和加法的 overflow、減法的 underflow 等處理機制，分析後予以強化 - [ ] 確保數值運算的 precision 和 accuracy，甚至要考慮 [Kahan summation algorithm](https://en.wikipedia.org/wiki/Kahan_summation_algorithm) $\to$ 你可以修改原有 fixed-point 的設計和實作方式 - [ ] 觀察原始程式碼的 `expression.[ch]` 裡頭 `vec` 相關的程式碼，提出對空間效率更好的實作。考慮到 `kcalc-fixed` 的使用情境往往是持續接受一系列字串輸入，應避免頻繁的記憶體配置和釋放 - [ ] 修正現有 `kcalc-fixed` 內部數值運算的錯誤，並留意到不合法的記憶體操作; - [ ] 指出原本 fixed-point 實作的限制，觀察 Linux 核心所用的 fixed-point 運算程式碼，提出針對 kcalc 的改進方案 - [x] 原有的測試程式不會檢驗輸出數值，應予以擴充，得以自動和預先輸入的期望數值進行比對; [TOC] ## 同學做的 - [ ] https://hackmd.io/@sysprog/2020-kcalc - [ ] https://hackmd.io/@lQBEPrYxRoOVV2MatUYusg/r1dvv-G38 - [ ] https://hackmd.io/@sammer1107/sysprog_hw6_kcalc - [x] https://hackmd.io/@YLowy/Sy9fENauv - [ ] https://hackmd.io/@RinHizakura/HyH6mtK_D - [x] https://hackmd.io/@Rsysz/kcalc - [ ] https://hackmd.io/@guaneec/sysprog-2020q3-kcalc ## 修正現有 kcalc-fixed 內部數值運算的錯誤，並留意到不合法的記憶體操作; Logic 的 `<`, `>`, `==`, `!=`, `>=`, `<=` 是有問題的 ``` Testing 1<2 ... 0.0000000037252903 Testing 2>1 ... 0.0000000037252903 Testing 1==1 ... 0.0000000037252903 Testing 1!=2 ... 0.0000000037252903 Testing 1>=1 ... 0.0000000037252903 Testing 1<=1 ... 0.0000000037252903 ``` #### `kcalc` 使用 32-bit 整數來表達單精度浮點數(即 `float`) **Fixed-point representation** kcalc introduces a fixed-point number representation for fractional values: one 32-bit size is divided into mantissa (28 bits), sign (1 bit) and exponent (3 bits). ``` MSB 31 4 3 2 1 0 LSB +----------------------------------+------+----------+ | mantissa | sign | exponent | +----------------------------------+------+----------+ ``` ```c= #define MASK(n) (((n) > 0) << 4) enum compare_type { LOWER = 1, EQUAL = 2, GREATER = 4 }; static int compare(uint64_t a, uint64_t b) { int flags = 0; if (a < b) flags |= LOWER; if (a > b) flags |= GREATER; if (a == b) flags |= EQUAL; return flags; } ``` #### Logic `!(not)` `expression.c` ```c=342 static uint64_t not(uint64_t a) { // return !a; return 1000; } ``` ``` 1.0000023283064370 Testing !0 ... 0.0000002328306437 Testing !1 ... ``` **改為 `return 0`：** ```c=342 static uint64_t not(uint64_t a) { // return !a; return 0; } ``` ``` Testing !0 ... 0.0000000000000000 Testing !1 ... 0.0000000000000000 ``` ## 擴充數值運算器的功能 ### 嘗試增加 user_func_add 首先註冊 user_func_add ```c= static struct expr_func user_funcs[] = { {"add", user_func_add, user_func_nop_cleanup, 0}, {NULL, NULL, NULL, 0}, }; ``` 實作 `user_func_add`: ```c= noinline uint64_t user_func_add(struct expr_func *f, vec_expr_t args, void *c) { uint64_t a = expr_eval(&args.buf[0]); uint64_t b = expr_eval(&args.buf[1]); return a + b; } ``` 測試結果： ``` Testing add(1,2) ... OK: 3.000000 == 3.000000 Testing add(3,4) ... OK: 7.000000 == 7.000000 Testing add(3.14, 2.11) ... OK: 5.250000 == 5.250000 ``` ### Sigma 實作 根據[`MathEX` 基本的使用方式：](https://hackmd.io/@sysprog/2020-kcalc#-數學式計算器) ```c= static struct expr_func user_funcs[] = { {"sigma", user_func_sigma, user_func_nop_cleanup, 0}, {"add", user_func_add, user_func_nop_cleanup, 0}, {NULL, NULL, NULL, 0}, }; ``` 實作 `user_func_sigma`: ```c= noinline uint64_t user_func_sigma(struct expr_func *f, vec_expr_t args, void *c) { int64_t lower = expr_eval(&args.buf[2]); int64_t upper = expr_eval(&args.buf[3]); /* return if bad function call */ if (args.len != 4) return NAN_INT; if (lower > upper) return NAN_INT; uint64_t sum = 0; int64_t *i = (int64_t *) args.buf[0].param.var.value; int max_iter = INT_MAX; for (*i = lower; *i <= upper && max_iter >= 0; *i += (1LL << 32), max_iter--) { uint64_t tmp = expr_eval(&args.buf[1]); if (tmp == NAN_INT || tmp == INF_INT) return tmp; uint64_t new_sum; if (__builtin_add_overflow(sum, tmp, &new_sum)) { pr_info("calc: sigma, overflow occur\n"); return INF_INT; } sum = new_sum; } if (max_iter <= 0) return NAN_INT; return sum; } ``` 使用方式: 需輸入 4 個 augument ，分別為 notation, expression, start, stopㄅ e.g. `sigma(n, n*3, 1, 10)` 等價於數學式 $\sum\limits_{n = 1}^{10}{n\times3}$ 測試結果： ``` Testing sigma(n, n*3, 1, 10) ... OK: 165.000000 == 165.000000 Testing sigma(n, n, 1, 10) ... OK: 55.000000 == 55.000000 Testing sigma(i, i**2, 1, 10) ... OK: 385.000000 == 385.000000 Complete ``` ### Sqrt 實作 使用 `@RinHizakura` 同學的[程式碼](https://hackmd.io/@RinHizakura/HyH6mtK_D#sqrt) ```c= noinline uint64_t user_func_sqrt(struct expr_func *f, vec_expr_t args, void *c) { int64_t ix0 = expr_eval(&vec_nth(&args, 0)); /* for 0 or NAN or INF, just return itself*/ if (ix0 == 0 || ix0 == NAN_INT || ix0 == INF_INT) return ix0; /* first, scale our number between 1 to 4 */ int lz = __builtin_clzll(ix0); /* for negative number */ if (lz == 0) return NAN_INT; /* for range from 0 to 30 */ if (lz <= 30) { if (lz & 1) lz++; ix0 >>= (30 - lz); } /* for range from 31 to 63 */ else { if (!(lz & 1)) lz++; ix0 <<= (lz - 31); } int64_t s0, q, t, r; /* generate sqrt(x) bit by bit */ q = s0 = 0; /* [q] = sqrt(x) */ r = 0x400000000; /* r = moving bit from right to left */ while (r != 0) { t = s0 + r; // t = s_i + 2^(-(i+1)) if (t <= ix0) { // t <= y_i ? s0 = t + r; // s_{i+1} = s_i + 2^(-i) ix0 -= t; // y_{i+1} = yi - t q += r; // q_{i+1} = q_{i}+ 2^(-i-1) } ix0 += ix0; r >>= 1; } if (lz < 31) return (q >> 1) << ((30 - lz) >> 1); return (q >> 1) >> ((lz - 31) >> 1); } ``` 測試結果： ``` Testing sqrt(10) ... OK: 3.162278 == 3.162278 Testing sqrt(20) ... OK: 4.472136 == 4.472136 Testing sqrt(30) ... OK: 5.477226 == 5.477226 Testing sqrt(40) ... OK: 6.324555 == 6.324555 Testing sqrt(50) ... OK: 7.071068 == 7.071068 ``` 超過範圍或是不合法的輸入，也可以判斷出來。 ``` Testing INT_MAX=2147483647, sigma(n, n, 1, INT_MAX) ... OK: INF_INT == INF_INT Testing sigma(m, n, 2, 1) ... OK: NAN_INT == NAN_INT ``` #### 有效的數字範圍 **整數：** `INT32` 可表示的上限 `INT_MAX=2147483647` ```bash= test_op 'INT_MAX=2147483647, sqrt(INT_MAX)' '46340.950001' test_op 'INT_MAX=2147483647, sqrt(INT_MAX+1)' '46340.950012' ``` 超過上限整數上限會有問題 ``` Testing INT_MAX=2147483647, sqrt(INT_MAX) ... OK: 46340.950001 == 46340.950001 Testing INT_MAX=2147483647, sqrt(INT_MAX+1) ... FAIL: 0.000000 != 46340.950012 ``` **精準度瑕疵** 觀察至小數點後 16 位 * scripts/test.sh: ```bash= test_op 'sqrt(2)' ``` ``` Testing sqrt(2) ... 1.4142135622478587 ``` * sqrt() function in C Library `math.h` ```c= double b = sqrt(2); ``` ``` >b 1.4142135623730951 ``` 觀察兩者： 1.414213562==2478587== # scripts/test.sh 1.414213562==3730951== # math.h 可以發現到精度不足的問題。 ## 新增功能，可以自動和預先輸入的期望數值進行比對; ```c= ➜ kcalc-fixed git:(master) scripts/test.sh [sudo] password for papd: filename: /home/papd/kcalc-fixed/calc.ko version: 0.1 description: Math expression evaluation author: National Cheng Kung University, Taiwan license: Dual MIT/GPL srcversion: 1F24815BA4159D060435AAD depends: retpoline: Y name: calc vermagic: 5.4.0-58-generic SMP mod_unload Testing x=5, x = x+1 ... FAIL: 6.000000 != 1.000000 Testing six=6, seven=7, six*seven ... OK: 42.000000 == 42.000000 Testing 小熊=6, 維尼=7, 小熊*維尼 ... OK: 42.000000 == 42.000000 Testing τ=1.618, 3*τ ... FAIL: 4.854000 != 4.854000 Testing $(τ, 1.618), 3*τ() ... FAIL: 4.854000 != 4.854000 Testing nop() ... �4h� Complete ``` ### 功能錯誤 ``` Testing τ=1.618, 3*τ ... FAIL: 4.854000 != 4.854000 Testing 111.2 ... FAIL: 111.200000 != 111.200000 ``` 明明 `printf` 出來一樣，卻判斷為 `FAIL`。 再測試看看，印到更後面： ```c= printf("k_result=%.16lf, ans = %.16f\n\n", k_result, ans); ``` 觀察到兩者實際是不一樣的 ``` inte=4, frac=-627065226, d_frac = -0.146000 result=3.8539999997834675, ans = 4.8540000000000001 inte=123, frac=528280976, d_frac = 0.123000 result=123.1229999997008093, ans = 123.1230000000000047 ``` 看了這篇[Back to Basic: 談浮點數的比較](https://novus.pixnet.net/blog/post/30654542)後，才知道使用 `==` 來判斷 `kcalc` 的計算結果和手動輸入的兩個浮點數，是很愚蠢的一件事。 可以改成讓兩束差在可接受的範圍內就好。 可以計算絕對誤差： ```c= if (fabs(a - b) <= 0.001) ``` 但是離 0 越遠，接鄰兩浮點數之間的差就越大。所有可正常表示的 float 當中，只有在 [-1, 1] 區間內，接鄰兩數差才會小於 FLT_EPSILON；若在 [1, 2] 區間，接鄰兩數差恰好為 FLT_EPSILON，負數亦然。隨著數字增加，FLT_EPSILON 很快就會小到不像話，到了 16.0f 時，接鄰兩數差將達到 FLT_EPSILON 的 16 倍；到了16777216.0f 時，接鄰兩浮點數差已經大於 1.0。 這就是為什麼在 [-1, 1] 區間外使用 DBL_EPSILON、FLT_EPSILON 並不合理，就如同死守公分做為誤差單位，對於精密機械而言大得誇張，但對於計算交通路程來說卻又小得可笑，還不如實際考慮需求而設定誤差門檻。 在數值範圍變化大的場合，更好的作法是隨著尺度調整誤差門檻，也就是運用相對誤差。其中一個簡單的想法概念碼如下： ```c= int AlmostEqual(float a, float b, float tolerance) { float absA = fabsf(a); float absB = fabsf(b); return fabsf(a - b) <= tolerance * fmaxf(absA, absB); } ``` ### 目前實作： ```c= if (fabs(result - test) > 0.01) ``` 在數值太大時會出錯： ``` Testing 500000000000.4 ... FAIL: 1783793664.400000 != 500000000000.400024 ```