2020q3 Homework6 (kcalc)

# 2020q3 Homework6 (kcalc) contributed by < `sammer1107` > ###### tags: `進階電腦系統理論與實作` > [作業說明](https://hackmd.io/@sysprog/2020-kcalc) # 完善測試機制原本的 test case 跑完也不確定結果到底正不正確，所以我決定先從改善 test 機制開始。 ```bash= NAN="NAN_INT" INF="INF_INT" test_op() { local expression=$1 local ans=$2 echo -n "Testing " ${expression} "... " echo -ne ${expression}'\0' > $CALC_DEV local ret=$($EVAL $(cat $CALC_DEV)) if [[ "$ans" != $NAN && "$ans" != $INF ]]; then ans=$(printf %lf\\n $ans) fi if [ "$ret" = "$ans" ]; then echo -e "\e[32mPASSED\e[0m" echo -e "==" $ans else echo -e "\e[31mFAILED\e[0m" echo -e "Got" $ret "instead of" $ans fi echo } ``` + test_op 除了接受第一個參數作為測試輸入，現在也接收第二個參數，作為預期結果。 + 如果逾期結果是一個輸字，則會利用 printf 命令確保與 eval 的輸出一致。 + 成功的話會出現綠色的 PASSED ，失敗的話則會出現紅色的 FAILED 使用範例： ```bash # multiply test_op '6*7' 42 ``` 輸出: ``` Testing 6*7 ... PASSED == 42.000000 ``` ## 失敗的 case 根據原先註解預期的輸出，我發現有幾個 test case 目前會失敗： 1. ``` Testing (3%0)|0 ... FAILED Got NAN_INT instead of 0.000000 ``` 2. ``` Testing x=5, x=(x!=0) ... FAILED Got 0.000000 instead of 1.000000 ``` ### 對 0 取餘數 + 位元 OR + 照理來講，與 0 做位元 OR ，應該不會改變值才對，可是 test case 卻預期會把 NAN_INT 變為 0。這是個奇怪的行為。 + 追查到原本的 MathEx 專案，也寫了一個一樣的 test case。而之所以會有這樣的答案，是因為 `|` 的操作會先將數字轉為 int 才做位元運算，而原先的轉換規則會將 NaN 轉為 0，才會造成 `0|0 = 0`。 + 如果想要維持一樣的行為，則需將 kcalc 的程式碼做更動： ```c case OP_BITWISE_OR: /* OK */ return GET_NUM(expr_eval(&e->param.op.args.buf[0])) | GET_NUM(expr_eval(&e->param.op.args.buf[1])); ``` 用 GET_NUM 來將小數轉為整數，如此 + 正常的數字會被截短，符合原本的行為 + nan 的整數部份為 0 ，會轉為 0。符合原本的行為 + inf 目前的設計是最高位一個 1 ，其他為 0 ，如此會維持原本的數值。 + 但原先 MathEx 的設計是轉成最大的整數 ==正負 INT_MAX== (雖然最小應該要是 INT_MIN，比 -INT_MAX 再小1) + inf 應該可以改良出有正負的 Inf，但在這裡先不做處理 + 同理其他 bitwise 操作也要做對應的更動 ### 邏輯判斷實作未更新第二個 case 邏輯判斷應該回傳 1 ，但卻回傳了 0。追查到**不等於**的實作，發現邏輯判斷 operator 都會用到 MASK 這個 Macro: ```c #define MASK(n) (((n) > 0) << 4) /* * LSB 4 bits for precision, 2^3, one for sign * MSB 28 bits for integer * LSB all 1 is NaN */ ``` 發現他是如果 n > 0 就設一個 1 在 bit 4，而下面還講解著為何是 4。看起來這是舊的實作所以應該把這個函式連同註解改掉，如果要回傳 1 ，我們就得要把 1 shift 32 bit： ```c #define BOOL(n) ((uint64_t)((n) != 0) << 32) ``` 由於 MASK 這個名字不太清楚，而且這個函式好像只有在邏輯判斷用到，所以改名為 BOOL。而運算的部份需要改成 ```c case OP_NE: /* OK */ return BOOL(!(compare(expr_eval(&e->param.op.args.buf[0]), expr_eval(&e->param.op.args.buf[1])) & (EQUAL))); ``` 其他的邏輯運算也需做對應的更改。 :::danger 改成 `BOOL` 會更難懂，mask 是 bitmask，但 bool 又是什麼？ :notes: jserv ::: :::info 因為我的認知裡面 masking 的作用是要 toggle/set/unset 某些 bit。但是原程式內的 MASK 用途卻沒有拿來作 mask，而比較像是用來產生數值。因為我們的比較運算子的行為很像 c 語言中把非 0 數值轉為 true 或 1 的作法，所以我才想說改成叫 BOOL (Boolean 值)。 > TODO: 思考更明確的變數名稱 :::  --- 如此目前的 test case 都正常了。 # 數值系統實作錯誤首先我嘗試釐清數值系統的格式，觀察 eval 的實作可以最容易看出來： ```c double result = (int) ipt.inte + ((double) ipt.frac / UINT32_MAX); ``` + 可以看到 `inte` 的部份被當成有號整數來看，而 `frac` 的部份則是無號小數。 + 所以當 `inte=12` 而 `frac=0.25` 時，他代表著 12.25。而當 `inte=-12` 而 `frac=0.25` 時，數值代表著 11.75 為了避免是 eval 有寫錯，我也參考了 `expr_parse_number` 的實作，以驗證理解的沒錯。 ## 比較數值大小 ```c 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; } ``` 舊有的實作直接在無號數下比較大小，這樣會造成正數小於負數。因此我們應該先把兩個都**先轉成有號數**數再來比較： ```c static int compare(uint64_t ua, uint64_t ub) { int flags = 0; int64_t a = ua, b = ub; if (a < b) flags |= LOWER; if (a > b) flags |= GREATER; if (a == b) flags |= EQUAL; return flags; } ``` 雖然我們把他分成小數與整數，但是整個 64 bit 一起看的話，事實上相當於一個 64 bit 二補數，再把數值除以 $2^{32}$。因此，我們**可以直接在有號數下做比較**。 ## isNan ```c static int isNan(int x) { return GET_FRAC(x) == GET_FRAC(NAN_INT); } ``` 這個 function 應該是要確認一個數字是不是 nan。但看起來有幾個問題： + 使用的地方都是傳入 uint64_t，但這裡會先被截短為 int + 就算沒有被截短，也應該整個 64 bit 皆比較，而不是只比較 frac。否則有可能 inte 非 0，卻被判定為 NaN 更改後如下： ```c static int isNan(uint64_t x) { return x == NAN_INT; } ``` # 增加新的 Operator ## Sigma ### function 介面設計：原本我想參照 [math expression evaluator](https://github.com/bugwheels94/math-expression-evaluator) 的設計，讓預設的 loop variable 為 `n`，如： ``` "Sigma(1,5,n)" => 15 "Sigma(1,5,1)" => 5 ``` 但在目前的架構下似乎很難實現，因為我們需要一個名為 `n` 的變數，而這些變數都儲存在 `calc()` 內的 `vars` 這個 linked list 內。但是，這個變數清單並不會被傳入 user function 內。如果我們想要修改 `n` 的值，必須從要加總的那個 `struct expr` 中不斷往下找變數，直到找到 `n` 並修改他，但是 `struct expr` 並沒有儲存變數名稱，所以也無從判別是否是 `n` 還是其他變數。後來我參考了 [RinHizakura](https://hackmd.io/@RinHizakura/HyH6mtK_D) 的設計，這個設計把 loop variable 也當成參數傳遞進去，如此自然而然就可以有個 handle 來修改 loop variable 了。這樣做還一個好處，那就是可以在 Sigma 內再放 Sigma 而不會有變數衝突的問題。 ``` "Sigma(n,1,5,n)" => 15 "Sigma(n,1,5,n**2)" => 55 ``` ### 實作 ```c noinline uint64_t user_func_sigma(struct expr_func *f, vec_expr_t args, void *ctx) { if (vec_len(&args) < 4) { return NAN_INT; } struct expr loop_var = vec_nth(&args, 0); int64_t start = expr_eval(&vec_nth(&args, 1)); int64_t end = expr_eval(&vec_nth(&args, 2)); struct expr n_expr = vec_nth(&args, 3); if (loop_var.type != OP_VAR) { return NAN_INT; } uint64_t sum = 0; for (int64_t i = start; i <= end; i += ((uint64_t) 1 << 32)) { *(loop_var.param.var.value) = *(uint64_t *) &i; uint64_t value = expr_eval(&n_expr); if (value == NAN_INT || value == INF_INT) { return value; } sum += value; } return sum; } ``` + 為了讓 user function 可以更容易使用 `OP_*`, `NAN_INT` 與 `INF_INT` 等定義，我把 `enum expr_type` 移動到了 `expression.h`中了。另外也把 `NAN_INT` 與 `INF_INT` 移動到 `fixed-point.c` 中，讓其他檔案可以 include `fixed-point.h` 中的 `extern const` 定義。 :::warning - [x] TODO: 確認是否需要引入 [compensated summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm) :notes: jserv ::: # 定點數、浮點數與 compensated sum ## 浮點數 + [Kahan summation algorithm](https://en.wikipedia.org/wiki/Kahan_summation_algorithm) 又稱 compensated sum, 是為了解決浮點數中作累加時會出現的錯誤。 + 會造成這個錯誤是因為假如要算 $a+b$ 的話，我們必須先把兩個數的指數變成一樣，再將尾數相加。我們以某種 8 bit mantissa 的浮點數做例子： \begin{aligned} a =& 1.00110011 \cdot 2^8\\ b =& 1.01010101 \cdot 2^{0}\\ a+b =& 1.00110011 \cdot 2^8\\ +&0.0000000100110011 \cdot 2^{8}\\ =& 1.00110100\color{red}{00110011} \cdot 2^{8} \end{aligned} 上式中可以看到因為兩個數的指數相差太大導致後面 8 bit 的結果都被 truncate 掉，這就是浮點數累加的誤差來源。反之，如果兩者都是 $2^{8}$ 以這個例子可以完美相加。 + **Kahan summation algorithm** 利用額外的一個變數將紅色的部份儲存起來，再下一次累加時補償回去，再把新的誤差儲存起來，如此重複下去。此作法可降低因為累加過程中的 order 差異過大而出現的誤差。 + 維基百科上也提供了令一種解法，叫作 pairwise summation。此方法的概念就是「**既然 order 差異會造成誤差，那就選擇 order 相近的相加囉**」。假設有 8 個數要相加： \begin{aligned} a_1 + a_2 + \ldots + a_{8} \end{aligned} 照順序加的話最後會出現大數加小數的情況，會出現誤差。但如果我們這樣加： \begin{aligned} \left((a_1 + a_2) + (a_1 + a_2)\right) + \left((a_1 + a_2) + (a_1 + a_2)\right) \end{aligned} 每個加號左右的數 order 都差不多（都經歷過差不多數目的累加），這樣也可以降低因相加所產生的誤差。 + **總結浮點數的部份**，可以看到浮點數是 non-associative 的，也就是相加的順序會影響相加的結果。 ## 定點數 + 在我上面的分析中，我提到我們目前使用的數值格式相當於 64 bit 的二補數再除以 $2^{32}$。而二補數一個重要的性質是他是一個 **abelian group**，也就是加法的順序永遠可以調換而不會影響相加結果。也就是，定點數上並不需要引進 compensated sum 這類技巧。 ## 實驗 ### 測試程式 + Sigma 的補償實作 ```c for (int64_t i = start; i <= end; i += ((uint64_t) 1 << 32)) { *(loop_var.param.var.value) = *(uint64_t *) &i; uint64_t value = expr_eval(&n_expr); if (value == NAN_INT || value == INF_INT) { return value; } uint64_t y = value - c; uint64_t t = sum + y; c = (t - sum) - y; sum = t; } ``` + 浮點數 (naive) ```c float sum=0; for (int j = 0; j < s; j++){ sum += (1.0f + 1.0f/1024); } printf("%.10f\n", sum); ``` + 浮點數 (compensated) ```c float sum=0.0, c=0.0; for (int j = 0; j < s; j++){ float y = (1.0f + 1.0f/1024) - c; float t = sum + y; c = (t - sum) - y; sum = t; } assert(sum == (1.0f + 1.0f/1024)*s); printf("%.10f\n", sum); ``` ### 結果 1. 這個實驗是將 $\frac{1}{3^n}$ 加總 $3^n$ 次。橫軸為 $n$，如果結果為一就代表是無誤差。 + 可以看到紅色線，浮點數再沒有做補償的情況下，誤差相當大。而做了補償之後（綠線），表現最好。至於 kcalc 定點數的實作無論有無補償，結果都一樣，如同預測的。 + 但是定點數兩個實作中都存在誤差，這是因為 $\frac{1}{3^n}$ 無法用有限的二進位表示，所以一定有 round-off error 存在，而這個誤差會在加總愈多次之後逐漸被放大。 ![](https://i.imgur.com/aQuCoG0.png) 2. 第二個實驗將 $1+\frac{1}{1024}$ 加總 $2^n$ 次。 y 軸為結果除以正解的值， 1 代表結果準確。 + 可以看到在這個測試下，浮點數在 $n=14$ 之後就開始有錯誤，這是因為浮點數的 mantissa 為 23 bit。而 $\frac{1}{1024}$ 在第 10 個 bit。因此大約從那裡開始就會發生 $\frac{1}{1024}$ 被遺忘的情況。 + 定點數因為沒有了 round-off error，得到 exact 的結果。 + 補償後的浮點數也可以得到 exact 的結果。因為這個數字本來就可以被浮點數 exact 表達。 ![](https://i.imgur.com/IfOKtWg.png) --- 結論： **定點數的確不需要作 compensated sum**

Read more

2020q3 Homework1 (quiz1)

2020q3 Homework3 (quiz3)

2020q3 Homework5 (quiz5)

2020q3 Homework5 (render)