# 2016q3 Homework1 (compute-pi)
contributed by <`HaoTse`>

---

## SIMD指令集
- AVX (Advanced Vector Extensions) 是 Intel 一套用來作 Single Instruction Multiple Data 的指令集。
- 我的電腦 CPU 為 intel i5-4200U
- 使用時需`#include <immintrin.h>`，並在 Makefile 中的編譯選項加入`-mavx`
- [Intel Intrinsics Guide](https://software.intel.com/sites/landingpage/IntrinsicsGuide/#expand=1,4629)

---

## OpenMP
在 raytracing 中使用過的 OpenMP 也可以應用在 compute-pi 當中。

不過在 compute-pi 中每次迭代都會讀取 pi 值並更新 pi 值，不同次的迴圈中並不相互獨立，有兩個方法可以避免：
1. 可以使用 atmoic、critical，防止變數同時被多個執行緒修改，但速度就會慢上許多，因為一次只允許一個執行續進行存取。
2. 使用 reduction，語法為 `reduction(<op>:<variable>)`。
> 作業採用方法2 [name=鄭皓澤]

---

## Wall-clock time vs. CPU time

----

### Wall-clock time
- 顧名思義就是掛鐘時間，也就是現實世界中實際經過的時間，是由 kernel 裡的 xtime 來紀錄，系統每次啟動時會先從設備上的 RTC 上讀入 xtime。
- xtime 精確度到**微秒**
- [xtime vs. jiffies](http://www.coctec.com/docs/linux/show-post-186188.html)
- System time：指系統上的時間，開機時它會讀取 RTC 來設定，可能過程中還會有時區換算等之類的設置。一般說來 System time 就是我們執行 date 命令看到的時間。
- Wall-clock time 可能會和 NTP 伺服器、時區、日光節約時間同步或使用者自己調整。所以「通常」不建議拿來量測程式時間，因為它不是一個穩定的時間，用專業一點的用語講，**Wall-clock time 不一定是單調遞增 (monotonic)**。

----

### CPU time
- 程序在 CPU 上面運行消耗 (佔用) 的時間，**clock()** 就是很典型用來計算 CPU time 的時間函式。
- **注意**. 如果有多條 threads，那 CPU time 算的是「每條 thread 的使用時間加總」，所以 CPU time 可能比 Wall-clock time 還大。
- 執行時間除了 CPU time 以外，可能還包括 I/O time、 communication channel delay、synchronization overhead…等等。

----

### time 指令
今日的 UNIX 系統內事實上存在兩個 time 指令，一個是系統所提供的 /usr/bin/time，這個 time 指令是 System V 版本所提供的工具，這是為了 Bourne shell 使用者所寫的工具；而另一個則是 C shell 內建指令 time。
> 因此 `$ make check` 輸出格式與直接 `$ time ./程式名稱` 不同。[name=鄭皓澤]

1. real time : 也就是 Wall-clock time，當然若有其他程式同時在執行，必定會影響到。
2. user time : 表示程式在 user mode 佔用所有的 CPU time 總和。
3. sys time : 表示程式在 kernel mode 佔用所有的 CPU time 總和 ( 直接或間接的系統呼叫 )。

----

### timer 比較

#### clock()
- `#include <time.h>`
- 回傳程式執行後佔用的 CPU 時間，單位為 **tick** 。

#### time()
- `#include <time.h>`
- 回傳日曆時間，也就是自 1970 年 1 月 1 日後經過的總秒數。

#### clock_gettime()
- `#include <time.h>`
- 編譯時連結 `librt.so`
gcc編譯時加入 `-Wl,-lrt` 選項，或連結時加入 `-lrt` 選項。
- 函式原型

	int clock_gettime(clockid_t clk_id, struct timespec *tp);

`struct timespec* tp` 是函式回傳結果，`struct timespec` 宣告如下

	struct timespec {
        time_t   tv_sec;        /* seconds */
        long     tv_nsec;       /* nanoseconds */
	};

- 精確到奈秒等級。

#### gettimeofday()
- `#include<sys/time.h>`
- 函式原型

	int gettimeofday(struct  timeval*tv,struct  timezone *tz )

---

## 分析
首先直接將 Makefile gencsv 中的 for loop 改為 `seq 1000 5000 1000000`，得到下圖


這張圖長的有點醜，這裡研究了 Libre Office 很久，最後直接使用 gnuplot 畫出來，如下圖
> 之後才發現 Libre Office 中的 Chart Type 要選 XY Scatter。[name=鄭皓澤]

> 注意. gnuplot 讀入 .csv檔時，","後面要加上空白或將他改為空白。[name=鄭皓澤]

----

### 信賴區間

在 trace `benchmark_clock_gettime.c` 時，發現裡面有許多 for loop 去跑計算 pi 的 function，回想起觀看[廖健富學長筆記](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I)時老師提到的**信賴區間**，將這裡的 for loop 視為取樣，並將信賴區間紀錄起來。

計算信賴區間函式如下

double compute_ci(double *min, double *max, double data[SAMPLE_SIZE])
{
    double mean = 0.0;
    double std_dev = 0.0; //standard deviation
    double std_err;
    int i = 0;

    //calculate mean value
    for(i = 0; i < SAMPLE_SIZE; i++)
    {
        mean += data[i];
    }
    mean /= (double)SAMPLE_SIZE;

    //calculate standard deviation
    for(i = 0; i < SAMPLE_SIZE; i++)
    {
        std_dev += (data[i] - mean) * (data[i] - mean);
    }
    std_dev = sqrt(std_dev / (double)SAMPLE_SIZE);

    //calculate standard deviation
    std_err = std_dev / sqrt((double)SAMPLE_SIZE);

    *min = mean - (1.96 * std_err);
    *max = mean + (1.96 * std_err);

    return mean;
}

> 如何使用信賴區間使得數據更好看還在研究了不少時間，當初機統應該好好學的QQ，下圖為使用信賴區間跑的成果

不過用 thread 的兩個版本還是有凸出，原因還在探討中。

----

### error rate
目前只能看出花費時間的比較，無法比較計算出來的 pi 的正確性，因此決定再畫一張圖做比較。首先定義 pi 在 double 型態下的正確值，在`<math.h>`中有定義了M_PI可以用做比較。加入以下 macro

#define compute_error(pi) ( (pi > M_PI) ? pi - M_PI : M_PI - pi)

- 輸出結果
這的結果與我預期的不一樣，我原本認為用同一個算式的 error rate 應該都相同，但在這張圖中我卻發現使用 AVX + unroll looping 的版本在 1000、11000、21000...中的數字差異甚大，原本我認為是我更動到了程式碼，一直找不到原因下，我將原本的程式碼 clone 下來，將 `time_test.c` 中的 N 改為 1000，`$ make check` 後發現 AVX + unroll looping 版本仍然輸出差異比別人，至於為什麼會這樣還在思考中。

---

## Leibniz

原本的 `baseline` 算法為
$$
\pi = 4\int_0^1\frac{1}{1+x^2}dx
$$

### Leibniz 級數

$$
\pi = 4\sum_{i=1}^\infty\frac{(-1)^n}{2n+1}
$$
- 證明
$\frac{1}{1+x^2} = 1 - x^2 + x^4 - \dots + (-1)^nx^{2n} + \frac{(-1)^{n+1}x^{2n+2}}{x^2+1}$

等號兩邊做定積分

$\frac {\pi}{4} = \int_0^1\frac{1}{1+x^2}dx = 1 - \frac{1}{3} + \frac{1}{5} - \dots + \frac {(-1)^n}{2n+1} + (-1)^{n+1}\int_0^1\frac{x^{2n+2}}{x^2+1}dx$

當 $n\to \infty$ 時，最後一項趨近於0

$\int_0^1\frac{x^{2n+2}}{x^2+1}dx < \int_0^1x^{2n+2}dx = \frac{1}{2n+3} \to 0$

$\frac {\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \dots + \frac {(-1)^n}{2n+1} = \sum_{i=1}^\infty\frac{(-1)^n}{2n+1}$

故得證$\pi = 4\sum_{i=1}^\infty\frac{(-1)^n}{2n+1}$

----

### 結果分析
將 Leibniz 公式轉成程式碼如下

double compute_pi_leibniz(size_t N)
{
    double pi = 0.0;
    for(size_t i = 0; i < N; i++)
    {
        int sign = (i % 2) ? -1 : 1;
        pi += (sign / (2.0 * (double)i +1.0));
    }

    return pi * 4.0;
}

得到結果
- runtime
這裡看出 Leibniz 的速度比 baseline 還快，我猜測是因為 Leibniz 的乘除運算較少。
- error_rate

---

## 參考資料
- [SIMD wiki](https://zh.wikipedia.org/wiki/%E5%8D%95%E6%8C%87%E4%BB%A4%E6%B5%81%E5%A4%9A%E6%95%B0%E6%8D%AE%E6%B5%81)
- [time 指令輸出格式](http://linux.vbird.org/linux_basic/0320bash/csh/no3-8-07.html)
- timer
- [c timer 整理](http://edisonx.pixnet.net/blog/post/52113788-%5Bc%5D-%E8%A8%88%E6%99%82%E5%99%A8%E6%95%B4%E7%90%86)
- [clock_gettime](https://starforcefield.wordpress.com/2012/08/06/c%E8%AA%9E%E8%A8%80%EF%BC%9A%E4%BD%BF%E7%94%A8clock_gettime%E8%A8%88%E7%AE%97%E7%A8%8B%E5%BC%8F%E7%A2%BC%E7%9A%84%E6%99%82%E9%96%93%E9%9C%80%E6%B1%82/)
- [gettimeofday](http://nosleep.pixnet.net/blog/post/205120138-%E7%A8%8B%E5%BC%8F%E9%96%8B%E7%99%BC-%7C-%5Blinux%5D%5Bc%5D-%E4%BD%BF%E7%94%A8-gettimeofday()-%E5%87%BD%E5%BC%8F%E8%A8%88%E7%AE%97)
- [gnuplot](http://www.codelast.com/%E5%8E%9F%E5%88%9B-gnuplot%E8%B0%83%E6%95%99%E8%AE%B0/)
- [standard error](https://zh.wikipedia.org/wiki/%E6%A0%87%E5%87%86%E8%AF%AF)
- [LaTex](http://mohu.org/info/symbols/symbols.htm)
- [信賴區間的迷思](http://yenchic-blog.logdown.com/posts/177267-confidence-interval-of-myth)
- [廖健富學長筆記](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I)

###### tags: `HaoTse` `sysprog21` `compute-pi` `SIMD AVX` `timer` `Leibniz`