# 2017q1 Homework3 (compute-pi)
contributed by <`xdennisx`>
### Reviewed by `baimao8437`
* 每次新增 method 都有新的分析圖形，但不能只丟圖上來就好，應稍微解釋圖片中的現象，如: Nilakantha vs Leibniz & Euler 的 error rate 圖形看起來一模一樣，沒有其他文字解釋的話看不出差異。
* 對於資料的抖動的情況，可以嘗試使用95%信賴區間對資料進行篩選。
* 由於此作業對於計算 $\pi$ 會使用多種 method ，使得每新增一種，就得在 makefile 中多寫很多行，可嘗試簡化 makefile 已提升測試效率，可參考[神之修改](https://github.com/0140454/compute-pi/commit/89ab12a144d0ce6d2d3c498e50277fe81c42b8d4)。
---
## 原始圖形
將於本 result_clock_gettime.csv 的圖形用 gunplot 輸出
OpneMP(4) 一開始都比大家高？
error rate
## baseline
原始公式：
```c=
double compute_pi_baseline(size_t N)
{
double pi = 0.0;
double dt = 1.0 / N; // dt = (b-a)/N, b = 1, a = 0
for (size_t i = 0; i < N; i++) {
double x = (double) i / N; // x = ti = a+(b-a)*i/N = i/N
pi += dt / (1.0 + x * x); // integrate 1/(1+x^2), i = 0....N
}
return pi * 4.0;
}
```
## OpenMP
指定要用幾個thread來執行
- [參考資料](https://kheresy.wordpress.com/2006/09/22/%E7%B0%A1%E6%98%93%E7%9A%84%E7%A8%8B%E5%BC%8F%E5%B9%B3%E8%A1%8C%E5%8C%96%EF%BC%8Dopenmp%EF%BC%88%E4%BA%94%EF%BC%89-%E8%AE%8A%E6%95%B8%E7%9A%84%E5%B9%B3%E8%A1%8C%E5%8C%96/)
### reduction
支援的運算元有 +, , –, &, ^, |, &&, ||，它讓各個執行緒針對指定的變數擁有一份有起始值的複本（起始值是運算元而定，像 +, – 的話就是 0， 就是 1），平行化計算時，以各自複本做運算，等到最後再以指定的運算元，將各執行緒的複本合在一起。
## AVX SIMD
AVX 是一組 SIMD(Single Instruction Multiple Data)的指令集，擴展了原有的 MMX 及 Intel Streaming SIMD Extensions(Intel SSE)
支援 AVX 的硬體上會有 16 個 256-bit 的 YMM(YMM0-YMM15) 暫存器，還有一個 32-bit 的控制暫存器 MXCSR，舊有的 XMM 暫存器可以視為取 YMM 的 lower half，如下圖。
* _mm256d : 它並不是一種暫存器，是指可以用來載入到 AVX 暫存器的 “Data type”，double precision, 64bit
* _mm256_set1_pd(1.0)：將參數浮點數值放到 _mm256 變數的所有位置。
* _mm256_set_pd(dt * 3, dt * 2, dt * 1, 0.0)：將dt * 3, dt * 2, dt * 1, 0.0這些參數依序放入_mm256 變數，參數順序和放進去的次序相反。
```c=
double compute_pi_avx(size_t N)
{
double pi = 0.0;
double dt = 1.0 / N;
register __m256d ymm0, ymm1, ymm2, ymm3, ymm4;
ymm0 = _mm256_set1_pd(1.0);
ymm1 = _mm256_set1_pd(dt);
ymm2 = _mm256_set_pd(dt * 3, dt * 2, dt * 1, 0.0);
ymm4 = _mm256_setzero_pd(); // sum of pi
for (int i = 0; i <= N - 4; i += 4) {
ymm3 = _mm256_set1_pd(i * dt); // i*dt, i*dt, i*dt, i*dt
ymm3 = _mm256_add_pd(ymm3, ymm2); // x = i*dt+3*dt, i*dt+2*dt, i*dt+dt, i*dt+0.0
ymm3 = _mm256_mul_pd(ymm3, ymm3); // x^2 = (i*dt+3*dt)^2, (i*dt+2*dt)^2, ...
ymm3 = _mm256_add_pd(ymm0, ymm3); // 1+x^2 = 1+(i*dt+3*dt)^2, 1+(i*dt+2*dt)^2, ...
ymm3 = _mm256_div_pd(ymm1, ymm3); // dt/(1+x^2)
ymm4 = _mm256_add_pd(ymm4, ymm3); // pi += dt/(1+x^2)
}
double tmp[4] __attribute__((aligned(32)));
_mm256_store_pd(tmp, ymm4); // move packed float64 values to 256-bit aligned memory location
pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
return pi * 4.0;
}
```
## 時間處理與 time 函式使用
### clock()
### clock_gettime()
## Wall-clock time vs. CPU time
### Wall-clock time
顧名思義就是掛鐘時間，也就是現實世界中實際經過的時間，是由 kernel 裡的 xtime 來紀錄，系統每次啟動時會先從設備上的 [RTC](https://zh.wikipedia.org/wiki/%E5%AF%A6%E6%99%82%E6%99%82%E9%90%98) 上讀入 xtime。
```c
struct timespec xtime __attribute__ ((aligned (16)));
struct timespec {
__kernel_time_t tv_sec; /* seconds */
long tv_nsec; /* nanoseconds */
};
```
這個值是自 1970-01-01 起經歷的秒數，另外 Linux 核心也會決定每秒要發生幾次中斷 (透過常數 HZ) ，每次 timer interrupt，就會去更新 `xtime`。
另外注意，有時系統要與網絡中某個節點時間同步，那這個 Wall-clock time 可能會和 [NTP 伺服器](http://linux.vbird.org/linux_server/0440ntp.php)、[時區](https://zh.wikipedia.org/wiki/%E6%97%B6%E5%8C%BA)、[日光節約時間](https://zh.wikipedia.org/wiki/%E5%A4%8F%E6%97%B6%E5%88%B6)同步或使用者自己調整。所以「通常」不建議拿來量測程式時間，因為它不是一個穩定的時間，用專業一點的用語講，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....等等。
- **real** time : 也就是 Wall-clock time，當然若有其他程式同時在執行，必定會影響到。
- **user** time : 表示程式在 user mode 佔用所有的 CPU time 總和。
- **sys** time : 表示程式在 kernel mode 佔用所有的 CPU time 總和 ( 直接或間接的系統呼叫 )。
## 其他計算 pi 的方式
### Leibniz formula for π
考慮以下分解：
對兩邊從0到1去做積分
當時，除積分項以外的項收斂到萊布尼茨級數。同時，積分項收斂到0：
所以這便證明了萊布尼茨公式。
```c=
double compute_pi_leibniz(size_t N)
{
double pi = 0.0;
int tmp = 1;
for (size_t i =0; i < N; i++) {
pi += tmp / (2.0 * (double)i + 1.0);
tmp *= (-1);
}
pi *= 4;
return pi;
}
```
### Euler
通式
$\frac{π^2}{6}$ = $\sum_{i=1}^n \frac{1}{n^2}$
```c=
double compute_pi_euler(size_t N)
{
double pi = 0.0;
for (size_t i = 1; i <= N; i++) {
pi += (1 / pow(i, 2.0));
}
pi *= 6;
return sqrt(pi);
}
```
### Nilakantha 級數
π = 3 + 4/(2*3*4) - 4/(4*5*6) + 4/(6*7*8) - 4/(8*9*10) + 4/(10*11*12) - (4/(12*13*14) ...
通式大概長這樣
π = 3+$\sum_{i=2}^n \frac{4*(-1)^n}{(2n)(2n-1)(2n-2)}$
```c=
double compute_pi_nilakantha(size_t N)
{
double pi = 0.0;
int tmp = 1;
for(size_t i=2; i <= N; i++){
pi += tmp / ((2.0 * (double)i)*(2.0*(double)i-1)*(2.0*(double)i-2));
tmp *= (-1);
}
pi *= 4.0;
pi += 3;
return pi;
}
```
[某個資料](http://zh.wikihow.com/%E8%AE%A1%E7%AE%97%E5%9C%86%E5%91%A8%E7%8E%87-Pi)說這個比萊布尼茨快，結果好像沒有，可能是環境或是其他因素
## Reference
- [廖建富學長筆記](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I)
- [shelly4132](https://hackmd.io/s/HJrLU0dp)
- [萊布尼茲公式](https://zh.wikipedia.org/wiki/%CE%A0%E7%9A%84%E8%8E%B1%E5%B8%83%E5%B0%BC%E8%8C%A8%E5%85%AC%E5%BC%8F)
- [pi 計算方式](http://zh.wikihow.com/%E8%AE%A1%E7%AE%97%E5%9C%86%E5%91%A8%E7%8E%87-Pi)
- [SIMD 指令集](https://software.intel.com/sites/landingpage/IntrinsicsGuide/#techs=AVX&expand=4602)
Sign in with Wallet
Connect another wallet