Arbitrarily accurate rational approximations of $\pi$ with explicit error bound

Question

I would like a method for computing $\pi$ for which I have an explicit bound on the error, and that has the property that if I do the calculation up to a certain accuracy, but then need an arbitrarily finer accuracy later, I can pick up the calculation process where I left off rather than having to start over from scratch. In light of the second requirement, I don't want to have to deal with, say, square roots, because whatever accuracy I calculate the roots to in order to reach some desired accuracy for $\pi$, there is some finer accuracy for which the root approximations already used would have needed to be more accurate. I believe I'll want to stick to algorithms for which the operations required are limited to arithmetic on rationals in order to avoid such issues.

The Maclaurin series for $4\arctan(1)$ is an example of an algorithm that obeys these properties: $$\pi=4\arctan(1)=4\sum_{n=0}^\infty\frac{(-1)^{n-1}}{2n+1}$$ Provided I save the value of $n$ for the last term I added, it's easy to resume the calculation, and I always know that the magnitude of the error is less than that of the next term.

However, I know that this algorithm converges very slowly. What contenders are there for faster algorithms that obey the desired constraints? It sounds like an unbounded spigot algorithm would work, but if I'm not mistaken, those are subject to more constraints than I care about, so maybe it's possible to do better.

Answer 1

Below follows an analysis that will allow you to compute $\pi$ rapidly using only basic arithmetic operations. If the initial guess for $\pi$ is rational, then all numbers are rational. The function $x \rightarrow \sin(x)$ appears early on, but a suitable approximation is derived at the end.

Fixed point iterations allow you to refine an approximation. Consider the iteration given by $$x_{n+1} = g(x_n),$$ where $$g(x) = x + \sin(x)$$ and $$x_0 \approx \pi$$ will be chosen later. It is easy to see that $\pi$ is a fixed point for $g$, i.e., $$g(\pi) = \pi.$$ We will now show that $x_n$ convergences to $\pi$ cubically, provided that $x_0$ is sufficiently close to $\pi$. Let $x \in \mathbb{R}$ be given. By Taylor's formula, there exists at least one $\xi$ between $x$ and $\pi$, such that $$ g(x) - \pi = g(x) - g(\pi) = g'(\pi)(x - \pi) + \frac{g''(\pi)}{2}(x - \pi)^2 + \frac{g^{(3)}(\xi)}{6}(x-\pi)^3.$$ However, since $$g'(x) = 1 +\cos(x), \quad g''(x) = - \sin(x),$$ this reduces to $$ g(x) - \pi = \frac{g^{(3)}(\xi)}{6}(x-\pi)^3.$$ It follows that $$ |\pi - g(x)| \leq \frac{1}{6}|\pi-x|^3,$$ because $$g^{(3)}(t) = - \cos(t)$$ is bounded by unity everywhere. In terms of the fixed point iteration we have $$|\pi - x_{n+1}| \leq \frac{1}{6} |\pi - x_n|^3.$$ By induction, we discover that $$ |\pi - x_n| \leq \left(\frac{1}{6}\right)^{m(n)} |\pi - x_0|^{3^n}, \quad m(n) = \frac{3^n - 1}{2}.$$ Convergence is assured if $x_0$ is chosen such that $$|\pi - x_0|^3 < 1.$$

This leaves us with the task of computing $x \rightarrow \sin(x)$ for $x$ close to $\pi$. For this problem we will utilize the trigonometric identity $$\sin(3\theta) = 3 \sin(\theta) - 4 \sin(\theta)^3.$$ Let $x$ be close to $\pi$ and let $k$ be a large integer, such $\theta_0 = x/3^k$ is close to $0$ and $$ \sin(\theta_0) \approx \theta_0 - \frac{1}{6} \theta_0^3 $$ is an acceptable approximation. Now let $\theta_j = 3 \theta_{j-1}$. Then by design $\theta_k = 3^k \theta_0 = x$ and since $$ \sin( \theta_j) = 3 \sin(\theta_{j-1}) - 4 \sin(\theta_{j-1})^3$$ we have an iteration for computing $x \rightarrow \sin(x)$.

Answer 2

I am not sure if you like Ramanujan’s rapidly-converging series, since there is one square root...

