Computing pi
===

###### tags: `sysprog-hw` `SIMD`

# 作業要求

* 在 GitHub 上 fork [compute-pi](https://github.com/sysprog21/compute-pi)，嘗試重現實驗，包含分析輸出
* 注意：要增加圓周率計算過程中迭代的數量，否則時間太短就看不出差異
* 詳閱廖健富提供的[詳盡分析](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I)，研究 error rate 的計算，以及為何思考我們該留意
* 用 [phonebok](/s/S1RVdgza) 提到的 gnuplot 作為 `$ make gencsv` 以外的輸出機制，預期執行 `$ make plot` 後，可透過 gnuplot 產生前述多種實做的效能分析比較圖表
* 可善用 [rayatracing](/s/B1W5AWza) 提到的 OpenMP, software pipelining, 以及 loop unrolling 一類的技巧來加速程式運作
* 詳閱前方 "Wall-clock time vs. CPU time"，思考如何精準計算時間
* 除了研究程式以外，請證明 [Leibniz formula for π](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80)，並且找出其他計算圓周率的方法
* 請詳閱[許元杰的共筆](https://embedded2015.hackpad.com/-Ext1-Week-1-5rm31U3BOmh)
* 將你的觀察、分析，以及各式效能改善過程，並善用 gnuplot 製圖，紀錄於「[作業區](https://hackmd.io/s/H1B7-hGp)」
----
# 實作

:::danger
如果在跑`make check`時發現下面的回應`Illegal instruction (core dumped)`代表電腦不支援AVX SIMD，7年前的電腦QQ
- Tip：
為甚麼 Linux 在發生 Illegal instruction 後仍然可以繼續執行？
CPU 會去捕捉 undefined instruction 然後回應上面的訊息，來做到 exception handler
[UD2](http://x86.renejeschke.de/html/file_module_x86_id_318.html)
:::

稍微修改一下 Makefile 就可以避免 AVX 的使用，執行`make check`跑出結果：

time ./time_test_baseline
N = 400000000 , pi = 3.141593
19.86user 0.01system 0:19.93elapsed 99%CPU (0avgtext+0avgdata 1640maxresident)k
0inputs+0outputs (0major+81minor)pagefaults 0swaps
time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
24.01user 0.11system 0:14.49elapsed 166%CPU (0avgtext+0avgdata 1648maxresident)k
0inputs+0outputs (0major+86minor)pagefaults 0swaps
time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
24.82user 0.17system 0:15.83elapsed 157%CPU (0avgtext+0avgdata 1628maxresident)k
0inputs+0outputs (0major+87minor)pagefaults 0swaps

Makefile 中的 topic gencsv 會決定圖形的取樣個數

gencsv: default
        for i in `seq 100 5000 1000000`; do \
                printf "%d " $$i;\
                ./benchmark_clock_gettime $$i; \
        done > result_clock_gettime.csv 

:::info
其中 seq 的用法為：
> seq [起始數字] [區間] [終止數字]
> eg:
>> seq 10 10 100
>> 數值將會是 10 20 30 ... 100
:::

gnuplot script

set xlabel 'N'
set ylabel 'Time(sec)'
set style fill solid
set title 'compute_pi time'
set term png enhanced font 'Verdana,10'
set output 'time.png'

plot "result_clock_gettime.csv" using 1:2 title 'baseline' with linespoints, \
     "result_clock_gettime.csv" using 1:3 title 'openmp_2' with linespoints, \
     "result_clock_gettime.csv" using 1:4 title 'openmp_4' with linespoints

第一次輸出的圖形(25倍時間)，可明顯看出資料分佈抖動的很厲害
> 註：我的電腦只有兩個threads，所以OpenMP_2 跟 OpenMP_4才會相同
1. 這邊有陷阱～
2. 要調整gettime()的時間精度
3. 用信賴區間(Confidence interval)消去極端值

// Baseline
    clock_gettime(CLOCK_ID, &start);
    for(i = 0; i < loop; i++) {
        compute_pi_baseline(N);
    }
    clock_gettime(CLOCK_ID, &end);
    printf("%lf ", (double) (end.tv_sec - start.tv_sec) +
           (end.tv_nsec - start.tv_nsec)/ONE_SEC);

> 這段 code 之中的時間其實是跑了25次後的時間長度！( loop = 25 )
----
加入信賴區間後
> 只取樣10次資料，抖動還是很大
----
> 取樣了100次資料
> 
----
> 取樣1000次資料
>
----
>將X軸取logscale
>
----
>Error with log scale (printf精度小數下6位)
----
>Error with log scale (printf精度小數下15位)
>- [ ] 為什麼是取15位？取幾位就夠了？
>
----
>Error with log scale (printf精度小數下15位)（seq 100 5000 1000000）
>
:::info
這邊的 avx_unrolling 的誤差是因為 code 之中

for (int i = 0; i <= N - 16; i += 16) {
        ymm14 = _mm256_set1_pd(i * dt);

        /*省略中間*/
		
        ymm6 = _mm256_add_pd(ymm6, ymm10);
        ymm7 = _mm256_add_pd(ymm7, ymm11);
        ymm8 = _mm256_add_pd(ymm8, ymm12);
        ymm9 = _mm256_add_pd(ymm9, ymm13);
    }

此處的 loop 執行一次等同於 baseline 版本執行16次
也因此我們的資料點（seq 100 5000 1000000）都沒有被16整除，導致 pi 的數值有少算的部分，才形成誤差，而 AVX 的版本是因為每個 for loop 只等同於 baseline 版本執行了4次，剛好可以整除資料點，才沒有算錯
ex：
N=20，上述迴圈只會執行一次，而導致資料少了20 - 16 = 4筆
:::

----
>Error with logscale（seq 103 5000 1000000)
>
:::info
可以看出當 N 不被4整除時，AVX就會產生誤差
:::
----
>經過修正後，成功將AVX的pi算回正確值！（seq 103 5000 1000000）
>
----
>修正AVX-unrolling（seq 103 5000 1000000）
>
----
:::success
已將修正後的程式碼推到github上了
:::

# 信賴區間

## 95%信賴區間

此處的運用，要了解自己的資料是透過取樣去預測母體或者對資料母體作運算，兩者的 95% 信賴區間的算法不同！

### 直接運算母體
- Lower Endpoint = avg(X) − 1.96 * σ
- Upper Endpoint = avg(X) + 1.96 * σ
- σ 為母體標準差
- n 為樣本資料數

### 從母體取樣本
- Lower Endpoint = avg(X) − 1.96 * σ/√n
- Upper Endpoint = avg(X) + 1.96 * σ/√n
- σ 為母體標準差
- n 為樣本資料數

## 標準差公式

----
# 參考資料
- [Leibniz formula for π](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80)
- [積分計算圓周率π](http://book.51cto.com/art/201506/480985.htm)
- [AVX SIMD](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I)
- http://sci.tuomastonteri.fi/programming/sse
- http://gruntthepeon.free.fr/ssemath/
- https://software.intel.com/sites/landingpage/IntrinsicsGuide/
- [Gnuplot](http://people.duke.edu/~hpgavin/gnuplot.html)
- [Confidence Intervals](https://www.maplesoft.com/support/help/Maple/view.aspx?path=MathApps/ConfidenceIntervals)
- [sqrt()](https://linux.die.net/man/3/sqrt)