A Rational Approximation of the Two-Term Machin-Like Formula for π

1. Preliminaries

In 1876, English astronomer and mathematician John Machin demonstrated an efficient method to compute digits of $\pi$ by using his famous discovery [1,2,3,4]

$$\frac{\pi }{4}=4arctan\left(\frac{1}{5}\right)-arctan\left(\frac{1}{239}\right).$$

(1)

In particular, due to relatively rapid convergence of this formula, he was the first to compute more than 100 digits of $\pi$. Nowadays, the equations of kind

$$\frac{\pi }{4}=\sum _{j=1}^J{A}_{j}arctan\left(\frac{1}{B_j}\right),$$

where $A_j$ and $B_j$ are either integers or rational numbers, are named after him as the Machin-like formulas for $\pi$ [1,2,3,4].

Theorem 1 below shows the arctangent formula (2) for $\pi$. We can use this equation as a starting point to generate the-two term Machin-like formula for $\pi$ of kind

$$\frac{\pi }{4}=2^{k-1}arctan\left(\frac{1}{x}\right)+arctan\left(\frac{1}{y}\right),$$

where $k,x\in \mathbb{N}$ and $y\in \mathbb{Q}$. Further, we will show the significance of the multiplier $2^{k-1}$ in this formula.

Theorem 1. 

There is a formula for π at any integer $k\ge 1$ [5]:

$$\frac{\pi }{4}=2^{k-1}arctan\left(\frac{\sqrt{2-a_{k-1}}}{a_k}\right),$$

(2)

where $a_0=0$ and $a_k=\sqrt{2+a_{k-1}}$ are nested radicals.

Proof of Theorem 1. 

Since

$$cos\left(\frac{\pi }{2^2}\right)=\frac{1}{2}\sqrt{2}=\frac{1}{2}{a}_1,$$

from the identity

$$cos\left(2x\right)=2{cos}^2\left(x\right)-1,$$

it follows, by induction, that

$$cos\left(\frac{\pi }{2^3}\right)=\frac{1}{2}\sqrt{2+\sqrt{2}}=\frac{1}{2}{a}_2,$$

$$cos\left(\frac{\pi }{2^4}\right)=\frac{1}{2}\sqrt{2+\sqrt{2+\sqrt{2}}}=\frac{1}{2}{a}_3,$$

$$cos\left(\frac{\pi }{2^5}\right)=\frac{1}{2}\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2}}}}=\frac{1}{2}{a}_4,$$

$$\begin{gathered} ⋮ \\ cos\left(\frac{\pi }{2^{k+1}}\right)=\frac{1}{2}\underset{k\mathrm{square}\mathrm{roots}}{\underbrace{\sqrt{2+\sqrt{2+\sqrt{2+\cdots +\sqrt{2}}}}}}=\frac{1}{2}{a}_k. \end{gathered}$$

Therefore, using

$$sin\left(\frac{\pi }{2^{k+1}}\right)=\sqrt{1-{cos}^2\left(\frac{\pi }{2^{k+1}}\right)},$$

we obtain

$$\frac{\pi }{2^{k+1}}=arctan\left(\frac{\sqrt{1-{cos}^2\left(\frac{\pi }{2^{k+1}}\right)}}{cos\left(\frac{\pi }{2^{k+1}}\right)}\right)=arctan\left(\frac{\sqrt{2-a_{k-1}}}{a_k}\right)$$

and Equation (2) follows.  □

Lemma 1 and its proof below show how Equation (2) is related to the well-known limit (3) for $\pi$.

Lemma 1. 

There is a limit such that [6]

$$\pi =\underset{k\to \infty }{lim}{2}^k\sqrt{2-\underset{k-1\mathrm{square}\mathrm{roots}}{\underbrace{\sqrt{2+\sqrt{2+\sqrt{2+\cdots +\sqrt{2}}}}}}}.$$

(3)

Proof of Lemma 1. 

Let

$$\underset{k\to \infty }{lim}{a}_k=x.$$

Then, we can write

$$\underset{k\to \infty }{lim}{a}_k=\underset{k\to \infty }{lim}\sqrt{2+a_{k-1}}=\sqrt{2+{lim}_{k\to \infty }{a}_k}$$

or

$$x=\sqrt{2+x}$$

or

$$x^2-x-2=0.$$

Solving this quadratic equation leads to two solutions $x=2$ and $x=-1$. Since $a_k$ cannot be negative, we came to conclusion that

$$\underset{k\to \infty }{lim}{a}_k=2.$$

Consequently, this gives

$$\underset{k\to \infty }{lim}\sqrt{2-a_{k-1}}=\sqrt{2-{lim}_{k\to \infty }{a}_{k-1}}=\sqrt{2-{lim}_{k\to \infty }{a}_k}=0.$$

Since the formula (2) is valid for any arbitrarily large k, we can write

$$\frac{\pi }{4}=\underset{k\to \infty }{lim}{2}^{k-1}arctan\left(\frac{\sqrt{2-a_{k-1}}}{a_k}\right)$$

As $arctan\left(x\right)\to x$ when $x\to 0$, we have

$$\begin{array}{cc}\hfill \frac{\pi }{4}& =\underset{k\to \infty }{lim}{2}^{k-1}\frac{\sqrt{2-a_{k-1}}}{a_k}=\underset{k\to \infty }{lim}{2}^{k-1}\sqrt{2-a_{k-1}}/\underset{k\to \infty }{lim}{a}_k\hfill \\ & =\underset{k\to \infty }{lim}{2}^{k-2}\sqrt{2-a_{k-1}}=\underset{k\to \infty }{lim}{2}^{k-1}\sqrt{2-a_k}\hfill \end{array}$$

and Equation (3) follows. This completes the proof. □

Lemma 2 and its proof below show how Equation (2) can be transformed into the two-term Machin-like formula for $\pi$ [7].

Lemma 2. 

In the following equation

$$\frac{\pi }{4}=2^{k-1}arctan\left(\frac{1}{⌊a_k/\sqrt{2-a_{k-1}}⌋}\right)+arctan\left(\frac{1}{{\beta }_k}\right),$$

(4)

the value ${\beta }_k$ is always a rational number.

Proof of Lemma 2. 

It is convenient to define

$${\alpha }_k=⌊a_k/\sqrt{2-a_{k-1}}⌋.$$