$$\frac1\pi=\frac{2\sqrt2}{9801}\sum_{k=0}^\infty \frac{(4k)!(1103+26390k)}{(k!)^4 396^{4k}}$$

If you really hate the $\sqrt2$, try some approximation of it, e.g. $$\sqrt2\approx \frac{577}{408}$$

ADDED:

To make the higher order terms look like a geometric series, I can derive a very loose error bound.

Using Stirling’s approximation: $$\frac{(4k)!}{(k!)^4}\approx \frac1{\sqrt2}\frac1{\pi^{3/2}}\frac{4^{4k}}{k^{3/2}}<\frac1{\sqrt2}\frac1{\pi^{3/2}}\frac{4^{4k}}{k}$$ Also, $$1103+26390k<26391k$$

Thus, the error $E$ for $\frac1\pi$ is bounded by $$E<C\sum^{\infty}_k\frac1{99^k}=\frac{C}{98\cdot99^{k-1}}$$ where $k$ is the number of terms you take and $$C=\frac{17594}{3267\pi^{3/2}}\approx 1$$ Or, $$E\approx\frac1{99^k}\approx100^{-k}$$

Answer 3

I've decided to use the BBP formula, linked by Alex D in the comments: $$\pi=\sum_{k=0}^\infty\frac{R(k)}{16^k}\text{ where }R(k)=\frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6}$$ I get an error bound by noting that $$R(k)>\frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+4}-\frac{1}{8k+4}=\frac{4}{8k+1}-\frac{4}{8k+4}>0$$ so each term of the series is positive and the sequence of partial sums increases monotonically. Taking for granted that the series does in fact converge to $\pi$, this means that every partial sum is an underestimate, and the error of the $n$th partial sum is equal to its tail. $$R'(k)=\frac{16}{(8x+4)^2}+\frac{8}{(8x+5)^2}+\frac{8}{(8x+6)^2}-\frac{32}{(8x+1)^2}=-\frac{122880x^5+354816x^4+405504x^3+228616x^2+63328x+6853}{(2x+1)^2(4x+3)^2(8x+1)^2(8x+5)^2}$$ which is negative for all $k\ge0$, so $R(k)$ strictly decreases over that range, and we have this for the tail of the $n$th partial sum: $$\sum_{k=n+1}^\infty\frac{R(k)}{16^k}<R(1)\sum_{k=n+1}^\infty\frac{1}{16^k}=\frac{106}{819}\left(\sum_{k=0}^\infty\frac{1}{16^k}-\sum_{k=0}^n\frac{1}{16^k}\right)=\frac{106}{819}\left(\frac{1}{1-\frac{1}{16}}-\frac{1-\frac{1}{16^n}}{1-\frac{1}{16}}\right)=\frac{106}{16^{n-1}*12285}$$

Answer 4

If you want to find a formula that converges to $\pi$ faster than Leibnitz's formula$$\pi=4\sum\limits_{k\geq0}\frac {(-1)^k}{2k+1},$$then just perform an Euler Transform on the infinite sum to get a faster formula for $\pi$ in terms of the factorial and double factorial.$$\pi=2\sum\limits_{k\geq0}\frac {k!}{(2k+1)!!}$$

First off, the transform gives us$$\frac {\pi}4\sum\limits_{n\geq0}\frac 1{2^{n+1}}\sum\limits_{m=0}^n\binom nm\frac {(-1)^m}{2m+1}$$The inner sum can be simplified by expanding the right-hand side's inner finite sum and realizing that it takes the form $2^n\frac {n!}{(2n+1)!!}$.$$\begin{align*}\sum\limits_{m=0}^0\binom 0m\frac {(-1)^m}{2m+1} & =2^0\cdot\frac {0!}{1!!}\\\sum\limits_{m=0}^1\binom 1m\frac {(-1)^m}{2m+1} & =2^1\cdot\frac {1!}{3!!}\\\sum\limits_{m=0}^2\binom 2m\frac {(-1)^m}{2m+1} & =2^2\cdot\frac {2!}{5!!}\\\cdots\qquad\cdots\qquad\cdots & \qquad\cdots\end{align*}$$Going on and simplifying the fraction, it's evident that we can rewrite the formula as$$\frac {\pi}2=\sum\limits_{n\geq0}\frac {n!}{(2n+1)!!}$$