2016q3 Homework1

作業要求

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

圓周率之證明

參考: wiki
這邊也有一個寫的很詳細的網站可以參考

一開始很難看得懂圓周率是怎麼得知的,雖然數學還了很多給數學老師了,但還是很好奇pi為什麼公式長這樣

為什麼一開始就能寫出

一直往下去查才找到為什麼現在就寫在這邊記錄一下這神奇的證明

Leibniz's Proof

p.s.在寫這段的時候有查過可以用https的API來插入公式 , 有興趣請搜尋MathJax

OAT是一個單位圓的半圓,他的面積就是
現在我們需要考慮的是OPQT這塊地面積假設為 C,就是OTA扇形的面積減去OTA三角形面積的部分
再來我們考慮OPQ(由直線OP,直線OQ,弧PQ來建構而成),當PQ的位置越來越接近的時候PQ就會越趨近於一條直線設為ds
這樣的話我們就能將OPQ視為三角形來處理,再來延伸PQ的直線QR與通過圓心的直線OR垂直,然後我們就能求出三角形的面積(底*高/2)為

另外可以用相似的概念在PQ線段上還有一個小三角形與三角形OPQ相似所以得到

再來假設P點在水平線上的座標為 x , 然後要用點想像力 , 想像 C 就是無限個小三角行的組合所以就可以將 C寫成

在來使用分部積分法計算 C的值就得到

再來請觀察一下整個圖形可以得到x,y的三角函數值
其中x用cos的2倍角公式可以替換出來
再來為了方便計算進行一連串的代換,讓x與y方便計算

參考:Tangent is Sine divided by Cosine、Secant is Reciprocal of Cosine、 Difference of Squares of Secant and Tangent
這樣就可以把 x 與 y 連接起來了

另外還能寫成等比級數之和

最後就可以求出 C 的值

到這邊就快結束了 , 我們需要記得的是 C 代表著OPQT這塊面積要再加上三角形OTA才可以得到整個半圓的面積 , 最後可以證明出 , 真的佩服

最後觀察泰勒展開式

當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

解釋部分參考

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

接下來反覆執行以下的步驟，直到 a{0}、b{0} 之間的誤差達到所需精度：

$$\begin{array}{rl}{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{array}$$
再利用以下算式求出 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;
}