2020q3 Homework6 (kcalc)

contributed by < sammer1107 >

tags: `進階電腦系統理論與實作`

作業說明

完善測試機制

原本的 test case 跑完也不確定結果到底正不正確，所以我決定先從改善 test 機制開始。























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

使用範例：

# multiply
test_op '6*7' 42

輸出:

Testing  6*7 ... PASSED
== 42.000000

失敗的 case

根據原先註解預期的輸出，我發現有幾個 test case 目前會失敗：

Testing  (3%0)|0 ... FAILED
Got NAN_INT instead of 0.000000

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 的程式碼做更動：

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:

#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：

#define BOOL(n) ((uint64_t)((n) != 0) << 32)

由於 MASK 這個名字不太清楚，而且這個函式好像只有在邏輯判斷用到，所以改名為 BOOL。而運算的部份需要改成

case OP_NE: /* OK */
    return BOOL(!(compare(expr_eval(&e->param.op.args.buf[0]),
                          expr_eval(&e->param.op.args.buf[1])) &
                      (EQUAL)));

其他的邏輯運算也需做對應的更改。

改成 BOOL 會更難懂，mask 是 bitmask，但 bool 又是什麼？

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

jserv

因為我的認知裡面 masking 的作用是要 toggle/set/unset 某些 bit。但是原程式內的 MASK 用途卻沒有拿來作 mask，而比較像是用來產生數值。因為我們的比較運算子的行為很像 c 語言中把非 0 數值轉為 true 或 1 的作法，所以我才想說改成叫 BOOL (Boolean 值)。

TODO: 思考更明確的變數名稱

如此目前的 test case 都正常了。

數值系統實作錯誤

首先我嘗試釐清數值系統的格式，觀察 eval 的實作可以最容易看出來：

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 的實作，以驗證理解的沒錯。

比較數值大小

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

舊有的實作直接在無號數下比較大小，這樣會造成正數小於負數。因此我們應該先把兩個都先轉成有號數數再來比較：

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

static int isNan(int x)
{
    return GET_FRAC(x) == GET_FRAC(NAN_INT);
}

這個 function 應該是要確認一個數字是不是 nan。但看起來有幾個問題：

使用的地方都是傳入 uint64_t，但這裡會先被截短為 int
就算沒有被截短，也應該整個 64 bit 皆比較，而不是只比較 frac。否則有可能 inte 非 0，卻被判定為 NaN

更改後如下：

static int isNan(uint64_t x)
{
    return x == NAN_INT;
}

增加新的 Operator

Sigma

function 介面設計：

原本我想參照 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 的設計，這個設計把 loop variable 也當成參數傳遞進去，如此自然而然就可以有個 handle 來修改 loop variable 了。這樣做還一個好處，那就是可以在 Sigma 內再放 Sigma 而不會有變數衝突的問題。

"Sigma(n,1,5,n)" => 15
"Sigma(n,1,5,n**2)" => 55

實作

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 定義。

TODO: 確認是否需要引入 compensated summation
Image Not Showing Possible Reasons
- The image file may be corrupted
- The server hosting the image is unavailable
- The image path is incorrect
- The image format is not supported
Learn More →
jserv

定點數、浮點數與 compensated sum

浮點數

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 00110011 \cdot 2^{8} \end{aligned}$
上式中可以看到因為兩個數的指數相差太大導致後面 8 bit 的結果都被 truncate 掉，這就是浮點數累加的誤差來源。反之，如果兩者都是
$2^{8}$ 以這個例子可以完美相加。
Kahan summation algorithm 利用額外的一個變數將紅色的部份儲存起來，再下一次累加時補償回去，再把新的誤差儲存起來，如此重複下去。此作法可降低因為累加過程中的 order 差異過大而出現的誤差。
維基百科上也提供了令一種解法，叫作 pairwise summation。此方法的概念就是「既然 order 差異會造成誤差，那就選擇 order 相近的相加囉」。假設有 8 個數要相加：

$\begin{array}{r} a_{1} + a_{2} + \dots + a_{8} \end{array}$
照順序加的話最後會出現大數加小數的情況，會出現誤差。但如果我們這樣加：

$\begin{array}{r} ((a_{1} + a_{2}) + (a_{1} + a_{2})) + ((a_{1} + a_{2}) + (a_{1} + a_{2})) \end{array}$
每個加號左右的數 order 都差不多（都經歷過差不多數目的累加），這樣也可以降低因相加所產生的誤差。
總結浮點數的部份，可以看到浮點數是 non-associative 的，也就是相加的順序會影響相加的結果。

定點數

在我上面的分析中，我提到我們目前使用的數值格式相當於 64 bit 的二補數再除以
$2^{32}$ 。而二補數一個重要的性質是他是一個 abelian group，也就是加法的順序永遠可以調換而不會影響相加結果。也就是，定點數上並不需要引進 compensated sum 這類技巧。

實驗

測試程式

Sigma 的補償實作

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)

float sum=0;
for (int j = 0; j < s; j++){
    sum += (1.0f + 1.0f/1024);
}
printf("%.10f\n", sum);

浮點數 (compensated)

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

結果

這個實驗是將
$\frac{1}{3^{n}}$ 加總
$3^{n}$ 次。橫軸為
$n$ ，如果結果為一就代表是無誤差。
- 可以看到紅色線，浮點數再沒有做補償的情況下，誤差相當大。而做了補償之後（綠線），表現最好。至於 kcalc 定點數的實作無論有無補償，結果都一樣，如同預測的。
- 但是定點數兩個實作中都存在誤差，這是因為
  $\frac{1}{3^{n}}$ 無法用有限的二進位表示，所以一定有 round-off error 存在，而這個誤差會在加總愈多次之後逐漸被放大。
第二個實驗將
$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 表達。

結論： 定點數的確不需要作 compensated sum