# 2017q1 Homework1(compute-pi)
contributed by < `baimao8437` >
### Reviewed by stanleytazi
+ 信賴區間的公式目前看到的都是 standard error 下去算
，而非standard deviation
[confidence interval](https://en.wikipedia.org/wiki/Confidence_interval)
+ 從比較圖來看好像都有看到OpenMP with 4 threads有比較大的抖動，有點好奇為什麼
+ default 設定跑25次來算，可以再試試看跑多次一點，看不同做法的時間會不會差的比較明顯
### Reviewed by `hunng`
- 利用多種方式計算並使用 AVX 實作
- 中間提到的不確定有沒有信賴區間，可以多次執行儲存結果後分析看看
- 而且有加上信賴區間比沒有加的抖動還要更大是不是出了什麼問題
### 花費時間 ~2/28~ ~-~ ~3/2~
* 讀書：9 hr
* 文件：4 hr
* coding: 3 hr
* ~~看著發呆：3 hr~~
## 目標設定
* 習慣 linux 系統基本操作
* 熟練 git 基本版本控管操作
* 繪圖工具 LibreOffice、gnuplot
* SIMD 實作練習
* 複習微積分？
* Hackmd 文件 LaTeX 數學式寫法
* 減少發呆時間...
## 前置作業
這個作業一開始我整個霧煞煞 沒了老溼的精美影片教學 頓失依靠
只好先讀神人們的共筆了解一些基本名詞與時間分析
## 開發環境
```
baimao@baimao-Aspire-V5-573G:~$ lscpu
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
型號： 69
Model name: Intel(R) Core(TM) i5-4200U CPU @ 1.60GHz
製程： 1
CPU MHz： 1711.242
CPU max MHz: 2600.0000
CPU min MHz: 800.0000
BogoMIPS: 4589.38
虛擬： VT-x
L1d 快取： 32K
L1i 快取： 32K
L2 快取： 256K
L3 快取： 3072K
NUMA node0 CPU(s): 0-3
```
## SIMD指令集
- AVX (Advanced Vector Extensions) 是 Intel 一套用來作 Single Instruction Multiple Data 的指令集。
- 我的電腦 CPU 為 intel i5-4200U
- 使用時需`#include <immintrin.h>`，並在 Makefile 中的編譯選項加入`-mavx`
## 效能分析
測試原版程式效能：
```clike
time ./time_test_baseline
N = 400000000 , pi = 3.141593
4.40user 0.00system 0:04.40elapsed 99%CPU (0avgtext+0avgdata 1808maxresident)k
0inputs+0outputs (0major+80minor)pagefaults 0swaps
time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
4.91user 0.00system 0:02.46elapsed 199%CPU (0avgtext+0avgdata 1772maxresident)k
0inputs+0outputs (0major+80minor)pagefaults 0swaps
time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
9.76user 0.00system 0:02.45elapsed 397%CPU (0avgtext+0avgdata 1792maxresident)k
0inputs+0outputs (0major+85minor)pagefaults 0swaps
time ./time_test_avx
N = 400000000 , pi = 3.141593
1.26user 0.00system 0:01.27elapsed 99%CPU (0avgtext+0avgdata 1816maxresident)k
0inputs+0outputs (0major+82minor)pagefaults 0swaps
time ./time_test_avxunroll
N = 400000000 , pi = 3.141593
1.22user 0.00system 0:01.22elapsed 99%CPU (0avgtext+0avgdata 1816maxresident)k
```
可以看到 user(user mode CPU time) 跟 system (kernel mode CPU time) 時間
也可以看到佔用了多少 CPU OpenMP 2 threads 最多佔用 200% CPU
OpenMP 4 threads 最多就佔用 400% CPU
由於我還沒做好算數學的心理準備 就先來畫畫好了~
想要建立出作業範例以 LibreOffice 建立的圖表（如下圖）：
先到 Makefile 的 gencsv 看到 seq 100 5000 25000
等於是 ```for(i=100; i<25000; i+=5000)```
我調成 ```seq 1000 200 400000``` 跑跑看將近兩千筆資料
然後把出來的 csv 以 LibreOffice Calc 開啟 選擇以逗號分隔
再把全部的資料選取好後 按下圖表 再做一些設定：
1. 圖表類型：設定為 線條圖 僅限線條 類型：平滑
2. 資料範圍：勾選 已欄表示的資料序列 + 第一欄作為標籤
3. 資料序列：把每條線名子改為對應 method
4. 圖表元素：x 軸為 N, y 軸為 Time(sec)
:::danger
這份報告之後可以拿來做研究所推甄、工作實習，甚至可作為論文的雛形，用與請精確 --jserv
:::
>Yes, Sir![name=baimao8437]
畫出來後:
預料之中 跟原圖差距甚遠
連一點指數上升的感覺都沒有 無法理解
而且原圖是 Baseline > OpenMP_2 > AVX > OpenMP_4 > AVX_unrolling
我的結果是 Baseline > OpenMP_4 > OpenMP_2 > AVX > AVX_unrolling
好像差了有點多 不過參考了別人的圖發現有人跟我一樣結果 就安心多了
補上 gnuplot perfomance + error rate:
## 信賴區間
可看出圖表顯示有一些極端值造成線條的抖動
因此參考[學長的筆記](https://hackmd.io/s/HySsF65p) 試圖建立95信賴區間來排除極端值
機統沒學好 信賴區間忘光光： 簡單來說就這三句話
```
計算平均值
計算標準差(standard deviation)
平均值的正負標準差乘以 2 就是 95% 信賴區間
```
標準差長這樣:
$$
s=\sqrt{{\frac{1}{n-1}\sum_{i=1}^n}(x_i-\bar x)^2}
$$
然後 makefile 中加入 ```-lm``` 一樣讓他跑 2000 筆資料
結果:
OpenMP 4 threads 還是抖了一下到讓我懷疑到底有沒有信賴區間...
補上 gnuplot:
可以看到 一樣的演算法 但是用不同的優化結果
AVX_unrolling 的優化會造成再某些 N 計算結果誤差比較大
其他的算法 error rate 皆相同
## Leibniz 公式證明
證明：
$$
\frac{\pi}{4}=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\cdots
$$
首先考慮
$$
1-x^2+x^4-x^6+x^8-\cdots=\frac{1}{1+x^2},\ |x|<1
$$
積分後
$$
x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\frac{x^9}{9}-\cdots=tan^{-1}x,\ |x|<1
$$
接著證明收斂
$$
\frac{1}{1+x^2}=1-x^2+x^4-x^6+x^8-\cdots+(-1)x^{2n}+\frac{(-1)^{n+1}x^{2n+2}}{{1+x^2}}
$$
兩邊積分後等於
$$
\frac{\pi}{4}=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\cdots+\frac{-1^n}{2n+1}+(-1)^{n+1}\int_0^1\frac{x^{2n+2}}{{1+x^2}}dx
$$
當$n\rightarrow\infty$時, 除積分項以外的項收斂到萊布尼茲級數, 同時積分項收斂至0
$$
\int_0^1\frac{x^{2n+2}}{{1+x^2}}dx<\int_0^1x^{2n+2}dx=\frac{1}{2n+3}\rightarrow0,\ when\ n\rightarrow\infty
$$
得證!
$$\pi=4\sum_{n=1}^\infty\frac{(-1)^n}{2n+1}$$
接著寫成 code:
```c=
double compute_pi_leibniz(size_t N)
{
double pi = 0.0;
int tmp = 1;
for (size_t i = 0; i < N; i++) {
pi += ((double)tmp / (2.0 * (double)i + 1.0));
tmp *= -1;
}
return pi * 4.0;
}
```
>第一次寫用了 pow(-1,i)去變換正負
>但太花費時間 所以用 tmp 去變換正負
然後 Leibniz 優化挑了 AVX 來實作:
```c=
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);//represent tmp
ymm1 = _mm256_set1_pd(1);
ymm2 = _mm256_set1_pd(2);
ymm4 = _mm256_setzero_pd(); // sum of pi
for (int i = 0; i <= N - 4; i += 4) {
ymm3 = _mm256_set_pd((double)i, (double)i + 1.0, (double)i + 2.0, (double)i + 3.0); // represent n
ymm3 = _mm256_mul_pd(ymm2, ymm3); // [2*n]
ymm3 = _mm256_add_pd(ymm1, ymm3); // 2*n[+1]
ymm3 = _mm256_div_pd(ymm0, ymm3); // [tmp/](2*n+1)
ymm4 = _mm256_add_pd(ymm4, ymm3); // pi += tmp/(2*n+1)
}
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;
}
```
原本很畏懼 AVX 但有了上面的範例就好懂多了 也比較容易寫出來
* gnuplot 圖表
有點緊密 不過可以看到 AVX_unrolling 還是比 Leibniz_avx 快一點點
但是 error rate 卻明顯 Leibniz_avx 優於 AVX_unrolling
從輸出 csv 來看這邊 因為 Leibniz & Leibniz_avx 的 error rate 又全部跟其他重疊了 所以只有一條線
## 其他算 $\pi$ 的方法 : Nilakantha method
$$
\pi=3+\frac{4}{2\times3\times4}-\frac{4}{4\times5\times6}+\frac{4}{6\times7\times8}-\frac{4}{8\times9\times10}+\cdots
$$
寫成 code 也不難 大概長這樣
```c=
double compute_pi_nilakantha(size_t N)
{
double pi = 0.0;
int tmp = 1;
double start = 0.0;
for (size_t i = 1; i < N; i++) {
start = (double)i * 2.0;
pi += (double)tmp / ((start) * (start + 1.0) * (start + 2.0));
tmp *= -1;
}
pi *= 4;
pi += 3;
return pi;
}
```
但是如果要用 AVX 優化的話 就要稍微細心一點了 改到頭昏腦脹
```c=
double compute_pi_nilakantha_avx(size_t N)
{
double pi = 0.0;
register __m256d ymm0, ymm1, ymm2, ymm3, ymm4, ymm5, ymm6;
ymm0 = _mm256_set_pd(1.0, -1.0, 1.0, -1.0);
ymm1 = _mm256_set1_pd(1);
ymm2 = _mm256_set1_pd(2);
ymm3 = _mm256_set1_pd(4);
ymm4 = _mm256_setzero_pd();//_mm256_set1_pd(4);
for (int i = 0; i < N - 4; i += 4) {
ymm6 = _mm256_set_pd((double)i + 1.0, (double)i + 2.0, (double)i + 3.0, (double)i + 4.0); //n
ymm6 = _mm256_mul_pd(ymm2, ymm6); //start = i*2
ymm5 = _mm256_add_pd(ymm6, ymm1); //(start+1)
ymm6 = _mm256_mul_pd(ymm6, ymm5); //start*(start+1)
ymm5 = _mm256_add_pd(ymm5, ymm1); //(start+2)
ymm6 = _mm256_mul_pd(ymm6, ymm5); //start*(start+1)*(start+2)
ymm6 = _mm256_div_pd(ymm3, ymm6); // 4/start*(start+1)*(start+2)
ymm6 = _mm256_mul_pd(ymm6, ymm0);
ymm4 = _mm256_add_pd(ymm4, ymm6);
}
double tmp[4] __attribute__((aligned(32)));
_mm256_store_pd(tmp, ymm4);
pi += tmp[0] + tmp[1] + tmp[2] + tmp[3];
pi += 3.0;
return pi;
}
```
我相信應該有更簡潔的寫法 但是至少這個結果是對的
來看圖表
效能來看最好的三個是 :
(時間少) AVX > AVX_unrolling > Leibniz_avx (時間多)
這個圖表看起來跟上面的 error 圖看起來一樣
因為 Nilakantha 的兩個算法 錯誤率都非常接近 0
錯誤率低的嚇人 我把 N 調小 單獨看 Nilakantha
```
N,error rate N,error rate
5,0.000608 55,0.000000
10,0.000079 60,0.000000
15,0.000023 65,0.000000
20,0.000010 70,0.000000
25,0.000005 75,0.000000
30,0.000003 80,0.000000
35,0.000002 85,0.000000
40,0.000001 90,0.000000
45,0.000001 95,0.000000
50,0.000001 100,0.000000
```
只要算到 50 幾項 就幾乎等於 $\pi$ 了
## 參考資料：
* [HankLo 的筆記](https://hackmd.io/s/HySsF65p)
* [HaoTse 的筆記](https://hackmd.io/s/ryHY0T_T)
* [廖建富學長筆記](https://hackpad.com/Hw1-Extcompute_pi-geXUeYjdv1I)
* [萊布尼茲公式](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)
* [Nilakantha 公式](https://scratch.mit.edu/projects/19839990/)
* [勃興大神的 makefile](https://github.com/0140454/compute-pi/blob/5300af6db6e952b044a68f0aaa9ae2a501d36bf2/Makefile)
Sign in with Wallet
Connect another wallet