(5)

Using the identity

$$arctan\left(\frac{1}{x}\right)=\frac{1}{2i}ln\left(\frac{x+i}{x-i}\right)$$

and taking into consideration that

$$\frac{\pi }{4}=arctan\left(1\right),$$

Equation (4) can be represented as

$$\frac{\pi }{2}i=ln{\left(\frac{{\alpha }_k+i}{{\alpha }_k-i}\right)}^{2^{k-1}}ln\left(\frac{{\beta }_k+i}{{\beta }_k-i}\right)$$

or

$$i={\left(\frac{{\alpha }_k+i}{{\alpha }_k-i}\right)}^{2^{k-1}}\left(\frac{{\beta }_k+i}{{\beta }_k-i}\right).$$

(6)

It is not difficult to see by substitution that the following formula [7]

$${\beta }_k=\frac{2}{{\left(\frac{{\alpha }_k+i}{{\alpha }_k-i}\right)}^{2^{k-1}}-i}-i$$

(7)

is a solution of the Equation (6). Since k is a positive integer greater than 1, we can see that the real and imaginary parts of the expression

$${\left(\frac{{\alpha }_k+i}{{\alpha }_k-i}\right)}^{2^{k-1}}={\left(\frac{{\alpha }_k^2-1}{{\alpha }_k^2+1}+i\frac{2{\alpha }_k}{{\alpha }_k^2+1}\right)}^{2^{k-1}}$$

cannot be an irrational number if ${\alpha }_k$ is an integer. This means that the real and imaginary parts of the value ${\beta }_k$ must be both rational. Since the first arctangent term of Equation (4) is a real number, the second arctangent term is also a real number. Therefore, we conclude that the imaginary part of the value ${\beta }_k$ is equal to zero. This completes the proof. □

This is not the only method to generate the two-term Machin-like formulas for $\pi$ of kind (4). Recently, Gasull et all. proposed a different method to derive the same equation (see Table 2 in [8]). Lemma 3 and its proof below show how the Equation (4) can be represented in a trigonometric form [7].

Lemma 3. 

Equation (4) can be expressed as

$$\frac{\pi }{4}=2^{k-1}arctan\left(\frac{1}{{\alpha }_k}\right)+arctan\left(\frac{1-sin\left(2^{k-1}arctan\left(\frac{2{\alpha }_k}{{\alpha }_k^2-1}\right)\right)}{cos\left(2^{k-1}arctan\left(\frac{2{\alpha }_k}{{\alpha }_k^2-1}\right)\right)}\right)$$

(8)

Proof of Lemma 3. 

Define

$${\kappa }_1=\frac{{\alpha }_k^2-1}{{\alpha }_k^2+1}$$

and

$${\lambda }_1=\frac{2{\alpha }_k}{{\alpha }_k^2+1}$$

such that

$${\beta }_k=\frac{2}{{\left({\kappa }_1+i{\lambda }_1\right)}^{2^{k-1}}-i}-i,$$

(9)

in accordance with Equation (7). Then, using de Moivre’s formula we can write the complex number in polar form as

$$\begin{array}{cc}\hfill {\left({\kappa }_1+i{\lambda }_1\right)}^{2^{k-1}}=& {\left({\kappa }_1^2+{\lambda }_1^2\right)}^{2^{k-2}}\left(cos\left(2^{k-1}\mathrm{Arg}\left({\kappa }_1+i{\lambda }_1\right)\right)\right)\hfill \\ & +isin\left(2^{k-1}\mathrm{Arg}\left({\kappa }_1+i{\lambda }_1\right)\right).\hfill \end{array}$$

Thus, substituting this equation into Equation (9), we obtain

$${\beta }_k=\frac{cos\left(2^{k-1}\mathrm{Arg}\left({\kappa }_1+i{\lambda }_1\right)\right)}{1-sin\left(2^{k-1}\mathrm{Arg}\left({\kappa }_1+i{\lambda }_1\right)\right)}=\frac{cos\left(2^{k-1}\mathrm{Arg}\left(\frac{{\alpha }_k+i}{{\alpha }_k-i}\right)\right)}{1-sin\left(2^{k-1}\mathrm{Arg}\left(\frac{{\alpha }_k+i}{{\alpha }_k-i}\right)\right)}.$$

(10)

Using the relation

$$\mathrm{Arg}\left(x+iy\right)=arctan\left(\frac{y}{x}\right),x>0,$$

we can write

$$\mathrm{Arg}\left({\kappa }_1+i{\lambda }_1\right)=arctan\left(\frac{2{\alpha }_k}{{\alpha }_k^2-1}\right).$$

(11)

Consequently, from the identities (10) and (11), we get

$${\beta }_k=\frac{cos\left(2^{k-1}arctan\left(\frac{2{\alpha }_k}{{\alpha }_k^2-1}\right)\right)}{1-sin\left(2^{k-1}arctan\left(\frac{2{\alpha }_k}{{\alpha }_k^2-1}\right)\right)}$$

(12)

and Equation (8) follows.  □

It should be noted that computation of the constant ${\beta }_k$ by using Equation (7) is not optimal. Specifically, at lager values of the integer k its application slows down the computation. Application of the Equation (12) for computation of the constant ${\beta }_k$ is also not desirable due to presence of the trigonometric functions. Theorem 2 and its proof show how this problem can be effectively resolved [7].

Theorem 2. 

The rational number ${\beta }_k$ is given by

$${\beta }_k=\frac{{\kappa }_k}{1-{\lambda }_k},$$

(13)

where the coefficients ${\kappa }_k$ and ${\lambda }_k$ can be found by a two-step iteration

