2019q3 Homework4 (compute-pi)

# 2019q3 Homework4 (compute-pi) contributed by < `colinyoyo26` > ## 開發環境 - Linux 4.15.0-47-generic Ubuntu 16.04.1 - CPU Architecture: x86_64 L1d cache： 32K L1i cache： 32K L2 cache： 256K L3 cache： 3072K NUMA node0 CPU(s): 0-3 CPU(s): 4 threads per core: 2 cpu cores : 2 ## Wall-clock time vs. CPU time - Wall-clock time - 實際經握的時間 - 可能與其他網路節點同步，所以不一定是 monotonic - CPU time - 程序在 CPU 上面運行消耗的時間 ## 實驗 :::info loop 設定在特定數值時編譯會產生為定義參考到的問題，有待釐清 ::: - 以下四種方法可以看到 OenpMP(4threads) 幾乎都會有波動 - 推測是我邊使用瀏覽器的關係，因為瀏覽器也用到多個執行緒，影響到使用 cpu core 的時間 - 實際測試不使用瀏覽器的情況，曲線如預期變得非常平滑 - 最不解的是為什麼 LEIBNIZ 使用 SIMD 的效能會最差 #### 離散積分 - loop = 25 ![](https://i.imgur.com/dCC2uAb.png) ![](https://i.imgur.com/2GxsoQ9.png) - loop = 110 ![](https://i.imgur.com/dwzuL90.png) ![](https://i.imgur.com/zZgBazb.png) #### LEIBNIZ - loop = 30 ![](https://i.imgur.com/Y7Iv3VK.png) ![](https://i.imgur.com/8ahoGU0.png) ![](https://i.imgur.com/MKwcrrc.png) ![](https://i.imgur.com/BGwiCU8.png) #### EULER - loop = 25 ![](https://i.imgur.com/p4vI9fi.png) ![](https://i.imgur.com/mq0nbmV.png) - loop = 100 ![](https://i.imgur.com/zItNGpg.png) ![](https://i.imgur.com/rB29VaW.png) ## 證明 Leibniz formula for π - 維基百科的推導比較 tricky 一點，我認為用泰勒展開比較直觀，推導如下 - 也許看起來泰勒展開逼近不夠嚴謹，但事實上這個函式的泰勒展開會等於原函式(解析型函數)，直接用長除法也會得到一樣的結果 - 泰勒展開式概念是用無窮級數去逼近函式 - 參考點的 N 階導數都等於原函式 $for$ $all$ $N = 0, 1, 2, 3......$ - 求每項係數的方法有點類似傅立葉係數，差在傅立葉是利用正交性消掉其他項，泰勒展開是微分到常數項，然後在帶入參考點消掉其他項 $according\ to\ Taylor\ expansion\\ \dfrac{1}{1+x^2} = \sum_{k=0}^{\infty}{(-x^2)^{i}}$ $then\\ \begin{split}\dfrac{\pi}{4} &= \tan^{-1}{(1)} \\ &= \int_{0}^{1}{\dfrac{dx}{1+x^2}}\\ &= \int_{0}^{1}{\sum_{k=0}^{\infty}{(-x^2)^{k}}{dx}}\\ &= \sum_{k=0}^{\infty}{\dfrac{(-1)^{k}(1)^{2k + 1}}{2k +1}} - \sum_{k=0}^{\infty}{\dfrac{(-1)^{k}(0)^{2k + 1}}{2k +1}}\\ &= \sum_{k=0}^{\infty}{\dfrac{(-1)^{k}}{2k +1}}\end{split}$ 延伸閱讀：[Leibniz formula for π](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80) ## RABINOWITZ AND WAGON’S SPIGOT ALGORITHM - This algorithm is a “spigot” algorithm: it pumps out digits one at a time and does notuse the digits after they are computed >這個演算法一次算出固定幾個 digits 計算後不再用到 - 基本概念 - 乘以基數的 n 次方把數字往左推 - 每次都輸出固定幾個 digits - 餘數部分留在公式內 - 接者在用相同邏輯運算剩下的 digits - 延伸閱讀: [ASpigot Algorithm for the Digits of Pi](http://www.cs.williams.edu/~heeringa/classes/cs135/s15/readings/spigot.pdf) - 以下參考: [Computing Pi in C](https://crypto.stanford.edu/pbc/notes/pi/code.html) - 可以搭配上面延伸閱讀第五頁的地方看 - 基本概念 - 數學式子先化成 Horner's form - 由又到左進行除法， quotion 進位到左邊括號， remainder留在原位 - 一直往左運算，就能得到左邊的 digits - 接者在用相同邏輯運算剩下的 digits $\dfrac{\pi}{2} = \dfrac{1}{2} \sum_{k=0}^{\infty}{\dfrac{(n!)^22^{n+1}}{(2n+1)!}}$ - 寫成 Horner's form $\dfrac{\pi}{2} = 1 + \dfrac{1}{3}{(1 + \dfrac{2}{5}(1+\dfrac{3}{7}(1+\dfrac{4}{9}(1+.....))))}$ - 根據以上式子 Dik T. Winter 寫成以下 C program ```cpp int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5; for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a, f[b]=d%--g,d/=g--,--b;d*=b);} ``` - 可以重寫成以下比較易讀的形式 - 我加了一些註解請搭配以下推導看 ```cpp #include <stdio.h> int main() { int r[2800 + 1]; int i, k; int b, d; int c = 0; for (i = 0; i < 2800; i++) { /* r[i] 為第 i 個左括號和第 i + 1 個左括號間的值 */ r[i] = 2000; } for (k = 2800; k > 0; k -= 14) { d = 0; i = k; for (;;) { d += r[i] * 10000; /* 取 bi */ b = 2 * i - 1; /* 更新 ri */ r[i] = d % b; /* 求 qi */ d /= b; i--; if (i == 0) break; /* qi * ai */ d *= i; } /*c + d / 10000 = q0*/ printf("%.4d", c + d / 10000); /* 存 r0 的值 */ c = d % 10000; } return 0; } ``` - 這個 program 的目標做以下運算 - 補充一下，以上 code 的邏輯一開始就把 2 乘進去，所以一開始才有 `r[i] = 2000;`，其中為3位數是為了精度的關係 $\pi = 2(1 + \dfrac{1}{3}{(1 + \dfrac{2}{5}(1+\dfrac{3}{7}(1+\dfrac{4}{9}(1+.....(1+\dfrac{2799}{{2\cdot2799}+1})))))})$ - 先從第一個 iteration 開始，看到 `d += r[i] * 10000` - 其中 `r[i]` 初始值為 2000 - 把兩邊同乘 ${10^7}$ 看 ${10^7}\pi = 2({10^7} + \dfrac{1}{3}{({10^7} + \dfrac{2}{5}({10^7}+\dfrac{3}{7}({10^7}+\dfrac{4}{9}({10^7}+....({10^7}+\dfrac{2799}{{2799}+1})))))})$ - 訂以下符號 $\dfrac{1}{3} = \dfrac{a_1}{b_1}$ $.....$ $\dfrac{2799}{2\cdot{2799}+1} = \dfrac{a_n}{b_n}$ - 根據除法原理 - 請由內到外對照每一項的除法 ${10^7} = q_n\cdot b_n + r_n where\ {0 <= r_n < {2\cdot b_n+1}}$ ${10^7 + a_n \cdot q_n} = q_{n-1}\cdot b_{n-1} + r_{n-1} where\ {0 <= r_{n-1} < {10^7 + q_n\cdot b_n}}$ > 1. 注意取下高斯 $r_n$ 去掉 > 2. 這邊要乘 $a_n$ 才對，原網頁寫錯了 - 所以把遞迴全部展開，並把下高斯去掉的項補回去可以得到 ${10^7}\pi = 2({10^7} + a_1\cdot q_1) + 2(0 + \dfrac{1}{3}{(r_1 + \dfrac{2}{5}(r_2+\dfrac{3}{7}({10^7}+\dfrac{4}{9}(r_3+....(r_n+\dfrac{2799}{{2799}+1})))))})$ - 我們就得到了前四項的值 $Let\ q_0 = 2({10^3} + \dfrac{a_1\cdot q_1}{10^4})$ ${10^3}\pi \approx q_0$ - 下個 iteration 用相同邏輯繼續求下去 $Let r_0 = 2({10^7} + a_1\cdot q_1)(mod 10^4)$ ${10^7}\pi - q_0 = r_0 + 2(0 + \dfrac{1}{3}{(r_1 + \dfrac{2}{5}(r_2+\dfrac{3}{7}({10^7}+\dfrac{4}{9}(r_3+....(r_n+\dfrac{2799}{{2799}+1})))))})$ > 接者用同樣方法求 ${10^7}\pi - q_0$ 的前四項 ## BBP - BBP 公式為 spitgot algorithm 可以求出 $\pi$ 十六進位下第 n 個小數 $\pi = 4 \sum_{k=0}^{\infty}{\dfrac{1}{16^k(8k+1)}}-2\sum_{k=0}^{\infty}{\dfrac{1}{16^k(8k+4)}}-\sum_{k=0}^{\infty}{\dfrac{1}{16^k(8k+5)}}-\sum_{k=0}^{\infty}{\dfrac{1}{16^k(8k+8)}}$ - 把前幾項推到整數位，以第一項為例 $\sum_{k=0}^{\infty}{\dfrac{16^{n-k}}{16^k(8k+1)}} = \sum_{k=0}^{n}{\dfrac{16^{n-k}}{16^k(8k+1)}} + \sum_{k=n+1}^{\infty}{\dfrac{16^{n-k}}{16^k(8k+1)}}$ - 其中最右項在加上前面項的餘數就是小數部份，進位的情況可以靠多計算後面幾項解決 $\sum_{k=0}^{n}{\dfrac{16^{n-k}(\mod 8k+1)}{16^k(8k+1)}} + \sum_{k=n+1}^{\infty}{\dfrac{16^{n-k}}{16^k(8k+1)}}$ - David H. Bailey 2006 年寫的 [code](https://www.experimentalmath.info/bbp-codes/) - 從中改寫並用 openMP 平行化做效能分析 - [github](https://github.com/colinyoyo26/compute-pi/tree/master/BBP) ![](https://i.imgur.com/lCOcVR5.png) - 補充其中的 modulo function 可以參考 [Modular exponentiation](https://en.wikipedia.org/wiki/Modular_exponentiation) $16^k \mod m$ - 其中從二進制思考 k 假設 k = 101b 則 $16^k \mod m = 16^{100b} \cdot 16^{00b} \cdot 16^{1b} \mod m$ - 進一步化簡 $16^k \mod m = ((16^{100b}) \mod m) \cdot ((16^{1b}) \mod m) \mod m$ - 實作是從 LSB 往 MSB 檢查 k 每個 bit 是 0 or 1 ## 參考資料 - [twzjwang 的共筆](https://hackmd.io/@twzjwang/ByYgvi-qg?type=view) - [0xff07 的共筆](https://hackmd.io/EAiy_D9SQESm-TlFdGwl4g?both) - [laochanlam 的共筆](https://hackmd.io/OeG5AEGXRj-2MrN7I1ChHA?both) - [Leibniz formula for π](https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80) - [Makefile 教學](https://kezeodsnx.pixnet.net/blog/post/25125766-makefile?source=---------26------------------) - [Makefile中的 $ 和 $$](https://blog.csdn.net/qq_34369618/article/details/52973720) - [Pi Formulas, Algorithms and Computations](https://bellard.org/pi/) - [pi_bin](https://bellard.org/pi/pi_bin.pdf) - [Computing Pi in C](https://crypto.stanford.edu/pbc/notes/pi/code.html) - [Spigot algorithm](https://www.cut-the-knot.org/Curriculum/Algorithms/SpigotForPi.shtml#RW) - [簡易的程式平行化－OpenMP（二）語法說明](https://cometlc.pixnet.net/blog/post/2596845) - [BBP Code Directory](https://www.experimentalmath.info/bbp-codes/)