派值

這篇文章是關於$\pi$值的計算。

計算

九進制

考慮公式
$$
\pi = \frac{\sqrt{3}}{2}\left[ \Phi(\frac{1}{9},1,\frac{1}{4}) - \frac{1}{3}\Phi(\frac{1}{9},1,\frac{3}{4}) \right]
$$

其中$\Phi(z,s,u)$是Lerch Transcendent Function,寫成級數和的形式即為
$$
\pi = \sum_{n = 0}^\infty\frac{1}{9^n}\left[ \frac{16\sqrt{3}}{3}\frac{n + 1}{16n^2+16n+3} \right]
$$
這是九進制下的Bailey-Borweim-Plouffe Formula,類似地,可以用Spigot Algorithm去計算九進制下指定位$\pi$的值,這個算法的優點在於不需要計算該位之前的所有小數。

可以把上述$\pi$的九進制下的級數和表示$f$寫為
$$
f = \sum_{n = 0}^\infty \frac{f_n}{9^n}
$$
$f_n$是一個零到八之間的整數,表示九進制下小數點後第$n$位的值,假設只對該位小數感興趣,可以用這個公式
$$
f_n = \lfloor 9 \cdot \{f\cdot 9^{n-1}\} \rfloor
$$
其中$\lfloor \cdot \rfloor$表示向下取整,${\cdot}$表示取小數部分。這裡的計算需要用到開頭提到的公式,用$\tilde f$表示近似值。原則上$n$不需要計算到無窮階,假設只需要計算到第$\kappa$階為止,
$$
\tilde f = \sum_{n = 0}^\kappa\frac{1}{9^n}\left[ \frac{16\sqrt{3}}{3}\frac{n + 1}{16n^2+16n+3} \right]
$$
假設下面剩餘項之和不大於一個小量$\epsilon_n$
$$
\sum_{i = 1+\kappa}^\infty \frac{9^n}{9^i}\cdot \left[ \frac{16\sqrt{3}}{3}\frac{i + 1}{16i^2+16i+3} \right]< \epsilon_n
$$
估計左邊的上界
$$
\begin{equation}
\begin{aligned}
\sum_{i = 1+\kappa}^\infty \frac{9^n}{9^i}\cdot \left[ \frac{16\sqrt{3}}{3}\frac{i + 1}{16i^2+16i+3} \right] =& \sum_{i = 1+\kappa}^\infty \frac{9^n}{9^{i}}\frac{\sqrt{3}}{3}\frac{1}{i + \frac{3}{16}\frac{1}{i+1}}\\
<& \sum_{i = 1+\kappa}^\infty \frac{9^n}{9^{i}}\frac{\sqrt{3}}{3}\frac{1}{1+\kappa}\\
\leq& \frac{\sqrt{3}}{24}9^{n-\kappa}
\end{aligned}
\end{equation}
$$
令上界不大於該小量,則
$$
\ln\frac{\sqrt{3}}{24} + (n-\kappa)\ln9 \leq \ln\epsilon_n
$$
考慮$\kappa$是整數,可以取
$$
\kappa_n = \lceil\frac{n\ln9 - \ln(8\sqrt{3}) - \ln\epsilon_n}{\ln9}\rceil
$$
$\lceil\cdot\rceil$表示向上取整,只需要計算到$\kappa_n$那一階為止。此外小量可以取
$$
\epsilon_n = 1-\{9 \cdot \{\tilde f\cdot 9^{n-1}\}\}
$$
即小量與當前狀態有關,可以在程序迭代的每一步都計算一個新的$\epsilon_n$和$\kappa_n$,當算出來的$\kappa_n$不大於當前迭代所在的階時,則可以終止迭代而輸出結果。

十進制

現在需要將其轉化為十進制,用$f_k’$表示十進制下小數點後第$k$位的值(零到九之間的整數),因此也可以寫成
$$
f = \sum_{k = 0}^\infty \frac{f_k’}{10^k}
$$
的形式。若只對十進制下第$k$位小數感興趣,則類似有如下公式
$$
f_k’ = \lfloor 10 \cdot \{f\cdot 10^{k-1}\} \rfloor
$$
這裡需要將九進制下得到的結果轉化為十進制,同上面一樣用$\tilde f$表示近似值
$$
\tilde f = \sum_{k = 0}^\kappa \frac{f_k}{9^k}
$$
假設可以略去九進制級數和中$\kappa$階以上的項
$$
\sum_{i = 1+\kappa}^\infty \frac{f_i}{9^i}\cdot 10^k < \epsilon_k
$$
考慮左邊的上界
$$
\sum_{i = 1+\kappa}^\infty \frac{f_i}{9^i}\cdot 10^k < \sum_{i = 1+\kappa}^\infty \frac{8}{9^i}\cdot 10^k = 10^k \cdot 9^{-\kappa}
$$
要求上界不大於一個小量$\epsilon_k$。得到
$$
\kappa \geq \frac{k\ln10-\ln\epsilon_k}{\ln9}
$$
因為$\kappa$是整數,可以記
$$
\kappa_k = \lceil\frac{k\ln10-\ln\epsilon_k}{\ln9}\rceil
$$
即級數求和運算只需要計算到$\kappa_k$那階為止,九進制下高於該階的項對十進制下的第$k$位小數的貢獻可以略去。同樣小量可以取

$$
\epsilon_k = 1-\{10 \cdot \{\tilde f\cdot 10^{k-1}\}\}
$$
類似可以在迭代的每一步都計算一個新的$\epsilon_k$和$\kappa_k$,當算出來的$\kappa_k$不大於當前迭代所在的階時,可以終止迭代而輸出結果。

後記