$$\left\{\begin{array}{c}\hfill {\kappa }_n={\kappa }_{n-1}^2-{\lambda }_{n-1}^2\\ \hfill {\lambda }_n=2{\kappa }_{n-1}{\lambda }_{n-1},\end{array}\right.n=2,3,4,\dots ,k$$

(14)

with initial values

$${\kappa }_1=\frac{{\alpha }_k^2-1}{{\alpha }_k^2+1}$$

and

$${\lambda }_1=\frac{2{\alpha }_k}{{\alpha }_k^2+1}.$$

Proof of Theorem 2. 

We notice that the following power reduction

$$\begin{array}{cc}\hfill {\left({\kappa }_1+i{\lambda }_1\right)}^{2k-1}& =\stackrel{k-1\mathrm{powers}\mathrm{of}2}{\overbrace{{\left({\left({\left({\kappa }_1+i{\lambda }_1\right)}^2\right)}^{2\cdots }\right)}^2}}=\stackrel{k-2\mathrm{powers}\mathrm{of}2}{\overbrace{{\left({\left({\left({\kappa }_2+i{\lambda }_2\right)}^2\right)}^{2\cdots }\right)}^2}}\hfill \\ & =\stackrel{k-3\mathrm{powers}\mathrm{of}2}{\overbrace{{\left({\left({\left({\kappa }_3+i{\lambda }_3\right)}^2\right)}^{2\cdots }\right)}^2}}=\stackrel{k-n\mathrm{powers}\mathrm{of}2}{\overbrace{{\left({\left({\left({\kappa }_n+i{\lambda }_n\right)}^2\right)}^{2\cdots }\right)}^2}}\hfill \\ & ={\left({\left({\kappa }_{k-2}+i{\lambda }_{k-2}\right)}^2\right)}^2={\left({\kappa }_{k-1}+i{\lambda }_{k-1}\right)}^2={\kappa }_k+i{\lambda }_k,\hfill \end{array}$$

where the numbers ${\kappa }_n$ and ${\lambda }_n$ can be found by two-step iteration (14), leads to

$${\beta }_k=\frac{2}{{\kappa }_k-i{\lambda }_k-i}-i=\frac{2{\kappa }_k}{{\kappa }_k^2+{\left({\lambda }_k-1\right)}^2}+i\left(\frac{2\left(1-{\lambda }_k\right)}{{\kappa }_k^2+{\left({\lambda }_k-1\right)}^2}-1\right),$$

according to the Equation (10).

From the Lemma 2 it follows that the imaginary part of the value ${\beta }_k$ is equal to zero. Consequently, the equation above can be simplified as

$${\beta }_k=\frac{2{\kappa }_k}{{\kappa }_k^2+{\left({\lambda }_k-1\right)}^2}.$$

(15)

However, since

$$i\left(\frac{2\left(1-{\lambda }_k\right)}{{\kappa }_k^2+{\left({\lambda }_k-1\right)}^2}-1\right)=0$$

it follows that

$$\frac{2\left(1-{\lambda }_k\right)}{{\kappa }_k^2+{\left({\lambda }_k-1\right)}^2}=1\iff {\kappa }_k^2+{\left({\lambda }_k-1\right)}^2=2\left(1-{\lambda }_k\right).$$

Substituting this result into the denominator of Equation (15), we obtain Equation (13). This completes the proof.  □

There is an interesting relation between ${\beta }_k$ and nested radicals $a_{k-1}$ and $a_k$. Specifically, comparing Equation (4) with Equation (18) from the Lemma 4 below, one can see that

$$arctan\left(\frac{1}{{\beta }_k}\right)=-2^{k-1}arctan\left(\frac{\left\{a_k/\sqrt{2-a_{k-1}}\right\}}{1+⌊a_k/\sqrt{2-a_{k-1}}⌋\left(a_k/\sqrt{2-a_{k-1}}\right)}\right),$$

where

$$\left\{a_k/\sqrt{2-a_{k-1}}\right\}=a_k/\sqrt{2-a_{k-1}}-⌊a_k/\sqrt{2-a_{k-1}}⌋$$

denotes the fractional part.

Lemma 4. 

There is a limit such that

$$\underset{k\to \infty }{lim}arctan\left(\frac{1}{{\beta }_k}\right)=0.$$

Proof of Lemma 4 

It is not difficult to see that using change of the variable

$$y\to \frac{y}{1+\left(x+y\right)x}$$

in the trigonometric identity

$$arctan\left(x\right)+arctan\left(y\right)=arctan\left(\frac{x+y}{1-xy}\right)$$

leads to

$$arctan\left(x+y\right)=arctan\left(x\right)+arctan\left(\frac{y}{1+\left(x+y\right)x}\right).$$

(16)

Define for convenience the fractional part as

$${\gamma }_k=\left\{a_k/\sqrt{2-a_{k-1}}\right\}.$$

Then, from Equation (2) it follows that

$$\frac{\pi }{4}=2^{k-1}arctan\left(\frac{1}{{\alpha }_k+{\gamma }_k}\right)$$

Solving the equation

$$\frac{1}{{\alpha }_k+{\gamma }_k}=\frac{1}{{\alpha }_k}+z,$$

we get

$$z=-\frac{{\gamma }_k}{{\alpha }_k\left({\alpha }_k+{\gamma }_k\right)}.$$

Therefore, we have

$$\frac{\pi }{4}=2^{k-1}arctan\left(\frac{1}{{\alpha }_k}-\frac{{\gamma }_k}{{\alpha }_k\left({\alpha }_k+{\gamma }_k\right)}\right).$$

(17)

Using the identity (16), we can rewrite Equation (17) as

$$\frac{\pi }{4}=2^{k-1}\left(arctan\left(\frac{1}{{\alpha }_k}\right)-arctan\left(\frac{{\gamma }_k}{1+{\alpha }_k\left({\alpha }_k+\gamma \right)}\right)\right)$$

or

$$\begin{array}{cc}\hfill \frac{\pi }{4}& =2^{k-1}\left(arctan\left(\frac{1}{⌊a_k/\sqrt{2-a_{k-1}}⌋}\right)\right.\hfill \\ \hfill & \left.-arctan\left(\frac{\left\{a_k/\sqrt{2-a_{k-1}}\right\}}{1+⌊a_k/\sqrt{2-a_{k-1}}⌋\left(a_k/\sqrt{2-a_{k-1}}\right)}\right)\right).\hfill \end{array}$$

(18)

As the fractional part

$$\left\{a_k/\sqrt{2-a_{k-1}}\right\}<1$$

while

$$\underset{k\to \infty }{lim}\left(1+⌊a_k/\sqrt{2-a_{k-1}}⌋\left(a_k/\sqrt{2-a_{k-1}}\right)\right)=\infty ,$$

we conclude that

$$\begin{array}{cc}\hfill \underset{k\to \infty }{lim}& -2^{k-1}arctan\left(\frac{\left\{a_k/\sqrt{2-a_{k-1}}\right\}}{1+⌊a_k/\sqrt{2-a_{k-1}}⌋\left(a_k/\sqrt{2-a_{k-1}}\right)}\right)\hfill \\ & =\underset{k\to \infty }{lim}arctan\left(\frac{1-{\lambda }_k}{{\kappa }_k}\right)=\underset{k\to \infty }{lim}arctan\left(\frac{1}{{\beta }_k}\right)=0.\hfill \end{array}$$

 □

The Lemma 5 and its proof below show how to obtain a single-term rational approximation (33) for $\pi$ by truncating Equation (19).

Lemma 5. 

There is a limit such that

$$\frac{\pi }{4}=\underset{k\to \infty }{lim}\frac{2^{k-1}}{⌊a_k/\sqrt{2-a_{k-1}}⌋}.$$

(19)

Proof of Lemma 5. 

From the Lemma 4 it immediately follows that

$$\frac{\pi }{4}=\underset{k\to \infty }{lim}{2}^{k-1}arctan\left(\frac{1}{⌊a_k/\sqrt{2-a_{k-1}}⌋}\right).$$

(20)

According to Lemma 1

$$\underset{k\to \infty }{lim}\sqrt{2-a_{k-1}}=0$$

and

$$\underset{k\to \infty }{lim}{a}_k=2.$$

Therefore, we can infer that

$$\underset{k\to \infty }{lim}⌊a_k/\sqrt{2-a_{k-1}}⌋=\infty$$

or

$$\underset{k\to \infty }{lim}\frac{1}{⌊a_k/\sqrt{2-a_{k-1}}⌋}=0.$$

This limit implies that the argument of the arctangent function in Equation (20) tends to zero with increasing the integer k. Consequently, the limit (19) follows from the limit (20) and the proof is completed.  □

Motivated by recent publications [8,9,10,11,12,13] in connection to our works [5,7,14,15], we proposed a methodology for determination of the coefficients ${\alpha }_k$ without computing nested square roots of 2 (see Equation (5)) and developed an algorithm providing squared convergence per iteration in computing the digits of $\pi$. To the best of our knowledge, this approach is new and have never been reported.

2. Methodologies

2.1. Arctangent Function

There are different efficient methods to compute the arctangent function in the Machin-like formulas. We will consider a few of them.

In our previous publications [16] we show how the two-term Machin-like formula for $\pi$ represented in form (8) is used to derive an iterative formula

$${\theta }_{n,k}=\frac{1}{\frac{1}{{\theta }_{n-1,k}}+\frac{1}{2^k}\left(1-tan\left(\frac{2^{k-1}}{{\theta }_{n-1,k}}\right)\right)},k\ge 1,$$

(21)

where initial value can be taken as

$${\theta }_{1,k}=2^{-k}$$

such that

$$\frac{\pi }{4}=2^{k-1}\underset{n\to \infty }{lim}\frac{1}{{\theta }_{n,k}}.$$

(22)

This iterative formula provides a squared convergence in computing digits of $\pi$ (see Mathematica code in [16]).

Taking change of the variable ${\theta }_n\to 1/{\theta }_n$ in Equation (21) yields a more convenient form

$${\theta }_{n,k}={\theta }_{n-1,k}+2^{-k}\left(1-tan\left(2^{k-1}{\theta }_{n-1,k}\right)\right),k\ge 1.$$

Consequently, the limit (22) can be rearranged now as

$$\frac{\pi }{4}=2^{k-1}\underset{n\to \infty }{lim}{\theta }_{n,k}.$$

Comparing this limit with Equation (2), we can see that this iteration procedure results in

$$\underset{n\to \infty }{lim}{\theta }_{n,k}=arctan\left(\frac{1}{{\alpha }_k}\right)$$

and is used de facto for determination of the arctangent function. The detailed procedure showing how to implement the computation with high convergence rate is given in our recent publication [17].

Alternatively, we can transform the two-term into multi-term Machin-like formulas for $\pi$ consisting of only the integer reciprocals. In order to do this, we can use the following equation [14]

$$\frac{\pi }{4}=2^{k-1}arctan\left(\frac{1}{{\alpha }_k}\right)+\left(\sum _{m=1}^{M}arctan\left(\frac{1}{⌊{\mu }_{m,k}⌋}\right)\right)+arctan\left(\frac{1}{{\mu }_{M+1,k}}\right),$$

(23)

where

$${\mu }_{m,k}=\frac{1+⌊{\mu }_{m-1,k}⌋{\mu }_{m-1,k}}{⌊{\mu }_{m-1,k}⌋-{\mu }_{m-1,k}}$$

with an initial value

$${\mu }_{1,k}={\beta }_k.$$

For example, by taking $k=2$ in the Equation (23), we can construct

$$\begin{array}{cc}\hfill \frac{\pi }{4}& =2^{2-1}arctan\left(\frac{1}{{\alpha }_2}\right)+arctan\left(\frac{1}{{\mu }_{1,2}}\right)\hfill \\ \hfill & =2arctan\left(\frac{1}{2}\right)-arctan\left(\frac{1}{7}\right).\hfill \end{array}$$

(24)

This equation is commonly known as the Hermann’s formula for $\pi$ [18,19].

By taking $k=3$, the Equation (23) gives

$$\frac{\pi }{4}=2^{3-1}arctan\left(\frac{1}{{\alpha }_3}\right)-arctan\left(\frac{1}{{\mu }_{1,3}}\right)$$

representing the original Machin’s formula (1) for $\pi$.

The case $k=4$ requires more calculations to obtain the multi-term formula for $\pi$ consisting of only integer reciprocals. In particular, at $M=1$ we get

$$\begin{array}{cc}\hfill \frac{\pi }{4}& =2^{4-1}arctan\left(\frac{1}{{\alpha }_4}\right)+arctan\left(\frac{1}{{\mu }_{1,4}}\right)\hfill \\ & =8arctan\left(\frac{1}{10}\right)-arctan\left(\frac{1758719}{147153121}\right)\hfill \end{array}$$

At $M=2$ and $M=3$ Equation (23) yields

$$\begin{array}{cc}\hfill \frac{\pi }{4}& =2^{4-1}arctan\left(\frac{1}{{\alpha }_4}\right)+arctan\left(\frac{1}{⌊{\mu }_{1,4}⌋}\right)+arctan\left(\frac{1}{{\mu }_{2,4}}\right)\hfill \\ & =8arctan\left(\frac{1}{10}\right)-arctan\left(\frac{1}{84}\right)-arctan\left(\frac{579275}{12362620883}\right)\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill \frac{\pi }{4}& =2^{4-1}arctan\left(\frac{1}{{\alpha }_4}\right)+arctan\left(\frac{1}{⌊{\mu }_{1,4}⌋}\right)+arctan\left(\frac{1}{⌊{\mu }_{2,4}⌋}\right)+arctan\left(\frac{1}{{\mu }_{3,4}}\right)\hfill \\ & =8arctan\left(\frac{1}{10}\right)-arctan\left(\frac{1}{84}\right)-arctan\left(\frac{1}{21342}\right)\hfill \\ & -arctan\left(\frac{266167}{263843055464261}\right),\hfill \end{array}$$

respectively.

Repeating this procedure over and over again, at $M=5$ we end up with 7-term Machin-like formula for $\pi$ consisting of only integer reciprocals

$$\begin{array}{cc}\hfill \frac{\pi }{4}& =2^{4-1}arctan\left(\frac{1}{{\alpha }_4}\right)+\left(\sum _{m=1}^{5}arctan\left(\frac{1}{⌊{\mu }_{m,4}⌋}\right)\right)+arctan\left(\frac{1}{{\mu }_{6,41}}\right)\hfill \\ \hfill & =8arctan\left(\frac{1}{10}\right)-arctan\left(\frac{1}{84}\right)-arctan\left(\frac{1}{21342}\right)\hfill \\ \hfill & -arctan\left(\frac{1}{991268848}\right)-arctan\left(\frac{1}{193018008592515208050}\right)\hfill \\ \hfill & -arctan\left(\frac{1}{197967899896401851763240424238758988350338}\right)\hfill \\ \hfill & -arctan\left(\frac{1}{\mathsf{\Omega }}\right),\hfill \end{array}$$

(25)

where

$$\begin{array}{cc}\hfill \mathsf{\Omega }=& 11757386816817535293027775284419412676799191500853701\dots \hfill \\ & 8836932014293678271636885792397,\hfill \end{array}$$

is the largest integer (see the Mathematica code in [17] that validates this seven-term Machin-like formula for $\pi$).

As we can see, Hermann’s (24), Machin’s (1) and derived (25) formulas for $\pi$ belong to the same generic group as all of them can be constructed from the same equation-template (23).

Our empirical results show that two arctangent series expansions can be used for computation with rapid convergence. The first equation is Euler’s expansion series given by [20]

$$arctan\left(x\right)=\sum _{n=0}^{\infty }\frac{2^{2n}{\left(n!\right)}^2}{\left(2n+1\right)!}\frac{x^{2n+1}}{{\left(1+x^2\right)}^{n+1}}.$$

(26)

The second equation is [7,21]

$$arctan\left(x\right)=2\sum _{n=1}^{\infty }\frac{1}{2n-1}\frac{g_n\left(x\right)}{g_n^2\left(x\right)+h_n^2\left(x\right)},$$

(27)

where

$$g_n\left(x\right)=g_{n-1}\left(x\right)\left(1-4/x^2\right)+4h_{n-1}\left(x\right)/x$$

and

$$h_n\left(x\right)=h_{n-1}\left(x\right)\left(1-4/x^2\right)-4g_{n-1}\left(x\right)/x$$

with initial values

$$g_1\left(x\right)=2/x$$

and

$$h_1\left(x\right)=1.$$

Computational test we performed shows that Equation (27) is faster in convergence than Equation (26). Recently, we proposed the generalization of the arctangent series expansion (27) [21].

2.2. Tangent Function

Generally, transformation of the two-terms to multi-terms formulas for $\pi$ with integer reciprocals is not required. In particular, we can use Newton–Raphson iteration method [15]. For example, both arctangent terms in the two-term Machin-like formula (4) for $\pi$ can be computed directly by using the following iterative formulas

$${\sigma }_n={\sigma }_{n-1}-\left(1-{\left(\frac{2tan\left({\sigma }_{n-1}/2\right)}{1+{tan}^2\left({\sigma }_{n-1}/2\right)}\right)}^2\right)\left(\frac{2tan\left({\sigma }_{n-1}/2\right)}{1-{tan}^2\left({\sigma }_{n-1}/2\right)}-\frac{1}{{\alpha }_k}\right)$$

and

$${\tau }_n={\tau }_{n-1}-\left(1-{\left(\frac{2tan\left({\tau }_{n-1}/2\right)}{1+{tan}^2\left({\tau }_{n-1}/2\right)}\right)}^2\right)\left(\frac{2tan\left({\tau }_{n-1}/2\right)}{1-{tan}^2\left({\tau }_{n-1}/2\right)}-\frac{1}{{\beta }_k}\right)$$

with initial values

$${\sigma }_1=\frac{1}{{\alpha }_k}$$

and

$${\tau }_1=\frac{1}{{\beta }_k}$$

such that

$$arctan\left(\frac{1}{{\alpha }_k}\right)=\underset{n\to \infty }{lim}{\sigma }_n$$

and

$$arctan\left(\frac{1}{{\beta }_k}\right)=\underset{n\to \infty }{lim}{\tau }_n.$$

Since the convergence of the Newton–Raphson iteration is quadratic [25], with proper implementation of the tangent function we may achieve an efficient computation.

The tangent function can be expanded as

$$\begin{array}{cc}\hfill tan\left(x\right)& =\sum _{n=1}^{\infty }\frac{{\left(-1\right)}^{n-1}{2}^{2n}\left(2^{2n}-1\right)B_{2n}}{\left(2n\right)!}{x}^{2n-1}\hfill \\ \hfill & =x+\frac{x^3}{3}+\frac{2x^5}{15}+\frac{17x^7}{315}+\frac{62x^9}{2835}+\cdots ,\hfill \end{array}$$

(28)

where $B_{2n}$ are the Bernoulli numbers, defined by the contour integral

$$B_n=\frac{n!}{2\pi i}\oint \frac{z}{e^z-1}\frac{dz}{z^{n+1}}.$$

However, application of Equation (28) is not desirable since the computation of the Bernoulli numbers $B_{2n}$ itself is a big challenge [22,23,24,26].

In order to resolve this problem, we proposed the following limit [17]

$$tan\left(x\right)=\underset{n\to \infty }{lim}\frac{2p_n^2\left(x\right)}{q_n\left(x\right)},$$

(29)

where

$$p_n\left(x\right)=p_{n-1}\left(x\right)+r_{n-1}\left(x\right),$$

$$q_n\left(x\right)=q_{n-1}\left(x\right)+2^{2n-1}{r}_{n-1}\left(x\right),$$

such that

$$r_n=\frac{{\left(-1\right)}^n}{\left(2n+1\right)!}{x}^{2n+1},$$

with initial values

$$p_0\left(x\right)=0,$$

$$q_0\left(x\right)=0.$$

Specifically, it has been shown that at $k=27$ application of Equation (29) results in more than 17 digits of $\pi$ per increment n (see the Mathematica codes in [17]).

3. Algorithmic Implementation

In our recent publication we have shown that [16]

$${\beta }_k\approx \frac{2}{1-tan\left(\frac{2^{k-1}}{{\alpha }_k}\right)}.$$

Consequently, in accordance with Equation (4), we obtain the following approximation

$$\frac{\pi }{4}\approx 2^{k-1}arctan\left(\frac{1}{{\alpha }_k}\right)+\frac{1}{2}\left(1-tan\left(\frac{2^{k-1}}{{\alpha }_k}\right)\right).$$

At sufficiently large k the value $1/{\alpha }_k<<1$. Therefore, according to the the Lemma 5, in this equation we can replace the first arctangent term by a rational number $2^{k-1}/{\alpha }_k$. This gives

$$\frac{\pi }{4}\approx \frac{2^{k-1}}{{\alpha }_k}+\frac{1}{2}\left(1-tan\left(\frac{2^{k-1}}{{\alpha }_k}\right)\right).$$

Unfortunately, we cannot compute efficiently the tangent function in this approximation since its argument $2^{k-1}/{\alpha }_k$ is not a small number as it tends to $\pi /4$ with increasing k. However, again by taking into consideration that $1/{\alpha }_k<<1$ from which it follows that

$$\frac{1}{{\alpha }_k}\approx arctan\left(\frac{1}{{\alpha }_k}\right),$$

we can write

$$\frac{\pi }{4}\approx \frac{2^{k-1}}{{\alpha }_k}+\frac{1}{2}\left(1-tan\left(2^{k-1}arctan\left(\frac{1}{{\alpha }_k}\right)\right)\right).$$

Now we can take advantage from the fact that the multiplier $2^{k-1}$ is continuously divisible by 2. Therefore, we can use the trigonometric identity

$$tan\left(2arctan\left(x\right)\right)=\frac{2x}{1-x^2}$$

$k-1$ times over and over again. Thus, this leads to the following iterative formula

$${\eta }_n\left(x\right)=\frac{2{\eta }_{n-1}\left(x\right)}{1-{\eta }_{n-1}^2\left(x\right)}$$

(30)

with an initial value

$${\eta }_1\left(x\right)=\frac{2x}{1-x^2}$$

such that

$${\eta }_{k-1}\left(\frac{1}{{\alpha }_k}\right)=tan\left(2^{k-1}arctan\left(\frac{1}{{\alpha }_k}\right)\right).$$

Since the left side of the equation above provides an exact value without (tangent and arctangent) trigonometric functions, we can regard this equation

$$\frac{\pi }{4}\approx \frac{2^{k-1}}{{\alpha }_k}+\frac{1}{2}\left(1-{\eta }_{k-1}\left(\frac{1}{{\alpha }_k}\right)\right)$$

(31)

as a rational approximation.

This rational approximation of the two-term Machin-like formula for $\pi$ can be used in an algorithm providing a quadratic convergence. This can be achieved with help of the Theorem 3.

Theorem 3. 

There is a formula for $1/\pi$ with binary output

$$\frac{1}{\pi }=\underset{k\to \infty }{lim}\sum _{n=1}^k\frac{1}{{10}^{n+1}}\left({\alpha }_{n}mod2\right)={\left[0.01010001011111001100\dots \right]}_2,$$

(32)

Proof of Theorem 3. 

The proof is related to the parity of the integer ${\alpha }_k$. According to the Lemma 5, we can write

$$\frac{\pi }{4}\approx \frac{2^{k-1}}{{\alpha }_k}$$

(33)

or

$$\frac{1}{\pi }\approx \frac{{\alpha }_k}{4\times 2^{k-1}}.$$

Consequently, if the integer ${\alpha }_k$ is even, then

$$\frac{{\alpha }_k}{4\times 2^{k-1}}-\frac{{\alpha }_{k-1}}{4\times 2^{k-2}}=0_2.$$

However, if the integer ${\alpha }_k$ is odd, then

$$\frac{{\alpha }_k}{4\times 2^{k-1}}-\frac{{\alpha }_{k-1}}{4\times 2^{k-2}}=\underset{k+1\mathrm{zeros}}{\underbrace{0.00\dots 00}}{1}_2.$$

This means that ${\alpha }_k$ contributes a binary digit $\underset{k+1\mathrm{zeros}}{\underbrace{0.00\dots 00}}{1}_2$ to the previous value ${\alpha }_{k-1}$ if and only if ${\alpha }_k$ is odd. This completes the proof.    □

Consider an example. There are four consecutive values ${\alpha }_4=10$, ${\alpha }_5=20$ and ${\alpha }_6=40$ and ${\alpha }_7=81$. Since the first three values are even, we have

$$\frac{{\alpha }_4}{4\times 2^{4-1}}=\frac{{\alpha }_5}{4\times 2^{5-1}}=\frac{{\alpha }_6}{4\times 2^{6-1}}=0.{0101}_2.$$

However, since ${\alpha }_7=81$ is odd, we obtain

$$\frac{{\alpha }_7}{4\times 2^{7-1}}=0.{0101}_2+0.{00000001}_2=0.{01010001}_2.$$

Consider how number of digits of $\pi$ can be doubled without computing square roots for the nested radicals $a_k$. We can take, for example, $k=7$. This yields

$$\begin{array}{cc}\hfill {\alpha }_7& =⌊\frac{a_7}{\sqrt{2-a_{7-1}}}⌋\hfill \\ \hfill & =⌊\frac{\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2}}}}}}}}{\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2}}}}}}}}⌋=81.\hfill \end{array}$$

