# 2017 q1 Homework1 (compute-pi)
contributed by <`zmke`>

### Reviewed by <`rayleigh0407`>
* 建議可以證明一下 Leibniz formula

### Reviewed by <`twzjwang`>
- 執行時同時使用其他程式，為什麼 omp_4 受到的影響最大?
- 可以探討作業要求中時間函式部分提到的問題、延伸思考

## 開發環境
OS: Ubuntu 16.04 LTS
Architecture: x86_64
CPU 作業模式： 32-bit, 64-bit
Byte Order: Little Endian
CPU(s): 4
Model name: Intel(R) Core(TM) i5-4200H CPU @ 2.80GHz
L1d 快取： 32K
L1i 快取： 32K
L2 快取： 256K
L3 快取： 3072K

## AVX
* Advanced Vector Extensions
* Intel 開發的 Single Instruction Multiple Data 指令集，在一個 cycle 內可以同時對多個數據執行同樣的指令
* ymm0~ymm15 共 16 個暫存器，每個暫存器內可以放256 bits
* e.g. 8個 32-bit float; 4個 64-bit double

## 時間處理與 time 函式使用
* NTP (Network Time Protocol)
* 透過網路封包交換，將電腦之間的時鐘校正和同步
* clock()
* return CPU time
* 回傳值是 `clock_t` type ，要除以 `CLOCKS_PER_SEC` 才能取得以秒為單位的時間
* On a 32-bit system where CLOCKS_PER_SEC equals 1000000 this function will return the same value approximately every 72 minutes. （ 超過72分鐘會 wrap around ）
* 在 muti threads 的情況下， CPU time 是每條 thread 的加總
* clock_gettime()
* `int clock_gettime(clockid_t clk_id, struct timespec *tp);`
* retrieve and set 的類型由 clk_id 決定
* CLOCK_REALTIME : system-wide clok ，會回傳 real time 也就是 wall-clock ，如果使用者任意更動系統時間，會影響測量結果
* CLOCK_MONOTONIC : 不受 discontinuous jumps 影響
* 精度為 nanoseconds
* time()
* wall-clock
* return 1970-01-01 00:00:00 +0000 (UTC) 到目前為止的時間，精度為 seconds

以前寫作業的時候曾經遇過 multi thread 時間測量的問題，一開始是用 clock() ，發現測出來的時間跟預期相差很大，現在終於知道原因了

## 原始版本
* 修改 Makefile

gencsv: default
	for i in `seq 500 500 250000`; do \
		printf "%d," $$i;\
		./benchmark_clock_gettime $$i; \
	done > result_clock_gettime.csv	
    
plot: gencsv
        gnuplot script.gp

從 N=500 開始，間隔 500 測試到 N=250000
`$ make plot`
起伏有點誇張，尤其是 omp_4，發現可能的原因是我邊等輸出邊用 youtube 聽音樂，把網頁關起來再試一次
這次在等輸出的時候沒有對電腦進行任何操作，明顯比上一次穩定很多

## Leibniz

double compute_pi_leibniz(size_t N){
	double pi = 0.0;
	for(size_t i = 0; i < N; i++)
		pi += pow(-1.0, i) / (2*i + 1);
	return pi * 4.0;
}

omp 版本

double compute_pi_leibniz_openmp(size_t N, int threads){
	double pi = 0.0;
	#pragma omp parallel num_threads(threads)
    {
        #pragma omp for reduction(+:pi)
        for(size_t i = 0; i < N; i++)
			pi += pow(-1.0, i) / (2*i + 1);
    }
    return pi * 4.0;
}

`$ make plot`
發現 Leibniz 求圓周率明顯比其他方法慢，猜測是 `pow(-1.0, i)` 太浪費時間，修改版本如下

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

omp

double compute_pi_leibniz_openmp(size_t N, int threads){
	double pi = 0.0;
	#pragma omp parallel num_threads(threads)
    {
        #pragma omp for reduction(+:pi)
        for(size_t i = 0; i < N; i++)
			pi += (i % 2) ? (-1.0 / (2*i + 1)) : (1.0 / (2*i + 1));
    }
    return pi * 4.0;
}

模仿原本的 avx，練習把 Leibniz 改成用 AVX 計算

double compute_pi_leibniz_avx(size_t N){
	double pi = 0.0;
	register __m256d ymm0, ymm1, ymm2, ymm3, ymm4;
	ymm0 = _mm256_set_pd(1.0, -1.0, 1.0, -1.0);		//pow(-1.0, i)
	ymm1 = _mm256_set1_pd(2.0);
	ymm2 = _mm256_set1_pd(1.0);
	ymm4 = _mm256_setzero_pd();						//sum of pi

	for(int i = 0; i <= N-4; i += 4){
		ymm3 = _mm256_set_pd(i, i+1.0, i+2.0, i+3.0);
		ymm3 = _mm256_mul_pd(ymm3, ymm1);			//2*i
		ymm3 = _mm256_add_pd(ymm3, ymm2);			//2*i+1
		ymm3 = _mm256_div_pd(ymm0, ymm3);			//pow(-1.0, i) / (2*i+1)
		ymm4 = _mm256_add_pd(ymm4, ymm3);			//pi += pow(-1.0, i) / (2*i+1)
	}
	double tmp[4] __attribute__((aligned(32)));
	_mm256_store_pd(tmp, ymm4);
	pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
	return pi * 4.0;
}

`$ make plot`
發現把 `pow(-1.0, i)` 改掉之後， leibinz_omp4 變成執行速度最快的版本

## Compute Error
以 arccos(-1) 為基準計算誤差

    // Baseline
	result = compute_pi_baseline(N);
	error = (result - PI > 0) ? result - PI : PI - result;
	printf("%lf,", error);

	// Leibniz
	result = compute_pi_leibniz(N);
	error = (result - PI > 0) ? result - PI : PI - result;
	printf("%lf\n", error);

`$ make plot_error`
兩條線幾乎重疊，表示在同樣的 N 精度差距不大

把判斷正負的判斷式移除後再畫一次，發現兩個求圓周率方法分別是從上下逼近
`$ make plot_error`

## Reference
[AVX是什麼？AVX指令集技術與應用解析](http://www.expreview.com/13236.html)
[TempoJiJi 筆記](https://hackmd.io/s/Hy7tssw6)
[clock()](https://linux.die.net/man/3/clock)
[clock_gettime()](https://linux.die.net/man/2/clock_gettime)
[NTP](http://linux.vbird.org/linux_server/0440ntp.php)
[Leibniz formula for pi](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80)