# 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


- loop = 110


#### LEIBNIZ
- loop = 30




#### EULER
- loop = 25


- loop = 100


## 證明 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)

- 補充其中的 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/)