# 2017q1 Homework1 (compute-pi)
contributed by <`Hong`>
###### tags: `sysprog2017` `compute-pi` `time` `gnuplot` `king1224`
### Reviewed by `kintaxxs`
- commit message 建議以一行72字元去做斷行
### Reviewed by `changyuanhua`
* 可以嘗試把 Leibniz 加上 AVX SIMD 和 Loop Unrolling 方法，這樣可以加快 real time 和 user time
* 可以回答一下延伸思考的問題
* 分析的圖要記得補喔
* 下面關於 gnuplot 製圖的內容有個小地方寫錯了，你的 makefile 中的 gencsv 應該是 100 ~ 100000，而不是 100 ~ 10000
* 在做圖部份，你可以試著把 makefile 中的 ` seq 100 1 100000 ` 中的 1 值改大一些，像是 100 ，這樣的話就不會因重疊覆蓋而看不到其他數據
## 開發環境
```
OS: 16.04.1 Ubuntu
Architecture: x86_64
CPU 作業模式： 32-bit, 64-bit
Byte Order: Little Endian
CPU(s): 4
On-line CPU(s) list: 0-3
每核心執行緒數： 2
每通訊端核心數： 2
Socket(s): 1
NUMA 節點： 1
供應商識別號： GenuineIntel
CPU 家族： 6
型號： 60
Model name: Intel(R) Core(TM) i5-4200H CPU @ 2.80GHz
製程： 3
CPU MHz： 2913.305
CPU max MHz: 3400.0000
CPU min MHz: 800.0000
BogoMIPS: 5587.41
虛擬： VT-x
L1d 快取： 32K
L1i 快取： 32K
L2 快取： 256K
L3 快取： 3072K
NUMA node0 CPU(s): 0-3
```
## 未優化版本
執行 `$ make check` 後可以得到各種實做方式計算 pi 所需的時間
* 執行 `$ make check`
```
time ./time_test_baseline
N = 400000000 , pi = 3.141593
3.38user 0.00system 0:03.38elapsed 99%CPU (0avgtext+0avgdata 1788maxresident)k
0inputs+0outputs (0major+82minor)pagefaults 0swaps
time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
3.41user 0.00system 0:01.72elapsed 198%CPU (0avgtext+0avgdata 1724maxresident)k
0inputs+0outputs (0major+86minor)pagefaults 0swaps
time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
6.66user 0.00system 0:01.72elapsed 387%CPU (0avgtext+0avgdata 1788maxresident)k
0inputs+0outputs (0major+89minor)pagefaults 0swaps
time ./time_test_avx
N = 400000000 , pi = 3.141593
0.98user 0.00system 0:00.98elapsed 99%CPU (0avgtext+0avgdata 1760maxresident)k
0inputs+0outputs (0major+83minor)pagefaults 0swaps
time ./time_test_avxunroll
N = 400000000 , pi = 3.141593
0.94user 0.00system 0:00.94elapsed 100%CPU (0avgtext+0avgdata 1700maxresident)k
0inputs+0outputs (0major+83minor)pagefaults 0swaps
```
分別執行測試程式後，各自記錄下來，並加上 make check 得到的 CPU 使用率。
* baseline
```
$ time ./time_test_baseline
N = 400000000 , pi = 3.141593
real 0m3.391s
user 0m3.388s
sys 0m0.000s
99%CPU
```
* openmp_2
```
$ time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
real 0m1.914s
user 0m3.800s
sys 0m0.000s
198%CPU
```
* openmp_4
```$ time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
real 0m1.727s
user 0m6.496s
sys 0m0.004s
387%CPU
```
* avx
```
$ time ./time_test_avx
N = 400000000 , pi = 3.141593
real 0m1.015s
user 0m1.000s
sys 0m0.012s
99%CPU
```
* avxunroll
```
$ time ./time_test_avxunroll
N = 400000000 , pi = 3.141593
real 0m0.952s
user 0m0.952s
sys 0m0.000s
100%CPU
```
經過 raytracing 的作業之後，發現適當的 loop unrolling 可以加快一些效能，在這份作業中的 baseline 版本是一個很大的 for-loop，即使沒辦法完全展開，我想試試看只展開成兩倍程式碼是否也能提昇效能：
```clike=
double compute_pi_baseline(size_t N)
{
double pi = 0.0;
double dt = 1.0 / N; // dt = (b-a)/N, b = 1, a = 0
double x;
//將原本迴圈中第 i 次及第 2i 次合併成一次做
size_t n = N/2; // 2i/N 可以看成 i/(N/2)，即為 i/n
size_t i = N/2;
do {
x = (double) i / N;
pi += dt / (1.0 + x * x );
x = (double) i / n;
pi += dt / (1.0 + x * x );
} while(--i);
return pi * 4.0;
}
```
不僅將迴圈中的程式碼展開兩倍，還將 for-loop 改成 do-while loop 希望能得到更好的 branch 效果，但最後執行時的效能還是與原來的幾乎一樣。
* * baseline_unrolling
```
N = 400000000 , pi = 3.425387
real 0m3.406s
user 0m3.404s
sys 0m0.000s
```
在一個很大的 for-loop 中，branch 與否幾乎是完全一模一樣的，在加上 compiler的優化，可以使 branch prediction 達到很高的準確度，展開前就沒有浪費太多時間在 branch 了，展開後自然看不到效能的提昇。
## gnuplot 製圖
我在學習 gnuplot 的時候發現了這個[教學網站](http://www.phy.fju.edu.tw/~ph116j/gnuplot_tutorial.pdf)，有一個我覺得很方便的地方，就是語法可以縮寫
* 原始版本
`plot 'out.txt' using 1:2 with linespoints pointtype 5`
* 縮寫版本
`plot 'out.txt' using 1:2 w lp pt 5`
如果你還想要對你的圖表做更多的設定，如顏色、寬度...等，那麼縮寫絕對可以讓你省下大量的程式碼，尤其是很多筆資料的時候。
#### 接下來來測試 compute-pi 的效能
執行 `$ make gencsv` 後可以得到一個結果檔案為 result_clock_gettime.csv，針對不同次數的 for 迴圈測試執行時間。
```
$ cat result_clock_gettime.csv
100,0.000022,0.000058,0.009070,0.000021,0.000010
5100,0.001121,0.000596,0.001139,0.000529,0.000458
10100,0.002167,0.001388,0.003336,0.001279,0.000908
15100,0.003467,0.001662,0.003165,0.001524,0.001017
20100,0.004343,0.002223,0.002196,0.001965,0.001176
```
可以看到第一行是迴圈次數，後面分別是五種方法的執行時間(順序如未優化版本中測試的五個方法)，特別的是這邊不是用空白，而是用逗號隔開每筆資料，這時候要使用 gnuplot 來畫圖的話，則以下列程式碼告知 gnuplot：
```
set datafile separator ","
```
我將 Makefile 中的 gencsv 改為 100~10000 每個 loop 次數都做一次，再使用 gnuplot 分別畫出各個版本效能
* 完整的 100~10000 次 loop 的時間圖，但幾乎都重疊覆蓋了
* 從其中抓取一小段，可以看到有不少抖動 ( jitter ) 發生
* 完整的 `runtime.gp`
```
reset
set datafile separator ","
set ylabel 'time(sec)'
set xlabel 'loop(times)'
set style fill solid
set key center top
set boxwidth 0.5
set title 'perfomance comparison'
set term png enhanced font 'Verdana,10'
set output 'runtime.png'
plot [:][0:0.3]'result_clock_gettime.csv' using 1:2 with linespoints title 'baseline', \
'' using 1:3 with linespoints title 'openmp2', \
'' using 1:4 with linespoints title 'openmp4', \
'' using 1:5 with linespoints title 'avx', \
'' using 1:6 with linespoints title 'avxunrolling'
```
## SIMD (Single instruction, multiple data)
這邊我想特別介紹 SIMD，之前修多處理機平行程式設計課程時就有學過 SIMD，如同字面上的意思，SIMD 就是多個相同的運算但是不同的資料輸入，例如有一支程式需要從 1+2, 3+4, ....一直加到 1000，像這樣只有一種加法運算，但每次需要運算的數值不一樣，就可歸類為 SIMD，如果程式中有這類運算，就很可能再經過方法優化，但之前上課時沒學到該怎麼優化。
前面使用到的方法中有一種為 [AVX](https://zh.wikipedia.org/wiki/SSE)，AVX 即為一種指令集，簡單來說，照著指令集的語法撰寫程式碼，即可同時處理 4 筆相同運算的指令，之後照著[這裡](http://book.51cto.com/art/201506/480985.htm)的範例簡單解說一下：
```clike=
double compute_pi_omp_avx(size_t dt){
double pi = 0.0;
double delta = 1.0 / dt;
__m256d ymm0, ymm1, ymm2, ymm3, ymm4;
ymm0 = _mm256_set1_pd(1.0);
ymm1 = _mm256_set1_pd(delta);
ymm2 = _mm256_set_pd(delta * 3, delta * 2, delta, 0.0);
ymm4 = _mm256_setzero_pd();
for (int i = 0; i <= dt-4; i += 4) {
ymm3 = _mm256_set1_pd(i * delta);
ymm3 = _mm256_add_pd(ymm3, ymm2);
ymm3 = _mm256_mul_pd(ymm3, ymm3);
ymm3 = _mm256_add_pd(ymm0, ymm3);
ymm3 = _mm256_div_pd(ymm1, ymm3);
ymm4 = _mm256_add_pd(ymm4, ymm3);
}
double tmp[4] __attribute__((aligned(32)));
_mm256_store_pd(tmp, ymm4);
pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
return pi * 4.0;
}
```
* line 4：
* __m256d 是一種資料型態，用來宣告 AVX 所需的暫存器
* line 5:
* _mm256_set1_pd(1.0) 將暫存器中四筆資料都設成 1.0
* line 7:
* _mm256_set_pd 與上一個不同的是，這裡必須傳入四個參數，分別代表四筆資料的值
* line 8:
* _mm256_setzero_pd 將暫存器四筆資料都設成 0
* line 12~16:
* _mm256_add_pd, _mm256_mul_pd, _mm256_div_pd 分別就是將兩個傳入的參數中各自的四筆資料，一一對應做相加、相乘及相除等運算
* line 18:
* __attribute__可以用來設置一些屬性，而 aligned(32) 表示這個屬性設置為以 32 為單位對齊，[更多範例](http://huenlil.pixnet.net/blog/post/26078382-%5B%E8%BD%89%5Dgnu-c-__attribute__-%E6%A9%9F%E5%88%B6%E7%B0%A1%E4%BB%8B)
* line 19:
* _mm256_store_pd 則可以將 ymm 暫存器中的四筆資料，存回我們常用的 array。
## 證明 [Leibniz formula for π](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80)
先讓我們看一下[泰勒級數](https://zh.wikipedia.org/zh-tw/%E6%B3%B0%E5%8B%92%E7%BA%A7%E6%95%B0)：
arctan(x) = $\sum_{n=0}^∞ \frac{(-1)^nx^{2n+1}}{2n+1}$
如果將 x=1 代入是不是就跟萊布尼茲算 π 的公式一樣了！
因此我這邊證明 arctan(x) 的泰勒級數：
* 先將 arctan(x) 微分後可以得到 $\frac{1}{1+x^2}$
* 這是一個首項為 1，等比為 $-x^2$ 的無窮等比級數
* 即為 1 - x^2 + x^4 - ..... ＝ $\sum_{n=0}^∞ (-1)^{n+1}x^{2n}$
* 再積分回去即可得原來的 arctan(x) ＝ C + x - $\frac{x^3}{3}$ + $\frac{x^5}{5}$ - ..... = C + $\sum_{n=0}^∞ \frac{(-1)^nx^{2n+1}}{2n+1}$
* 先將 x=0 代入可得常數 C=0
* 得證 arctan(x) = $\sum_{n=0}^∞ \frac{(-1)^nx^{2n+1}}{2n+1}$
再將 x=1 代入，得 $\frac{π}{4}$ = $\sum_{n=0}^∞ \frac{(-1)^n}{2n+1}$
## 其他計算 π 的方法 (待補圖)
### Monte Carlo method
[Monte Carlo method](https://openhome.cc/Gossip/AlgorithmGossip/MathPI.htm) 主要是用於解與面積有關的題目，而 π 就等於一個單位圓的面積，做法為先畫出一個比所求面積更大的矩形，並在這矩形中亂數取點，並紀錄總共取的點數量與落在所求面積內的點數量，則可以得到總面積與所求面積的比例。
* 程式碼
```clike=
double compute_pi_Monte_Carlo(size_t N)
{
srand(time(NULL));
size_t IN = 0;
size_t n=N;
double x;
double y;
do {
x = (double)rand() / RAND_MAX;
y = (double)rand() / RAND_MAX;
if( x*x + y*y <=1) ++IN;
} while(--n);
return (double)IN/N*4.0;
}
```
* 執行結果
```
$ time ./time_test_Monte_Carlo
N = 400000000 , pi = 3.141577
real 0m8.357s
user 0m8.344s
sys 0m0.000s
```
這是一個可以廣泛應用在很多算面積問題的方法，實作上也不算困難，因此算是一個非常方便的演算法，但隨機總會有誤差值，即使把總點數量取的非常大，也還是會有些許誤差，因此此方法的準確性為本篇所有方法中最低，而大量的樣本數也會造成大量的時間花費。
### Leibniz formula
經過上一段的證明後可得 Leibniz formula for π，可以看到公式中每一項間是正負交錯，我試想能否將兩項合併以省略正負交錯的判斷，觀察後得一規律：
$1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9}$ ...
= $\frac{2}{1 \times 3} +\frac{2}{5 \times 7} + \frac{2}{9 \times 11} +$ ...
= $\frac{2}{1 \times 3} + \frac{2}{1 \times 3 + 32} + \frac{2}{5 \times 7 + 64} +$ ...
分母為 3+0, 3+32, 3+64, 3+96 ......
即為首項為 0、公差為 32 的等差級數和再加上 3 的值。
在程式碼中若以 $3+32\times0 , 3+32\times1 , 3+32\times2$... 來運算則會需要很多的乘法運算，因此我利用 tmp2 變數來暫存 0, 32, 64... 這些值，省掉乘法運算。
* 程式碼
```clike=
double compute_pi_Leibniz(size_t N)
{
size_t n = N/2;
size_t tmp = 3;
size_t tmp2 = 0;
double pi=0;
do {
tmp += tmp2;
tmp2 += 32;
pi += (double)1/ (double)tmp;
} while(--n);
return pi*8.0;
}
```
* 執行結果
```
$ time ./time_test_Leibniz
N = 400000000 , pi = 3.141593
real 0m0.843s
user 0m0.840s
sys 0m0.000s
```
經過優化後的 Leibniz formula 執行效果則變成本頁所有方法中最佳的一種
>如果再將兩項兩項合併後，會得到更有效率的算法嗎？
>如果可以再合併，那是否可以無限合併下去？
>[name=洪正皇][color=#50ddff]
>可以再合併一次，可將兩個除法減少為一個乘法
>無限合併待思考
>[name=洪正皇][color=#c0dd68]
## 如何精準計算時間
關於 time 有很多的分類，如 user time, system time, Wall-clock time 及 CPU time ...等，因應不同的需求，找到合適的方式計時，這邊針對要挑選 Wall-clock time 或是 CPU time 來測量整體程式時間做分析：
* Wall-clock time：真實經過的時間，而 Linux 的時鐘中有兩個全域變數
* xtime : 紀錄著自 1970/1/1 00:00:00 以來經過的秒數，但會隨著使用者調整時區，或是自動從網路上更新時間而被影響，因此不能保證永遠都是遞增的，即用此紀錄時間非常危險
* jiffies : 紀錄著從電腦開機到現在總共的時鐘中斷次數，永遠為遞增。
* CPU time：總共使用 CPU 的時間
## Error Rate
[廖健富的共筆](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I)中使用了 error rate來分析精準度：
```clike=
// Error rate calculation
#define M_PI acos(-1.0)
double pi = compute_pi(n);
double diff = pi - M_PI > 0 ? pi - M_PI : M_PI - pi;
double error = diff / M_PI;
```
不論是用何種離散的方式來計算 π，永遠不可能做到項次有無限多項，那麼如何才算精準呢？透過計算 error rate 的方式，設定可以容忍的誤差值，即可找出最少需要多少項次才能達到需要的精確度，就可以不用每次都取很大的 N，做很大量又耗時的運算。
## 參考資料：
[TempoJiJi 共筆](https://hackmd.io/s/Hy7tssw6#)
[jkrvivian 共筆](https://hackmd.io/s/rkGc0sOp#)
[泰勒展開](http://calculus.yuyumagic424.net/wp-content/uploads/2012/12/%E6%B3%B0%E5%8B%92%E5%B1%95%E9%96%8B.pdf)
[gnuplot 教學](http://www.phy.fju.edu.tw/~ph116j/gnuplot_tutorial.pdf)
[時間處理與 time 函式使用](https://charles620016.hackpad.com/ep/pad/static/OfmE372fliX)
Sign in with Wallet
Connect another wallet