之所以考慮小量是因為可能涉及到進位。此外在計算小數時,先取模再做除法,大數取模時注意利用這個性質
$$
a \cdot b \equiv a_0\cdot b_0 \mod m
$$
其中$a \equiv a_0 \mod m$以及$b \equiv b_0 \mod m$,特別在處理冪次時。

後續會補上關於$\pi$等式的證明。

補充

證明等式之前,考慮下面兩個積分
$$
I_1 = \int_0^1dx(1 + ax)^{-1/4}(1 + bx)^{-3/4}
$$
以及
$$
I_2 = \int_0^1dx(1 + ax)^{-3/4}(1 + bx)^{-1/4}
$$
令$t^4 = (1 + ax) / (1 + bx)$,則$x = (1 - t^4)/(-a + bt^4)$,$dx/dt = -4t^3(-a+b)/(-a + bt^4)^2$。以及$1 + ax = (-a+b)t^4/(-a+bt^4)$,$1 + bx = (-a+b)/(-a+bt^4)$。得到
$$
I_1(a,b) = \int_1^{\tau} dt \frac{4t^2}{a - bt^4}
$$
其中$\tau^4 = (1 + a) / (1 + b)$,以及
$$
I_2(a,b) = \int_1^\tau dt \frac{4}{a - bt^4}
$$
這裡考慮它們的線性組合
$$
I(a,b) = c_1I_1(a,b) + c_2I_2(a,b) = 4\int_0^\tau dt \frac{c_1t^2 + c_2}{a - bt^4}
$$
如果$ab > 0$,則
$$
\frac{c_1t^2 + c_2}{a - bt^4} = \frac{1}{a}\left[\frac{t_1}{1-\sqrt{b/a}t^2}+\frac{t_2}{1+\sqrt{b/a}t^2}\right]
$$
其中$t_1$和$t_2$是待定係數,代入可得
$$
\begin{align}
c_2 =& t_1 + t_2\\
c_1 =& \sqrt{b/a}(t_1 - t_2)
\end{align}
$$
現在挑選係數使得$t_1 = 0$,則$c_1 / c_2 = -\sqrt{b/a}$。這樣積分為
$$
I(a,b) = \frac{4t_2}{a\sqrt[4]{b/a}}\left[\text{atan}(\sqrt[4]{b/a}\tau) - \text{atan}(\sqrt[4]{b/a})\right]
$$
通過設定$a$和$b$的大小,就可以得到與$\pi$相關的結果,譬如可以取$-1/9$和$-1$。

之所以考慮這兩個積分,原因在於根據定義,Lerch Transcendent Function
$$
\Phi(z,s,a) \equiv \sum_{n = 0}^\infty\frac{z^n}{(a+n)^s}
$$
在某些情況下可以寫成Hypergeometric Function的形式
$$
\ _2F_1(a,b;c;z) \equiv \sum_{n=0}^\infty\frac{(a)_n(b)_n}{(c)_n}\frac{z^n}{n!} = 1 + \frac{ab}{1!c}z+\frac{a(a+1)b(b+1)}{2!c(c+1)}z^2 + …
$$
其中$(\cdot)_n$是Rising Factorial符號。它有積分表達式,但要求$\left|z\right| < 1$以及$0 < \mathcal{Re}(b) < \mathcal{Re}(c)$
$$
\ _2F_1(a,b;c;z) = \frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\int_0^1dt\frac{t^{b-1}(1-t)^{c-b-1}}{(1-tz)^a}
$$
其中$\Gamma(\cdot)$是伽馬函數。回到開頭的問題,根據上述討論可以得到
$$
\Phi(\frac{1}{9},1,\frac{1}{4}) = \sum_{n = 0}^\infty\frac{1}{9^n}\frac{1}{\frac{1}{4}+n} = 4\sum_{n=0}\frac{(\frac{1}{4})_n(1)_n}{(\frac{5}{4})_n}\frac{1}{n!} = 4\ _2F_1(\frac{1}{4},1;\frac{5}{4};\frac{1}{9})
$$
整理成積分的形式
$$
\Phi(\frac{1}{9},1,\frac{1}{4}) = \frac{4\Gamma(\frac{5}{4})}{\Gamma(1)\Gamma(\frac{1}{4})}\int_0^1dt(1-\frac{t}{9})^{-\frac{1}{4}}(1-t)^{-\frac{3}{4}} = I_1(-\frac{1}{9},-1)
$$
同理可得
$$
\begin{equation}
\begin{aligned}
\Phi(\frac{1}{9},1,\frac{3}{4}) =& \frac{4}{3}\ _2F_1(\frac{3}{4},1;\frac{7}{4};\frac{1}{9})\\
=& \frac{4\Gamma(\frac{7}{4})}{3\Gamma(1)\Gamma(\frac{3}{4})}\int_0^1dt(1-\frac{t}{9})^{-\frac{3}{4}}(1-t)^{-\frac{1}{4}}\\
=& I_2(-\frac{1}{9},-1)
\end{aligned}
\end{equation}
$$

現在回到$I(a,b)$,代入具體值後得到
$$
I(-\frac{1}{9},-1) = -\frac{6\pi}{\sqrt{3}}t_2
$$
令其等於$\pi$,解出$t_2 = -\frac{\sqrt{3}}{6}$,因此$c_1 = \frac{\sqrt{3}}{2}$,$c_2 = -\frac{\sqrt{3}}{6}$,則
$$
\pi = c_1I_1(-\frac{1}{9},-1) + c_2I_2(-\frac{1}{9},-1) = \frac{\sqrt{3}}{2}\Phi(\frac{1}{9},1,\frac{1}{4}) -\frac{\sqrt{3}}{6}\Phi(\frac{1}{9},1,\frac{3}{4})
$$
這就是文章開頭提到的線性組合。