# 2017q1 Homework01 (compute-pi)
contrubuted by <`hugikun999`>

### Reviewed by `heathcliffYang`
- 將不同的計算 pi 的方法實作很棒，希望可以引入更多方法的比較[Gauss - Legendre method](https://forums.anandtech.com/threads/calculating-pi-gauss-legendre-method.1519452/)、還有多執行緒版[Hyper PI](http://virgilioborges.com.br/hyperpi/)

## inplement
* cycle finding
對於一個f(x)，給定一初始值，請找到第一次開使重覆的地方。（即Xi = Xj）一旦出現這種情況就會週期性出現。

memset(void *str,int c,size_t z)

說明：將這個指標所指向的前z個字節設為"c";

回傳值：pointer。

* Memoization
一種最佳化的技巧，紀錄function call的結果，當有同樣的輸入發生時直接回傳cache的東西。

* [相關pi的演算法](https://en.wikipedia.org/wiki/Category:Pi_algorithms)

* 時間函數
（1）UTC：世界標準時間。
（2）calender time:日曆時間，即相對時間，對不同的編譯系統會有不同標準時間點。
（3）clock tick:計時單元，其時間長短由CPU決定。

* 時間精準度：
（1）clock_gettime():取得目前時間。（1ns）
（2）Windows:QueryPerformanceCounter():取得rdtsc的counter數。( < 1ns)
(3)clock():精確度只到ms。

* xtime：
（1）microsecond:微秒
（2）ms:毫秒
(3) 從1970/1/1開始計時，所謂的實時時鐘，精確度是us。

* jiffies:
(1)每一個tick是以系統頻率決定，若是200MHz則可以精確到5ns。
## AVX
* 延伸SSE指令集
[這裡](https://isite.tw/2016/01/24/14639)有歷代指令集的介紹。最大的特點是將XMM 128bit -> YMM 256bit，提高CPU在處理特定向量和浮點數運算的效率。
* 4個operands
通過VEX prefix實現4個operands，表示為dest、src1、src2、src3。雖然有4個operand，但速度不會到4x。VEX又可分為128和256版，分別集成escape prefix、SIMD prefix、REX prefix。
* 指令編碼格式
VEX的部份
（1）提供寄存器寻址。
（2）提供 escape prefix 与 SIMD prefix
（3）指示 operand 是 128 位还是 256 位
（4）提供 REX.RXB 的作用

* 實現不改變destination
* AVX512
再將256bit擴大到512bit。
* AVX相關語法
[這裡](https://software.intel.com/en-us/articles/introduction-to-intel-advanced-vector-extensions)是基本介紹每個字的意思。
[這裡](https://software.intel.com/sites/landingpage/IntrinsicsGuide/#expand=4635,198,100,112)是AVX還有其它SSE的版本一些相關用法，裡面有每個function的做法和對應的assembly code。
## compute_pi
在開使之前先用```time```測試了raytracing，結果發現由main.c的clock_gettime()會小於由```time```所測出來的real。

# Rendering scene
Done!
Execution time of raytracing() : 1.059423 sec

real	0m1.064s
user	0m8.364s
sys	0m0.028s

### Makefile
(一)
我在寫Makefile的時候發現我的tab後面加的東西沒有變成指令，查了一下覺得可能是因為我好像把tab變成四個空白。
[參考資料](http://koen.iteye.com/blog/602240)
（二）
在```Makefile```中的```gencsv```裡有寫到100、500、25000，經過測試知道100是N的初始值，500是每次增加的大小，25000是N的上限。

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

### gnuplot
這次使用的是.csv檔和上一次不一樣，按照上次的去做runtime.gp會一直出現```"runtime.gp", line 11: warning: Skipping data file with no valid points```這個錯勿，參考[Vivian Lin](https://hackmd.io/MYZgpgTAJmCsUFpYDNiwQFnMhBOAhmAGx5gYCMY+IUuRy5QA?view#)在其中多加一行```set datafile separator ","```就可以成功吃到資料了。

set xlabel 'N'
set ylabel 'Time(sec)'
set autoscale fix
set style fill solid
set datafile separator ","
set title 'runtime'
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, \
     'result_clock_gettime.csv' using 1:5 title 'AVX' with linespoints, \
     'result_clock_gettime.csv' using 1:6 title 'AVX+unroll' with linespoints, \

### file

* openMP
在編譯時記得要在gcc後加上參數```-fopenmp```，檔案內也要記得```#include <omp.h>```，不然就算加上一些openmp的語法也不會有作用。[這裡](https://kheresy.wordpress.com/2006/09/13/%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%8C%EF%BC%89%E8%AA%9E%E6%B3%95%E8%AA%AA%E6%98%8E/)有語法的介紹和一些範例。

#pragma omp parallel  num_threads(threads)

這個意思是指定threads，在（）內填入想指定的數目，通常thread越高就會跑越快，但是會有一個極限。

#pragma omp for private(x) reduction(+:pi)

for是將底下的for迴圈拆開來分給每個thread去執行。
private(x)是讓每個thread都有一份自己的x而不去影響到其它thread所存的東西，如果改到別的thread存的東西會造成錯誤。
reduction(+:pi)很類似private()，讓大家都有一份複本各自運算，不同的是最後再將所有的```pi```相加起來。

* AVX

    register __m256d ymm0, ymm1, ymm2, ymm3, ymm4;

宣告出4個256bit的register。```__m256d```是指4個64bit的double。

ymm0 = _mm256_set1_pd(1.0);

將4個double分別紀錄1.0進去。

ymm2 = _mm256_set_pd(dt * 3, dt * 2, dt * 1, 0.0);

將4個double分別照順序記錄dt*3、dt*2 、dt*1、dt。

ymm4 = _mm256_setzero_pd();

全部bit記錄為0。

ymm3 = _mm256_add_pd(ymm3, ymm2);

將ymm3加上ymm2存到ymm3中。add可以替換sub、mul、div之類的算術方法，也可換為store之類的功能。

double tmp[4] __attribute__((aligned(32)));

(1)__attribute__可以設置函數屬性（Function Attribute）、變量屬性（Variable Attribute）和類型屬性（Type Attribute）。
(2)
aligned(32)
該聲明將強制編譯器確保（盡它所能)分配空間時採用32字節對齊方式。取32byte因為一個double是8byte，須要4個即為32byte。

### 未優化
(一)
(1)```make check```

time ./time_test_baseline
N = 400000000 , pi = 3.141593
3.49user 0.00system 0:03.49elapsed 100%CPU (0avgtext+0avgdata 1780maxresident)k
0inputs+0outputs (0major+85minor)pagefaults 0swaps
time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
3.58user 0.00system 0:01.79elapsed 199%CPU (0avgtext+0avgdata 1788maxresident)k
0inputs+0outputs (0major+86minor)pagefaults 0swaps
time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
3.80user 0.00system 0:00.95elapsed 398%CPU (0avgtext+0avgdata 1784maxresident)k
0inputs+0outputs (0major+92minor)pagefaults 0swaps
time ./time_test_avx
N = 400000000 , pi = 3.141593
1.01user 0.00system 0:01.01elapsed 100%CPU (0avgtext+0avgdata 1756maxresident)k
0inputs+0outputs (0major+83minor)pagefaults 0swaps
time ./time_test_avxunroll
N = 400000000 , pi = 3.141593
0.98user 0.00system 0:00.98elapsed 100%CPU (0avgtext+0avgdata 1792maxresident)k
0inputs+0outputs (0major+84minor)pagefaults 0swaps

（2）由100～50000的圖


（二）
我更改了```Makefile```裡面的N，發現了N從0～20的區間都會跑很久，甚至是懷疑N=1就卡住了，從perf top 監看，發現一直```compute_pi_avx```處於熱點，目前還不確定是什麼原因，等看完AVX相關的文件應該就會比較了解。

## 優化
### 信賴區間
(一)
從圖中可以看到有些超出範圍太多，可以視為例外並沒有太多參考價值，故對這些數據取2個標準差可以使其有更精準的量測結果。
我在```benchmark_clock_gettime.c```中多寫了一個用來計算標準差的function，當中有多用到```math.h```這個library，所以在```Makefile```中必須多加上```-lm```才能執行。
（1）Makefile

default: computepi.o
    $(CC) $(CFLAGS) computepi.o time_test.c -DBASELINE -o time_test_baseline
    $(CC) $(CFLAGS) computepi.o time_test.c -DOPENMP_2 -o time_test_openmp_2
    $(CC) $(CFLAGS) computepi.o time_test.c -DOPENMP_4 -o time_test_openmp_4
    $(CC) $(CFLAGS) computepi.o time_test.c -DAVX -o time_test_avx
    $(CC) $(CFLAGS) computepi.o time_test.c -DAVXUNROLL -o time_test_avxunroll
    $(CC) $(CFLAGS) computepi.o benchmark_clock_gettime.c -o benchmark_clock_gettime -lm

注意```-lm```必須加在gcc的最後面，如果直接加在```CFLAGS```後面會錯。
（2）圖（sample = 10）

雖然次數的間隔比之前大，但是從圖中可以看到趨勢線已經比較少有極端值的存在了。

(3)圖（sample = 25）
這次取25個sample震盪幅度又變更小了，可是我的openmp4很奇怪，會有一些超乎常理的暴增。目前覺得是因為我的電腦只有2個threads，所以宣告成四個threads時會有互搶的問題，所以才有些會超高。另一個是有些openmp4比AVX還要低，AVX一次處理4個register，openmp4也像是一次跑4個的概念，但是卻比AVX還要快。

(4)圖 （sample = 25）

驗證在（三）中的想法，我將AVX中加上了openmp去嘗試，結果發現圖非常的醜，應該是不能用private去讓每個thread有自己的256bit register，可是也有一些低的AVX+openmp可以低到跟openmp4一樣了，可見最低值還是有下降。
另外我增跑了openmp8，雖然沒有比openmp4更低卻有比openmp4更穩定。

### Error check
（一）
除了花費時間外，從下面的圖中可以看出來由五種方法所算出來的正確性，剛開始N很小時所做出來的正確性很低，當N變大時正確性也提高。

（二）
AVX和AVX+unroll會有震盪是因為每次加4和每次加16，如果N不是4或著16的倍數就會造成少算的結果，所以會產生震盪。將N和每次的差都改為16的倍數就不會和其它演算法有差。


### [Machin-like formula](http://personalpages.to.infn.it/~zaninett/pdf/machin.pdf)
$4/\pi$ = $4*arctan1/5$ - $arctan1/239$

實作上沒有很困難的地方。
* arctan 實作

double arctan(double x,int N)
{
    double answer = 0.0;
    double tmp = 0.0;
    for (int i = 0; i < N; i++) {
        tmp = pow(-1,(double) i) * pow(x, (double) (2*i + 1)) / (double) (2*i +1);
        answer += tmp;
        printf("%e\n", tmp);
    }
        printf("answer: %e\n", answer);

    return answer;
}

* Machin-like 實作

double Machin_like(size_t N)
{
    double pi = 0.0;

    pi = 4 * arctan(0.2, N) - arctan(1/239, N);

    return pi * 4.0;
}

由結果可以得知Machin-like的收斂非常快，計算量也比其他的幾種方法更高，在同樣的 N 下所花費的時間式其他的好幾十倍。下面列出$archtan 1/5$ 和 $archtan 1/239$在前10次的運算的值，代表只要 `N = 10` 就可以計算到10的負15次方，甚至是10的負47次方的大小，對於其他的演算法可以說收斂速度非常快。
59265358979222781954376841
* arctan1/5 前10次運算值

2.000000e-01
-2.666667e-03
6.400000e-05
-1.828571e-06
5.688889e-08
-1.861818e-09
6.301538e-11
-2.184533e-12
7.710118e-14
-2.759411e-15

* arctan1/239 前10次運算值

4.184100e-03
-2.441659e-08
2.564723e-13
-3.207131e-18
4.366932e-23
-6.255045e-28
9.265822e-33
-1.405854e-37
2.171635e-42
-3.401624e-47

Machin-like 在 `N = 5` 之前就已經將值逼近到接近3.14了，相對於其他方法需要切割上千份才能夠有如此的精準度，但是從圖中看來逼近的值並不是正確的，這部份還需要修正。

上圖不會逼近正確的值，原因在於 $1/239$ 太小，傳入 `pow` 會造成運算上的不精確，傳值給 `pow` 時必須強轉成 double 才可以得到精確的值。

雖然說精確度比起其他更高，但是花費的時間也是其他所不能比較的。

Machin-like 的誤差幾乎看不出來的，而且是從 N 極小時就有如此的精準度。