# 2020q3 Homework1 (lab0)
contributed by < `OscarShiang` >

## 測試環境

$ cat /etc/os-release 
NAME="Ubuntu"
VERSION="20.04 LTS (Focal Fossa)"
$ cat /proc/version
Linux version 5.4.0-47-generic (buildd@lcy01-amd64-014) 
(gcc version 9.3.0 (Ubuntu 9.3.0-10ubuntu2))

## 作業要求
* [x] ==詳細閱讀 [C Programming Lab](https://www.cs.cmu.edu/afs/cs/academic/class/15213-s20/www/labs/cprogramminglab.pdf)== ，依據指示著手修改 `queue.[ch]` 和連帶的檔案，測試後用 Git 管理各項修改，記得也該實作 `q_sort` 函式。
* [x] 運用 Valgrind 排除 `qtest` 實作的記憶體錯誤，並透過 Massif 視覺化 "simulation" 過程中的記憶體使用量，需要設計對應的實驗
* [x] 研讀論文 [Dude, is my code constant time?](https://eprint.iacr.org/2016/1123.pdf)，解釋本程式的 "simulation" 模式是如何透過以實驗而非理論分析，達到驗證時間複雜度，需要解釋 [Student's t-distribution](https://en.wikipedia.org/wiki/Student%27s_t-distribution) 及程式實作的原理。注意：現有實作存在若干致命缺陷，請討論並提出解決方案;
<s>
* [x] 參照 [antirez/linenoise](https://github.com/antirez/linenoise)，強化 `qtest` 命令列功能，使其支援單行或多行編輯、歷史命令處理、命令自動補完 (例如輸入 `h` 之後按下 Tab 按鍵，自動接續補上 `elp`，成為完整的 `help` 命令)。除了整合程式碼外，應當說明 [antirez/linenoise](https://github.com/antirez/linenoise) 的運作原理，注意到 [termios](http://man7.org/linux/man-pages/man3/termios.3.html) 的運用
* 思考如何預設開啟 AddressSanitizer，卻又能滿足效能測試所需 (可能需要兩種以上的 build target)
* 指出現有程式的缺陷 (提示: 和 [RIO 套件](http://csapp.cs.cmu.edu/2e/ch10-preview.pdf) 有關)，嘗試強化並提交 pull request
</s>
> 已整合在 [sysprog21/lab0-c](https://github.com/sysprog21/lab0-c)

## 實作 queue.[ch]

:::info
:pencil: 作業筆記請參考 [2020q1 Homework1 (lab0)](/2b4MZK2DTZqMoe69dcpPhQ)

:computer: 實作的程式碼則在 [OscarShiang/lab0-c](https://github.com/OscarShiang/lab0-c) 的 `master` 分支
:::

## 導入 [antirez/linenoise](https://github.com/antirez/linenoise)

:::success
:desktop_computer: 導入 [antirez/linenoise](https://github.com/antirez/linenoise) 的實作可以參考 [OscarShiang/lab0-c](https://github.com/OscarShiang/lab0-c/tree/linux2020) 的 `linux2020` 分支
:::

## 研讀論文 〈[Dude, is my code constant time?](https://eprint.iacr.org/2016/1123.pdf)〉

使用〈[Dude, is my code constant time?](https://eprint.iacr.org/2016/1123.pdf)〉中所實作出的工具 `dudect` 來檢測演算法的時間複雜度是否屬於 $O(1)$ 的特點有以下幾個

- **經由實際測量而非理論來推導時間複雜度**
`dudect` 以重複測量程式執行時間建立樣本，利用 Welch's t-test 來分析樣本時間的分佈情形，進而推論程式執行的時間複雜度是否為常數時間。而不是通過公式推導的方式來驗證，進而降低了在實作上判斷時間複雜度為 $O(1)$ 的難度。
- **不會受到硬體架構差異的影響**
因為 `dudect` 在實作上是透過統計模型推斷時間複雜度是否為常數時間，所以不會因為使用不同的硬體進行檢測而有不同的結果。

### Welch's t-test

在判斷兩個 class 的樣本平均是否具有同質性，這邊使用到的是 Welch's t-test，在 [Welch's t-test - Wikipedia](https://en.wikipedia.org/wiki/Welch%27s_t-test) 中提到的：

> Welch's t-test is an adaptation of Student's t-test, and is more reliable when the two samples have unequal variances and/or unequal sample sizes.

在假設兩組樣本都具屬於常態分佈的前提下，可以使用 Welch's t-test 來檢測兩者的平均是否等值

Welch's t-test 所定義的 t value 如下

$$
t = \frac{\overline{x}_1 - \overline{x}_2}{\sqrt{\frac{s^2_1}{n_1} + \frac{s^2_2}{n_2}}}
$$

### 計算平均與變異數的時使用到的 Welford method

根據 [Algorithms for calculating variance - Wikipedia](https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm) 對於 Welford's online algorithm 的解釋：
> It is often useful to be able to **compute the variance in a single pass, inspecting each value $x_{i}$ only once**; for example, when the data are being collected without enough storage to keep all the values, or when costs of memory access dominate those of computation.

可以得知 Welford's online algorithm 與直接計算平均與變異數的差別在於不需要先將所有數值加總後才進行計算，而是逐次將數值加入並計算對應的 $\overline{x}_n$ 與 $s^2_n$。

其手法的原理大致如下：

假設現在我們有以下共 $n$ 個的樣本

$$
x_1, x_2, x_3, \dots, x_{n}
$$

一般我們計算平均的方式是

$$
\overline{x}_n = \frac{1}{n}\sum_{i = 1}^{n}x_i
$$

若我們要在此計算出來的平均上再加入新的樣本 $x_k$，則須這樣重新計算平均的數值

$$
\overline{x}_{n + 1} = \frac{n \times \overline{x}_n + x_k}{n + 1}
$$

對 $n$ 做 $n' = n + 1$ 的變換與化簡後則可以寫成

$$
\begin{align}
\overline{x}_{n'} &= \frac{(n' - 1) \ \overline{x}_{n' - 1} + x_k}{n'} \\\\
&= \overline{x}_{n' - 1} + \frac{\overline{x}_{n' - 1} - x_k}{n'}
\end{align}
$$

而這樣的做法也可以避免在加總所有實驗數據的時候發生 Overflow 。

但是若直接利用這樣的方式計算變異數的話，則可能會造成部分問題，因為若使用該方法計算 n - 1 個樣本的變異數 $\sigma^2_{n - 1}$ 來計算 n 的樣本的 $\sigma^2_n$ 時

$$
\begin{align}
s^2_n &= \frac{n - 2}{n - 1}s^2_{n - 1} + \frac{(x_n - \overline{x}_{n - 1})^2}{n} \\\\
&= s^2_{n - 1} + \frac{(x_n - \overline{x}_{n - 1})^2}{n} - \frac{s^2_{n - 1}}{n - 1}, \ \ \ \ (n > 1)
\end{align}
$$

而在上式在使用電腦計算時發生 Underflow 的情況而產生 Numerical Instability 的情況。

為了免除這樣的問題，所以 Welford method 在操作上定義差值平方的總合為 $M_{2, n}$

$$
M_{2, n} = \sum_{i = 1}^n (x_i - \overline{x}_n)^2
$$

並將上述 $(*)$ 處的等式兩側同乘 $n$ 後則變成

$$
M_{2, n} = M_{2, {n - 1}} + (x_n - \overline{x}_{n - 1})(x_n - \overline{x}_n)
$$

待我們計算出 $M_{2, n - 1}$ 後再換算回變異數即可

$$
s^2_n = \frac{M_{2, n}}{n}
$$

### 在 lab0-c 中的實作

根據上述論文的理論，在 `lab0-c/dudect` 中對應的實作流程如下：

1. **生成隨機字串**
因為在我們測試的 list 中是使用字串作為資料存放，所以在測試的一開始會先透過在`random.c` 中定義的 `randombyte` 來生成 10000 長度為 7 byte 的字串。

2. **亂數產生 list 長度**
在我們的測試程式中，預設會重複測試 10000 組資料，所以會使用前一個步驟中所使用的 `randombyte` 函式來產生每組測試資料的 list 長度。並利用 `randombit` 隨機將該組資料分配到 2 個 class 之中。

3. **依序初始化各個測試的 list 並測量函式的執行 cycle 數**
在每次進行測試的時候，程式會先根據上個階段產生的 list 長度生成一個隨機的 list，接著執行測試的函式（`q_size` 或是 `q_insert_tail`）進行測試。透過呼叫 Intel x86-64 架構的指令 `rdtsc` 取得 time stamp counter 的數值，以此來計算執行的 cycle 數。

4. **計算 $\overline{x},\ M_{2, n}$，並使用 Welch's t-test 檢驗兩個 class 是否有類似的分佈**
將我們在第 3 步測得的資料依照在第 2 步中分配到的 class 分別計算 $\overline{x}$ 與 $M_{2, n}$ ，使用上述的 Welford's online algorithm 來計算兩個 class 個別的 $\overline{x}$ 與 $M_{2, n}$。最後利用前一步取得的資料計算出 t 值

在 `lab0-c/dudect` 的實作中，若 $t < 10$ 則表示兩組資料具有相同的分佈，也就可以合理的推論我們所測量的函式執行時間應為常數時間，也就是 $O(1)$ 的複雜度。

### 實作上的問題

#### 使用 `rdtsc` 指令計算 cycle 數所造成的問題

在〈[Using the RDTSC Instruction for Performance Monitoring](https://www.ccsl.carleton.ca/~jamuir/rdtscpm1.pdf)〉這篇文章中有提到因為 Intel Pentium 系列的處理器有支援 out-of-order execution ，也就是亂序執行的功能，所以在使用該指令取得 cycle count 的時候可能我們想要測量的函式還沒有被執行到，而導致無效的測量。

而在 [多核时代不宜再用 x86 的 RDTSC 指令测试指令周期和时间](https://blog.csdn.net/solstice/article/details/5196544) 中也有提到使用 `rdtsc` 來計時或是 benchmarking 可能會也以下幾種潛在的問題：
- 因為在多核心的實驗環境下我們不能保證我們的程式在測試的過程中不會發生資料遷移的情形，並且每個 CPU 的 Time Stamp Counter 的數值可能會有誤差，進而導致 benchmarking 的結果可能會不夠準確或有誤。
- CPU 的頻率可能會有變動，像是在省電模式時頻率會較低，高效能模式時 CPU 的頻率會開到比較大。而這樣的變化會讓 cycle counter 記錄到的 cycle 刻度會不同，造成量測精度的問題。

#### 使用 `dropsize` 忽略部分資料

在 `dudect/constant.c` 中可以看到工具中有設定 `dropsize` 的數值為 20，但在套用統計模型的時候卻沒有看到 `dudect` 有針對測試的數據做後處理，而是直接將前 20 筆的資料捨棄掉，計算第 20 ~ 9999 的資料。

#### t-test Threshold 的選擇

在該篇論文中第三部分 Result 的 B. Memory Comparison 中可以看到在使用 Timing Leakage Detection 進行測試的時候，他們所使用的 t value 的 threshold 為 4.5 但是實際上在 `lab0-c/dudect` 中看的 threshold 則為 10

### 加入 pre-processing

在目前的實作中並沒有在論文中提到的後處理，而是單純以一階 t-test 來判斷是否有 timing leakage

而在實際執行時有時會遭遇這樣的情形：

$ ./dut_size
meas:    0.01 M, max t:  +11.19, max tau: 1.12e-01, (5/tau)^2: 2.00e+03.
meas:    0.01 M, max t:  +11.29, max tau: 1.12e-01, (5/tau)^2: 1.99e+03.
meas:    0.01 M, max t:  +11.25, max tau: 1.11e-01, (5/tau)^2: 2.02e+03.

...

meas:    0.01 M, max t:  +12.41, max tau: 1.13e-01, (5/tau)^2: 1.96e+03.
meas:    0.01 M, max t:  +12.45, max tau: 1.13e-01, (5/tau)^2: 1.97e+03.
meas:    0.01 M, max t:  +12.60, max tau: 1.14e-01, (5/tau)^2: 1.94e+03.
meas:    0.01 M, max t:  +12.42, max tau: 1.11e-01, (5/tau)^2: 2.02e+03.
meas:    0.01 M, max t:  +12.37, max tau: 1.10e-01, (5/tau)^2: 2.05e+03.
meas:    0.01 M, max t:  +12.36, max tau: 1.10e-01, (5/tau)^2: 2.07e+03.
meas:    0.01 M, max t:  +12.41, max tau: 1.10e-01, (5/tau)^2: 2.07e+03.
meas:    0.01 M, max t:  +12.48, max tau: 1.10e-01, (5/tau)^2: 2.06e+03.
meas:    0.01 M, max t:  +12.45, max tau: 1.09e-01, (5/tau)^2: 2.09e+03.
meas:    0.01 M, max t:  +12.76, max tau: 1.11e-01, (5/tau)^2: 2.01e+03.

...

meas:    0.03 M, max t:  +19.33, max tau: 1.16e-01, (5/tau)^2: 1.86e+03.
meas:    0.03 M, max t:  +19.34, max tau: 1.16e-01, (5/tau)^2: 1.86e+03.
meas:    0.03 M, max t:  +19.37, max tau: 1.16e-01, (5/tau)^2: 1.86e+03.
meas:    0.03 M, max t:  +19.41, max tau: 1.16e-01, (5/tau)^2: 1.86e+03.
meas:    0.03 M, max t:  +19.50, max tau: 1.16e-01, (5/tau)^2: 1.85e+03.
meas:    0.03 M, max t:  +19.49, max tau: 1.16e-01, (5/tau)^2: 1.86e+03.
meas:    0.03 M, max t:   +0.33, max tau: 1.96e-03, (5/tau)^2: 6.49e+06.
meas:    0.03 M, max t:   +0.33, max tau: 1.98e-03, (5/tau)^2: 6.41e+06.
meas:    0.03 M, max t:   +0.34, max tau: 1.99e-03, (5/tau)^2: 6.33e+06.
meas:    0.03 M, max t:   +0.34, max tau: 2.03e-03, (5/tau)^2: 6.04e+06.
^C

在演算法確定為 $O(1)$ 的情況下卻出現 $t > 10$ 的情況，而被判定為非常數時間。

推測其發生的狀況可能是因為 content switch 或是 interrupt 所造成，所以先導入前處理的方式來針對異常資料做剔除的動作。

在論文中提到兩種前處理的手法：
- 剔除超出百分位數的資料
- 使用高階 t-test 檢測

而在 dudect 中也有對應的實作，我將 lab0-c 的程式碼進行修改以導入上述兩個手法，但是在執行的時候發現在使用 percentile 剔除異常資料後，即使 `q_size` 已是常數時間的實作仍會被判定為 not constant time 的狀況。

因此我將 class 0 與 class 1 的測量結果，也就是 `exec_time` 的資料提取出來，從圖表進行分析：

#### Class 0 (fixed class) 測量 `q_size` 的結果


可以看到大部分的測量結果是落在 0 ~ 100 個 cycle 的區間中，但也可以看到有部份的資料是超過甚至遠大於 100 個 clock cycle 的情況，在 10000 筆資料中，有 154 筆的資料大於 100 個 cycle。

#### Class 1 (random class) 測量 `q_size` 的結果


在本圖中可以看到在 $O(1)$ 的實作下 fixed class 與 random class 的分佈其實是十分類似的。但在結果中可以看到 random class 的測量結果較 fixed class 的結果集中在約 15 ~ 35 個 cycle 的區間中。而在本結果中超過 100 個 cycle 的資料共有 156 筆。

為了比較 $O(1)$ 與 $O(n)$ 兩個實作的差異，我將 `q_size` 的函式修改為以下 $O(n)$ 的實作：

int q_size(queue_t *q) 
{
    if (!q)
        return 0;

    int size = 0;
    list_ele_t *curr = q->head;
    while (curr) {
        size++;
        curr = curr->next;
    }   
    return size;
}

並將 fixed class 與 random class 的執行結果繪成下圖：

#### Class 0 (fixed class) 測量非常數實作之 `q_size` 執行時間結果


可以看到上圖其實與上述常數時間實作的結果大致相同，因為 fixed class 的測試都是對一個長度為 0 的 list 呼叫 `q_size`，所以在測量時間上與常數時間的分佈不會有太大的差異。

#### Class 1 (random class) 測量非常數實作之 `q_size` 執行時間結果


從這個上圖的結果可以看到 $O(n)$ 版本的 `q_size` 在測量數值上不管是數量級或是 cycle 數的分佈上都與 fixed class 的版本或是 $O(1)$ 的情況相當的不同，所以在使用 t-test 來檢測的時候就可以很容易的得到 50 以上的 t value 進而推斷這樣的實作不是 $O(1)$ 的實作。

### 改用 `clock_gettime` 進行測量

由於上述提到使用 Intel 架構的 `rdtsc` 指令可能會有測量誤差的問題，因此我將 `cpucycle.h` 中用來測量的函式改以 `time.h` 的 `clock_gettime` 來進行測量。

>　原本想要用 `sys/time.h` 的 `gettimeofday` 來進行測量，但是因為 `q_size` 的執行時間太短，無法使用 microsecond 的精度來進行測量，因此改以 `clock_gettime` 來實作。 [color= orange][name= Oscar]

#### Class 0 (fixed class) $O(1)$ 實作之 `q_size` 測量結果


#### Class 1 (random class) $O(1)$ 實作之 `q_size` 測量結果


## 參考資料

1. [2020 年秋季 進階電腦系統理論與實作課程作業 —— lab0](https://hackmd.io/W0kSCHBkQcqCVDEw5NOgEA?both)
2. [OscarShiang/lab0-c](https://github.com/oscarshiang/lab0-c)
3. [Dude, is my code constant time?](https://eprint.iacr.org/2016/1123.pdf)
4. [Algorithms for calculating variance - Wikipedia](https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm)
5. [Using the RDTSC Instruction for Performance Monitoring](https://www.ccsl.carleton.ca/~jamuir/rdtscpm1.pdf)
6. [rdtsc指令，测量程序的运行速度](http://blog.chinaunix.net/uid-24774106-id-2779245.html)
7. [多核时代不宜再用 x86 的 RDTSC 指令测试指令周期和时间](https://blog.csdn.net/solstice/article/details/5196544)
8. [x86: unify/rewrite SMP TSC sync code](https://lwn.net/Articles/211051/)

###### tags: `sysprog2020`