(34)

However, it is not reasonable to compute the square roots of 2 so many times to obtain this number. Instead, we can simplify computation considerably by using the value of the $1/\pi$ in the binary form according to the Theorem 3. Thus, ignoring the first two initial zeros in the binary output of the Equation (32), we have a corresponding sequence

$${\left(s\right)}_k^{\infty }=\left(1,0,1,0,0,0,1,0,1,1,1,1,1,0,0,\dots \right).$$

(35)

This sequence can be obtained by using the built-in Mathematica directly using $1/\pi$. For example, the following code:

 

RealDigits[

  ImportString[ToString[BaseForm[N[1/\[Pi],20],2]],

    "Table"][[1]][[1]]][[1]][[1;;20]]

 

returns the first 20 digits from the sequence (35):

 

{1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0}

 

From this sequence, we choose the sub-sequence (say, up to seventh element)

$${\left(s\right)}_k^7=\left(1,0,1,0,0,0,1\right)$$

and apply it accordingly as

$${\alpha }_k=\left\{\begin{array}{ccc}2{\alpha }_{k-1},\hfill & \hfill & \mathrm{if}{\mathit{k}}^{\mathrm{th}}\mathrm{binary}\mathrm{digit}\mathrm{is}0\hfill \\ 2{\alpha }_{k-1}+1,\hfill & \hfill & \mathrm{if}{\mathit{k}}^{\mathrm{th}}\mathrm{binary}\mathrm{digit}\mathrm{is}1.\hfill \end{array}\right.$$

Explicitly, this step-by-step procedure results in

$$\begin{array}{cc}\hfill {\alpha }_1& =2{\alpha }_0+1=\left(2\times 0+1\right)=1,\hfill \\ \hfill {\alpha }_2& =2{\alpha }_1+0=\left(2\times 1+0\right)=2,\hfill \\ \hfill {\alpha }_3& =2{\alpha }_2+1=\left(2\times 2+1\right)=5,\hfill \\ \hfill {\alpha }_4& =2{\alpha }_3+0=\left(2\times 5+0\right)=10,\hfill \\ \hfill {\alpha }_5& =2{\alpha }_4+0=\left(2\times 10+0\right)=20,\hfill \\ \hfill {\alpha }_6& =2{\alpha }_5+0=\left(2\times 20+0\right)=40,\hfill \\ \hfill {\alpha }_7& =2{\alpha }_6+1=\left(2\times 40+1\right)=81.\hfill \end{array}$$

Thus, we can see how a very simple procedure can be used to determine the value of the rational number ${\alpha }_7=81$ without using a sophisticated Equation (34) consisting of 14 nested square roots of 2.

At ${\alpha }_7=81$ the corresponding Machin-like formula is

$$\begin{array}{cc}\hfill \frac{\pi }{4}& =2^{7-1}arctan\left(\frac{1}{{\alpha }_7}\right)+arctan\left(\frac{1}{{\beta }_7}\right)\hfill \\ \hfill & =64arctan\left(\frac{1}{81}\right)-arctan\left(\frac{\stackrel{111\mathrm{digits}}{\overbrace{2154947582\dots 4298183679}}}{\underset{113\mathrm{digits}}{\underbrace{4599489202\dots 6981324801}}}\right),\hfill \end{array}$$

(36)

where the constant

$${\beta }_7=-\frac{\stackrel{113\mathrm{digits}}{\overbrace{4599489202\dots 6981324801}}}{\underset{111\mathrm{digits}}{\underbrace{2154947582\dots 4298183679}}}$$

can be computed either by using Equation (7) or, more efficiently, by using Equation (13) based on two-step iteration (14).

The following Mathematica code:

 

(* String for long number \[Beta]_7 *)

strBeta7=

  ToString[StringJoin[

    "21549475820057881611210311984288158234143531212163819254",

      "1568712000964806160594022446140062110943660584298183679/",

        "459948920218008069525744651226752553899687099736076594466",

          "78719072620659988130828378620624183170066256006981324801"

            ]];

 

(* Verification *)

Print[\[Pi]/4==64*ArcTan[1/81]-ArcTan[ToExpression[strBeta7]]];

 

validates the Equation (36) by returning True.

Suppose that we do not know the sequence other than ${\left(s\right)}_1^7$. However, with help of the Equation (31) we can find other digits of $\pi$ in the iterative process. In particular, using Equation (30), we have

$${\eta }_{7-1}\left(\frac{1}{81}\right)=\frac{\stackrel{113\mathrm{digits}}{\overbrace{2310519339\dots 5639754240}}}{\underset{113\mathrm{digits}}{\underbrace{2288969863\dots 1341570561}}}\approx 1.00941448647564092749$$

Substituting this value into Equation (31), we can find a significantly better approximation of $\pi$.

The following is the Mathematica code:

 

Print["Equation (33) at k = 7: ",

  MantissaExponent[N[\[Pi]-4*(64/81),20]][[2]]//Abs,

    " digits of \[Pi]"];

Print["Equation (31) at k = 7: ",

  MantissaExponent[

    N[\[Pi]-4*(64/81+1/2*(1-1.00941448647564092749)),

      20]][[2]]//Abs," digits of \[Pi]"];

 

produces the following output:

 

Equation~(33) at k = 7: 1 digits of À

Equation~(31) at k = 7: 4 digits of À

 

Initial sequence ${\left(s\right)}_1^7$ helped us to find the value ${\alpha }_7=81$. Now due to higher accuracy of Equation (31) we can generate the sequence in which its upper index is doubled

$${\left(s\right)}_1^{14}=\left(1,0,1,0,0,0,1,0,1,1,1,1,1,0\right)$$

and with help of this sequence we can find the corresponding value ${\alpha }_{14}=10430$.

Unfortunately, doubling the upper index k does not always work. For example, if we attempt to double the upper index by using the initial sequence ${\left(s\right)}_1^8$, then we get ${\alpha }_{16}=41722$ instead of correct value ${\alpha }_{16}=41721$. Therefore, the upper index of the sequence should be slightly less than two.

The two-terms approximation (31) doubles the number of digits of $\pi$ as compared to the single-term approximation (33). This means that using the sequence ${\left(s\right)}_1^k$ we can obtain all sequences ${\left(s\right)}_1^{k+1}$, ${\left(s\right)}_1^{k+2}$, ${\left(s\right)}_1^{k+3}$, etc. up to ${\left(s\right)}_1^{k_0}$, where $k_0$ is an integer slightly smaller than $2k$. Our numerical results show that doubling value k does not always provide the correct sequence as a few binary digits at the end of the sequence ${\left(s\right)}_1^{2k}$ occasionally may not be correct. However, when we use the empirical equation

$$k_0=⌊\left(2-\frac{1}{32}\right)k⌋,$$

then the corresponding sequence ${\left(s\right)}_1^{k_0}$ is a sub-sequence of the infinite sequence (35) and, therefore, it is appeared to be correct. It is interesting to note that the number 32 in this equation is the largest integer that we found on the basis of our numerical results.

The following is a Mathematica code that shows number of digits of $\pi$ at given iteration number n and integer k:

 

Clear[str,sps,k,\[Gamma],\[Alpha],lst,\[Eta]]

 

(* String for conversion of 1/\[Pi] to sequence *)

str="ImportString[ToString[BaseForm[N[1/piAppr,k0],2]],

  \"Table\"][[1]][[1]]";

 

(* String for space separation *)

sps[n_]:=Module[{m=1,sps=" "},

  While[m<n,sps=StringJoin[sps," "];m++];If[m==n,sps]];

 

(* Converting number to string with length q *)

cnv2str[p_,q_]:=Module[{},StringTake[StringJoin[ToString[p],sps[q]],q]]

 

(* Defining \[Eta]-function *)

\[Eta][n_,x_,k_]:=Module[{K=k/1.5,y=x},y=N[(2*y)/(1-y^2),K];cntr = 1;

  While[cntr<n,y=(2*y)/(1-y^2);cntr++];y];

 

(* Define \[Alpha]_1, \[Alpha]_2 and \[Alpha]_3] *)

\[Alpha][1]=1;

\[Alpha][2]=2;

\[Alpha][3]=5;

 

(* Input values *)

k=3;\[Gamma]=\[Alpha][3];

 

(* Heading *)

Print["-------------------------------"];

Print["Iteration | k", sps[5], "| Digits of \[Pi]"];

Print["-------------------------------"];

 

n=1;

While[n<=12,

 

  intR=1/\[Gamma];

  k0=\[LeftFloor](2-1/32)*k\[RightFloor];

 

  piAppr=4*(2^(k-1)*intR+1/2*(1-\[Eta][k-1,intR,k0]));

 

  (* Extracting the sequence {1,0,1,0,0,0,1...} *)

  lst=RealDigits[ToExpression[str]][[1]][[1;;k0]];

 

  (* Main computation *)

  K=k+1;

  While[K<=k0,\[Gamma]=

    2*\[Gamma]+lst[[K]];\[Alpha][K]=\[Gamma];K++];k=k0;

 

  (* Aligned output" *)

  Print[cnv2str[n,5],sps[4]," | ",cnv2str[k,5]," | ",

    MantissaExponent[N[\[Pi]-piAppr,k0]][[2]]//Abs];

 

  n++];

 

This code generates the output:

 

------------------------------

Iteration | k    | Digits of À

------------------------------

1         | 5    | 1

2         | 9    | 2

3         | 17   | 4

4         | 33   | 9

5         | 64   | 20

6         | 126  | 38

7         | 248  | 75

8         | 488  | 149

9         | 960  | 293

10        | 1890 | 577

11        | 3720 | 1137

12        | 7323 | 2240

 

As we can see from the third column, the number of digits of $\pi$ doubles at each iteration.

The Mathematica code below shows the values of ${\alpha }_k$ at given k:

 

(*Heading*)

Print["-----------------"];

Print["k",sps[5],"| \[Alpha][k]"];

Print["-----------------"];

 

k=2;

While[k<=25,

 

  (*Aligned output" *)

  Print[cnv2str[k,5], " | ", \[Alpha][k]];

  k++];

 

This code returns the following output:

 

—————–

k     | α${}_k$

—————–

2     | 2

3     | 5

4     | 10

5     | 20

6     | 40

7     | 81

8     | 162

9     | 325

10    | 651

11    | 1303

12    | 2607

13    | 5215

14    | 10430

15    | 20860

16    | 41721

17    | 83443

18    | 166886

19    | 333772

20    | 667544

21    | 1335088

22    | 2670176

23    | 5340353

24    | 10680707

25    | 21361414

 

As we can see, all numbers ${\alpha }_k$ are the same as those reported in [9].

Since the constants ${\alpha }_k$ have been computed already, we can use them to validate the formula (32) for $1/\pi$ in the binary form. The Mathematica code below:

 

f[K_]:=N[Sum[(1/10^(k+1))*Mod[\[Alpha][k],2],{k,1,K}],K];

Print["-------------------------------------------------------------"]

Print[k,sps[2]," | Binary output"];

Print["-------------------------------------------------------------"]

 

k := 1;

While[k<=5,Print[10*k,sps[2],"| ",Subscript[f[10*k],2]];k++]

 

Print["-------------------------------------------------------------"]

 

Print["Built-in Mathematica:"]

Print["1/\[Pi]=",

  Subscript[StringJoin["[",

    StringSplit[ToString[BaseForm[N[1/Pi,15],2]]][[1]], "...]"],

    2]];

 

returns the output:

 

————————————————————-

k  | Binary output

————————————————————-

10 | 0.0001000101100${}_2$

20 | 0.00010001011111001100000${}_2$

30 | 0.000100010111110011000001101101100${}_2$

40 | 0.0001000101111100110000011011011100100111000${}_2$

50 | 0.00010001011111001100000110110111001001110010001000000${}_2$

————————————————————-

Built-in Mathematica:

$1/\pi$ = [0.010100010111110011000001101101110010011100100010000...]${}_2$

 

according to the Equation (32). The original binary representation of the number $1/\pi$ generated by built-in Mathematica is also shown for comparison (see also [27] for binary sequence for $1/\pi$).

4. Conclusions

We consider the properties of the two-term Machin-like formula for $\pi$ and propose its two-term rational approximation (31). Using this approach, we developed an efficient algorithm for computing digits of $\pi$ with squared convergence. The constants ${\alpha }_k$ in this approximation are computed without nested square roots of 2.