# 2016q3 Homework1
## 作業要求
- 請依照各作業需求，自 GitHub 網站 fork 個別專案，並將連結貼於下方「作業區」
- phonebook, raytracing, compute-pi, clz
- 程式碼縮排有明確要求，務必遵守，這是團隊合作必要的準備，詳情見 phonebook 裡頭的 astyle 段落
- 分項作業請建立個別新的 HackMD 頁面，作為開發紀錄，並錄製「進度解說簡介」影片，在 20 分鐘之內如: 「開發紀錄(phonebook) / YouTube 連結」
- 共筆示範 / 錄影示範
- 錄製影片時，記得放大瀏覽器和終端機的字體，並確保成果在 1024x768 解析度得以清晰觀看，發話請接近麥克風，確保收音狀況良好。在 Ubuntu Linux 下，可善用 RecordMyDesktop 這套軟體來錄製桌面
- HackMD 教學和作業原則
- 請學員於 9 月 28 日中午前填寫課程「表單」，讓我們及早認識你。

## 圓周率之證明
參考: [wiki](https://zh.wikipedia.org/wiki/%E5%9C%93%E5%91%A8%E7%8E%87
)
這邊也有一個[寫的很詳細的網站](http://episte.math.ntu.edu.tw/articles/mm/mm_13_3_05/index.html)可以參考

一開始很難看得懂圓周率是怎麼得知的,雖然數學還了很多給數學老師了,但還是很好奇pi為什麼公式長這樣
為什麼一開始就能寫出
一直往下去查才找到為什麼現在就寫在這邊記錄一下這神奇的證明
### Leibniz's Proof
p.s.在寫這段的時候有查過可以用https的API來插入公式 , 有興趣請搜尋[MathJax](https://www.mathjax.org/)
OAT是一個單位圓的半圓,他的面積就是
現在我們需要考慮的是OPQT這塊地面積假設為 **C**,就是OTA扇形的面積減去OTA三角形面積的部分
再來我們考慮OPQ(由直線OP,直線OQ,弧PQ來建構而成),當PQ的位置越來越接近的時候PQ就會越趨近於一條直線設為**ds**
這樣的話我們就能將OPQ視為三角形來處理,再來延伸PQ的直線QR與通過圓心的直線OR垂直,然後我們就能求出三角形的面積(底*高/2)為
另外可以用[相似](https://tw.answers.yahoo.com/question/index?qid=20050116000010KK01312)的概念在PQ線段上還有一個小三角形與三角形OPQ相似所以得到
再來假設P點在水平線上的座標為 **x** , 然後要用點想像力 , 想像 **C** 就是無限個小三角行的組合所以就可以將 **C**寫成
在來使用[分部積分法](https://zh.wikipedia.org/wiki/%E5%88%86%E9%83%A8%E7%A9%8D%E5%88%86%E6%B3%95)計算 **C**的值就得到
再來請觀察一下整個圖形可以得到x,y的三角函數值
其中x用cos的2倍角公式可以替換出來
再來為了方便計算進行一連串的代換,讓x與y方便計算
參考:[Tangent is Sine divided by Cosine](https://proofwiki.org/wiki/Tangent_is_Sine_divided_by_Cosine)、[Secant is Reciprocal of Cosine](https://proofwiki.org/wiki/Secant_is_Reciprocal_of_Cosine)、 [Difference of Squares of Secant and Tangent](https://proofwiki.org/wiki/Difference_of_Squares_of_Secant_and_Tangent)
這樣就可以把 **x** 與 **y** 連接起來了
另外還能寫成[等比級數之和](https://proofwiki.org/wiki/Sum_of_Geometric_Progression)
最後就可以求出 **C** 的值
到這邊就快結束了 , 我們需要記得的是 **C** 代表著OPQT這塊面積要再加上三角形OTA才可以得到整個半圓的面積 , 最後可以證明出 , 真的佩服
最後觀察[泰勒展開式](https://proofwiki.org/wiki/Leibniz%27s_Formula_for_Pi)
當x值帶一時就是我們半圓的面積

## 第一次執行

N=400000000

time ./time_test_baseline
N = 400000000 , pi = 3.141593
5.28user 0.00system 0:05.28elapsed 100%CPU (0avgtext+0avgdata 1808maxresident)k
0inputs+0outputs (0major+86minor)pagefaults 0swaps
time ./time_test_openmp_2
N = 400000000 , pi = 3.141593
5.36user 0.00system 0:02.68elapsed 199%CPU (0avgtext+0avgdata 1816maxresident)k
0inputs+0outputs (0major+91minor)pagefaults 0swaps
time ./time_test_openmp_4
N = 400000000 , pi = 3.141593
5.52user 0.00system 0:01.54elapsed 357%CPU (0avgtext+0avgdata 1724maxresident)k
0inputs+0outputs (0major+92minor)pagefaults 0swaps
time ./time_test_avx
N = 400000000 , pi = 3.141593
1.64user 0.00system 0:01.64elapsed 100%CPU (0avgtext+0avgdata 1768maxresident)k
0inputs+0outputs (0major+85minor)pagefaults 0swaps
time ./time_test_avxunroll
N = 400000000 , pi = 3.141593
1.83user 0.00system 0:01.83elapsed 99%CPU (0avgtext+0avgdata 1692maxresident)k
0inputs+0outputs (0major+84minor)pagefaults 0swaps

## 求圓周率
### 黎曼積分

double computePi_v1(size_t N)
{
    double pi = 0.0;
    double dt = 1.0 / N; // dt = (b-a)/N, b = 1, a = 0
    for (size_t i = 0; i < N; i++) {
        double x = (double) i / N; // x = ti = a+(b-a)*i/N = i/N
        pi += dt / (1.0 + x * x); // integrate 1/(1+x^2), i = 0....N
    }
    return pi * 4.0;
}

### Leibniz Formula
將上述Leibniz formula轉換成c code

double compute_pi_lei(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);
    }
    pi *= 4; 
    return pi;
}

### Euler

### [Gauss–Legendre algorithm](https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm)

解釋部分[參考](http://drele.nxthosts.net/2016/04/parallel_calpi_gauss_legendre/)

偶間找到的【平行運算】計算高精度圓周率(PI)－使用迅速收斂演算法「Gauss–Legendre Algorithm」
「Gauss–Legendre Algorithm」是以迅速收斂著稱
在效率方面有大幅度的提升並且在寫法方面也簡潔很多。
–
以下先簡略介紹一下「Gauss–Legendre」演算法：
這是一個設計用來計算 PI 的演算法，它以迅速收斂著稱
在 27 次迭代過後就可以求得精確到小數點後 1 億位的 PI
不過缺點就是吃的記憶體稍大了點，計算到 1 億位吃到 1 GB 左右的記憶體。
–
求 PI 的步驟分為三項
首先要設定迭代項初始值：

接下來反覆執行以下的步驟，直到 a{0}、b{0} 之間的誤差達到所需精度：
\begin{align} a_{n+1} & = \frac{a_n + b_n}{2}, \\ b_{n+1} & = \sqrt{a_n b_n}, \\ t_{n+1} & = t_n - p_n(a_n - a_{n+1})^2, \\ p_{n+1} & = 2p_n. \end{align}
再利用以下算式求出 PI 近似值：
此演算法具有二階收斂性，本質上說就是演算法每執行一步正確位數就會加倍。
–

程式方面由於這樣迭代完後就是高精度了,所以已double的位數來說只需要做4次就夠了真是太方便啦！！！！

double compute_pi_Gauss()
{
    double a = 1.0;
    double b = 1.0/sqrt ( 2.0 );
    double t = 1.0/4.0;
    double x = 1.0;
    double y;
    int i;
    for ( i = 0; i < 4; i++ ) //without N , 4 in enought
    {
        y = a;
        a = (a + b)/2.0;
        b = sqrt ( b*y );
        t -= x*(y - a)*(y - a);
        x *= 2;
    }
	double pi = (a + b)*(a + b)/(4*t);
	return pi;
}