Application of a New Iterative Formula for Computing π and Nested Radicals with Roots of 2

1. Introduction

Throughout many centuries computing digits of $\pi$ remained a big challenge [1,2,3,4]. However, in 1876 English astronomer and mathematician John Machin found an efficient method to resolve this problem. Historically, he was the first to calculate over 100 digits of $\pi$. In his approach, John Machin discovered and then used the following remarkable formula [1,2,3,4]

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

(1)

Nowadays, the identities of kind

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

(2)

are named after him as the Machin-like formulas for $\pi$.

The arctangent terms in the Machin-like formulas for $\pi$ can be computed by using the Maclaurin expansion series

$$arctan\left(x\right)=x-{\displaystyle \frac{x^3}{3!}}+{\displaystyle \frac{x^5}{5!}}-{\displaystyle \frac{x^7}{7!}}+\dots =\sum _{n=0}^{\infty }{\displaystyle \frac{x^{2k-1}}{(2k-1)!}},x\le \left|1\right|.$$

(3)

Since according to this equation, we get

$$arctan\left(x\right)=x+O\left(x^3\right),$$

the convergence of the Machin-like formulas (2) for $\pi$ is always better when the values $B_j$ are larger by absolute value.

To estimate efficiency of the Machin-like formulas for $\pi$, Lehmer introduced a measure, defined as [5,6,7]

$$\mu =\sum _{j=1}^J{\displaystyle \frac{1}{{log}_{10}\left(\left|B_j\right|\right)}}.$$

According to this formula, the smaller measure indicates the higher efficiency of the Machin-like formula $\pi$. Lehmer’s measure is smaller at larger absolute values of $B_j$ and smaller number of the summation terms J.

Generally, we should imply that Lehmer’s measure is valid only if all coefficients $B_j$ are integers. Otherwise if $B_j\notin \mathbb{Z}$, its fractional parts

$$\left\{B_j\right\}=B_j-\lfloor B_j\rfloor$$

may cause further computational complexities requiring more usage of the computer memory and extending considerably a run-time in computing digits of $\pi$[7]. As the fractional part is not desirable in computation of the digits of $\pi$, it may be more preferable to apply the Machin-like formulas where all coefficients $B_j$ are integers.

The Machin-like formulas (2) for $\pi$ can be validated by using the following relation

$$\prod _{j=1}^J{\left(B_j+i\right)}^{A_j}\propto (1+i).$$

(4)

The right side of this relation implies that the real part of the product must be equal to its imaginary part as follows

$$\Re \left\{\prod _{j=1}^J{\left(B_j+i\right)}^{A_j}\right\}=\Im \left\{\prod _{j=1}^J{\left(B_j+i\right)}^{A_j}\right\}.$$

For example, the original Machin formula (1) for $\pi$ can be readily validated by applying this relation

$${(5+i)}^4{(259+i)}^{-1}=2(1+i)$$

since the real and imaginary parts of the product are equal to the same number 2.

In our previous publication [8], we proposed a method for deriving the two-term Machin-like formula for $\pi$ in form

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

(5)

where $\alpha$ is an integer (see [9] for available values of $\alpha$)

$$\alpha =⌊{\displaystyle \frac{c_k}{\sqrt{2-c_{k-1}}}}⌋$$

that can be computed by using the nested radicals such that

$$c_0=0,c_k=\sqrt{2+c_{k-1}}$$

and

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

(6)

Recently, Gasull et al. proposed an alternative method to obtain the two-term Machin-like formula for $\pi$ of kind (5). In their publication, they suggested an iterative method based on $R_j(n,x)$ functions that can be defined as [7]

$$R_j(n,x)=tan(n\sigma +j\pi /4),x=tan\left(\sigma \right),n,j\in \mathbb{N}.$$

Application of formula (6) for determination of the coefficient $\beta$ is not optimal. In particular, as integer k increases the exponent $2^{k-1}$ in the denominator blows up very rapidly. As a consequence, this drastically delays computation. To resolve this problem, we proposed a two-step iteration procedure [8] that will be discussed in the next section. However, this iteration involves squaring that doubles the number of digits at each consecutive cycle of two-step iteration. Although as compared to equation (6) this iterative procedure is more efficient, it still remains computationally costly at larger values of k.

Motivated by recent publications [10,11,12,13,14,15,16], we propose a more efficient approach in determination of the coefficient $\beta$ in equation (6). Furthermore, with these results we show how the constant $\pi$ and nested radicals consisting of square roots of 2 can be obtained. Some numerical results with Mathematica codes are provided. To the best of our knowledge, this approach is new and has never been reported in scientific literature.

2. Preliminaries

As it has been mentioned above, equation (6) should be avoided for computation of the constant $\beta$ at $k\gg 1$. The following theorem shows how else the constant $\beta$ can be calculated.

Theorem 1. 

There is a formula such that

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

(7)

where

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

(8)

with initial values

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

and

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

Proof. 

From induction it follows that

$$\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\dots }\right)}^2}}=\stackrel{k-2\mathrm{powers}\mathrm{of}2}{\overbrace{{\left({\left({\left({\kappa }_2+i{\lambda }_2\right)}^2\right)}^{2\dots }\right)}^2}}\hfill \\ \hfill & =\stackrel{k-3\mathrm{powers}\mathrm{of}2}{\overbrace{{\left({\left({\left({\kappa }_3+i{\lambda }_3\right)}^2\right)}^{2\dots }\right)}^2}}=\stackrel{k-n\mathrm{powers}\mathrm{of}2}{\overbrace{{\left({\left({\left({\kappa }_n+i{\lambda }_n\right)}^2\right)}^{2\dots }\right)}^2}}\hfill \\ \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}$$

(9)

Therefore, from equations (6) and (9), we obtain

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

(10)

Since $\beta$ is a real number, the imaginary part of the equation above must be equal to zero. Consequently, we have

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

or

$$2(1-{\lambda }_k)={\kappa }_k^2+{({\lambda }_k-1)}^2$$

or

$${\kappa }_k^2=2(1-{\lambda }_k)-{({\lambda }_k-1)}^2.$$

(11)

Substituting equation (11) into the real part of equation (10), we obtain equation (7). This completes the proof.    □

In our earlier publication [17], we show how to generate efficiently the multi-term Machin-like formulas for $\pi$ by using the following equation-template

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

(12)

where

$${\theta }_{m,k}={\displaystyle \frac{1+⌊{\theta }_{m-1,k}⌋{\theta }_{m-1,k}}{⌊{\theta }_{m-1,k}⌋-{\theta }_{m-1,k}}},m\ge 2.$$

In the formula (12), we imply that $M\ge 0$ such that at $M=0$ the sum

$$\sum _{m=1}^{M}arctan\left({\displaystyle \frac{1}{⌊{\theta }_{m,k}⌋}}\right)$$

is equal to zero. Consequently, when $M=0$ the equation (12) is reduced to

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

(13)

where ${\gamma }_k=\alpha$ and ${\theta }_{1,k}=\beta$ according to equation (5) above. In fact, the equation (12) can be obtained from equation (5) together with the following identity [17]

$$arctan\left({\displaystyle \frac{1}{z}}\right)=arctan\left({\displaystyle \frac{1}{⌊z⌋}}\right)+arctan\left({\displaystyle \frac{⌊z⌋-z}{1+z⌊z⌋}}\right),z\notin [0,1).$$

Consider how the two well-known Machin-like formulas for $\pi$ can be derived by using equation (12). At $k=2$ and $M=0$ equation (12) results in

$${\displaystyle \frac{\pi }{4}}=2arctan\left({\displaystyle \frac{1}{2}}\right)-arctan\left({\displaystyle \frac{1}{7}}\right).$$

(14)

This equation is commonly known as the Hermann’s formula for $\pi$[18]. At $k=2$ and $M=0$ equation (12) leads to the original Machin formula (1) for $\pi$.

Application of the equation (12) may also be a useful technique to transform the arctangent term with the quotient ${\theta }_{1,k}=\beta$ into sum of arctangents with reciprocal integers. For example, at $k=4$ we get

$${\gamma }_4=⌊{\displaystyle \frac{c_4}{\sqrt{2-c_3}}}⌋=⌊{\displaystyle \frac{\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2}}}}}{\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2}}}}}}⌋=10.$$

Consequently, applying two-step iteration (8) we can find that

$${\theta }_{1,4}=-{\displaystyle \frac{147153121}{1758719}}.$$

Substituting these values into equation (12), results in

$${\displaystyle \frac{\pi }{4}}=8arctan\left({\displaystyle \frac{1}{10}}\right)-arctan\left({\displaystyle \frac{1758719}{147153121}}\right).$$

(15)

As it has been mentioned above, it is desirable to apply reciprocal integers rather than quotients. The equation (12) can be used to transform the quotients into reciprocal integers. For example, at $k=4$ and $M=1$ from equation (12) it follows that

$${\displaystyle \frac{\pi }{4}}=8arctan\left({\displaystyle \frac{1}{10}}\right)-arctan\left({\displaystyle \frac{1}{84}}\right)-arctan\left({\displaystyle \frac{579275}{12362620883}}\right).$$

At $k=4$ and $M=2$ equation (12) gives

$$\begin{array}{cc}\hfill {\displaystyle \frac{\pi }{4}}& =8arctan\left({\displaystyle \frac{1}{10}}\right)-arctan\left({\displaystyle \frac{1}{84}}\right)-arctan\left({\displaystyle \frac{1}{21342}}\right)\hfill \\ & -arctan\left({\displaystyle \frac{266167}{263843055464261}}\right).\hfill \end{array}$$

Incrementing the integer M over and over again, at $k=4$ and $M=5$ equation (12) finally yields the 7-term Machin-like formula

$$\begin{array}{cc}\hfill {\displaystyle \frac{\pi }{4}}& =8arctan\left({\displaystyle \frac{1}{10}}\right)-arctan\left({\displaystyle \frac{1}{84}}\right)-arctan\left({\displaystyle \frac{1}{21342}}\right)\hfill \\ \hfill & -arctan\left({\displaystyle \frac{1}{991268848}}\right)-arctan\left({\displaystyle \frac{1}{193018008592515208050}}\right)\hfill \\ \hfill & -arctan\left({\displaystyle \frac{1}{197967899896401851763240424238758988350338}}\right)\hfill \\ \hfill & -arctan\left({\displaystyle \frac{1}{\Omega }}\right),\hfill \end{array}$$

(16)

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

where all arctangent arguments are reciprocal integers.

As we can see, the Hermann’s (14), Machin’s (1) and the generated (16) formulas belong to the same generic group as all of them can be derived by using the same equation-template (12).

The following Mathematica code validates the equation (16).

(* Define long string *)

longStr = StringJoin["11757386816817535293027775284419412676",

    "7991915008537018836932014293678271636885792397"];

(* Computing the coefficient *)

coeff = (10 + I)^8*(84 + I)^-1*(21342 + I)^-1*(991268848 + I)^-1*

    (193018008592515208050+I)^-1*

    (197967899896401851763240424238758988350338 + I)^-1*

            (FromDigits[longStr]+I)^-1;

Print[Re[coeff] == Im[coeff]];

by returning the output True.

Once all quotients are transformed into reciprocal integers, we can compute the digits of $\pi$ by using the Maclaurin series expansion (3) of the arctangent function. However, our empirical results show that the following two expansion series

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

(17)

and

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

(18)

where

$$g_1\left(x\right)=2/x,h_1\left(x\right)=1,$$

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

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

can be used more efficiently for computing digits of $\pi$ due to more rapid convergence.

Equation (17) is known as Euler’s series expansion. It is interesting to note that this series expansion can be derived from the following integral [19]

$$arctan\left(x\right)={\int }_0^{\pi /2}{\displaystyle \frac{xsinu}{1+x^2}}{\displaystyle \frac{1}{1-\frac{x^2{sin}^{2}u}{1+x^2}}}du.$$

Equation (18) can be derived by substituting

$$f(x,t)={\displaystyle \frac{x}{2}}\left({\displaystyle \frac{1}{1+ixt}}+{\displaystyle \frac{1}{1-ixt}}\right)$$

into the identity

$${\int }_0^{1}f(x,t)dt={\left.2\sum _{m=1}^M\sum _{n=0}^{\infty }{\displaystyle \frac{1}{{\left(2M\right)}^{2n+1}(2n+1)!}}{\displaystyle \frac{{\partial }^{2n}}{\partial t^{2n}}}f(x,t)\right|}_{t={\displaystyle \frac{m-1/2}{M}}}$$

that we proposed and used earlier (see [20] and literature therein for more details). Computational test we performed shows that the arctangent series expansion (18) is significantly faster in convergence than the arctangent series expansion (17).

At $k=6$ with help of the two-step iteration (8), we obtain

$${\theta }_{1,6}=-{\displaystyle \frac{2634699316100146880926635665506082395762836079845121}{38035138859000075702655846657186322249216830232319}}.$$

Consequently, the two-term Machin-like formula for $\pi$ becomes

$$\begin{array}{cc}\hfill {\displaystyle \frac{\pi }{4}}& =32arctan\left({\displaystyle \frac{1}{40}}\right)\hfill \\ & -arctan\left({\displaystyle \frac{38035138859000075702655846657186322249216830232319}{2634699316100146880926635665506082395762836079845121}}\right).\hfill \end{array}$$

As we can see, this equation contains a quotient with large number of digits in numerator and denominator. We may attempt to reduce the number of digits by approximation. Unfortunately, the two-step iteration (8) is not efficient in approximating ${\theta }_{1,k}$. Any our attempt to approximate ${\theta }_{1,k}$ by two-step iteration (8) either do not provide a desired accuracy or completely diverge from the value ${\theta }_{1,k}$. This makes approximation inefficient and unpredictable especially at lager values of the integer k.

Moreover, the number of digits rapidly grows with increasing k. For example, at $k=27$ and $M=0$ equation (12) results in

$$\begin{array}{cc}\hfill {\displaystyle \frac{\pi }{4}}& =2^{27-1}arctan\left({\displaystyle \frac{1}{{\gamma }_{27}}}\right)+arctan\left({\displaystyle \frac{1}{{\theta }_{1,27}}}\right)\hfill \\ & =67108864arctan\left({\displaystyle \frac{1}{85445659}}\right)-arctan\left({\displaystyle \frac{9732933578\dots 4975692799}{2368557598\dots 9903554561}}\right),\hfill \end{array}$$

where argument in the second arctangent function contains $522,185,807$ digits in the numerator and $522,185,816$ digits in the denominator. In the recent publication, Gasull et al. showed that at $k=31$ number of digits in the numerator and denominator in the two-term Machin-like formula for $\pi$ of kind (5) are $9,647,887,023$ and $9,647,887,033$, respectively (see Table 2 in [7]). Consequently, the two-step iteration (8) is appeared to be impractical for approximating the constants ${\theta }_{1,k}$.

This problem can be effectively resolved by using a new method of iteration that will be shown in the next section.

3. An Iterative Formula to Compute $\pi$

A problem that appears with two-step iteration (8) is a rapidly growing number of digits. In particular, the number of digits in ${\kappa }_n$ and ${\lambda }_n$ doubles at each consecutive increment of the index n. This may also restrict the application of the two-step iteration (8). The following theorem shows a new iteration technique as an alternative to the double-step iteration (8).

Theorem 2. 

There is a two-term Machin-like formula

$${\displaystyle \frac{\pi }{4}}=arctan\left({\displaystyle \frac{1}{u}}\right)+arctan\left({\displaystyle \frac{u-1}{u+1}}\right).$$

Proof. 

The proof is straightforward and immediately follows from the following identity

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

(19)

Specifically, assuming

$$x={\displaystyle \frac{1}{u}}$$

and

$$y={\displaystyle \frac{u-1}{u+1}},$$

according to identity (19) we have

$${\displaystyle \frac{1/u+(u-1)/(u+1)}{1-1/u\left((u-1)/(u+1)\right)}}=1.$$

This completes the proof since

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

   □

Consider now the next theorem.

Theorem 3. 

There is a relation

$$2^{k-1}arctan\left({\displaystyle \frac{1}{{\gamma }_k}}\right)=arctan\left({\displaystyle \frac{1}{v_k}}\right),$$

where

$$v_n={\displaystyle \frac{1}{2}}\left(v_{n-1}-{\displaystyle \frac{1}{v_{n-1}}}\right),n=\{2,3,\dots k\}$$

such that

$$v_1={\gamma }_k.$$

Proof. 

Since

$$2^{0}arctan\left({\displaystyle \frac{1}{{\gamma }_1}}\right)=arctan\left({\displaystyle \frac{1}{v_1}}\right),$$

it follows that

$$2^{1}arctan\left({\displaystyle \frac{1}{{\gamma }_2}}\right)=arctan\left({\displaystyle \frac{1}{v_2}}\right),$$

where

$${\displaystyle \frac{1}{v_2}}={\displaystyle \frac{2/v_1}{1-1/v_1^2}}.$$

Similarly, at $k=3$ we get

$$2^{2}arctan\left({\displaystyle \frac{1}{{\gamma }_3}}\right)=arctan\left({\displaystyle \frac{1}{v_3}}\right),$$

where

$${\displaystyle \frac{1}{v_3}}={\displaystyle \frac{2/v_2}{1-1/v_2^2}}.$$

Therefore, by induction for an arbitrary integer k we can write

$$2^{k-1}arctan\left({\displaystyle \frac{1}{{\gamma }_k}}\right)=arctan\left({\displaystyle \frac{1}{v_k}}\right),$$

where

$${\displaystyle \frac{1}{v_k}}={\displaystyle \frac{2/v_{k-1}}{1-1/v_{k-1}^2}}.$$

The equation above can be rewritten as

$$v_k={\displaystyle \frac{1}{2}}\left(v_{k-1}-{\displaystyle \frac{1}{v_{k-1}}}\right).$$

(20)

and this proves the theorem.    □

Thus, from this and above theorems it follows that

$${\displaystyle \frac{\pi }{4}}=2^{k-1}arctan\left({\displaystyle \frac{1}{v_1}}\right)+arctan\left({\displaystyle \frac{v_k-1}{v_k+1}}\right).$$

(21)

where

$${\displaystyle \frac{v_k-1}{v_k+1}}={\displaystyle \frac{1}{{\theta }_1}}.$$

The iterative formula (20) and equation (21) are our main results. Despite simplicity, their significance can be demonstrated in computing digits of $\pi$ and nested radicals consisting of square roots of 2.

4. Approximation Methodologies

4.1. A Rational Approximation of $\pi$

Recently, we have reported a rational approximation for $\pi$ by approximating the two-term Machin-like formula (13) [21]. Here we show how a rational approximation for $\pi$ can be constructed by using the equation (21) based on iterative formula (20).

Consider the following theorem.

Theorem 4. 

There is a limit such that

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

Proof. 

This is Ramanujan’s nested radical. Since

$$c_k=\underset{k\mathrm{square}\mathrm{roots}}{\underbrace{\sqrt{2+\sqrt{2+\sqrt{2+\dots +\sqrt{2}}}}}}=\sqrt{2+c_{k-1}},$$

we can write

$$\underset{k\to \infty }{lim}{c}_k=\underset{k\to \infty }{lim}\sqrt{2+c_{k-1}}=\underset{k\to \infty }{lim}\sqrt{2+c_k}.$$

Let

$$z=\underset{k\to \infty }{lim}{c}_k.$$

Then, solving the equation

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

or

$$z^2-z-2=0,$$

we get two solutions $z_1=-1$ and $z_2=2$. Since the sequence

$$\{c_0,c_1,c_2,c_3,\dots \}=\left\{0,\sqrt{2},\sqrt{2+\sqrt{2}},\sqrt{2+\sqrt{2+\sqrt{2}}},\dots \right\}$$

monotonically increases and positive at any $c_n$ except $c_0=0$, we have to exclude the negative solution. Consequently, the theorem is proved.    □

From the theorem 4 we can make a conclusion that

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

(22)

since due to relation

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

leading to

$$\underset{k\to \infty }{lim}\sqrt{2-c_k}=\sqrt{2-{lim}_{k\to \infty }{c}_k}=0$$

the limit (22) follows.

Next, will we try to find a rational approximation for computing digits of $\pi$. This can be done by using equation (21). In order to approximate it, we need to consider two lemmas below.

Lemma 1. 

There is a limit such that

$${\displaystyle \frac{\pi }{4}}=\underset{k\to \infty }{lim}{2}^{k-1}arctan\left({\displaystyle \frac{1}{{\gamma }_k}}\right).$$

Proof. 

By definition

$${\displaystyle \frac{c_k}{\sqrt{2-c_{k-1}}}}-1<{\gamma }_k<{\displaystyle \frac{c_k}{\sqrt{2-c_{k-1}}}}.$$

Since

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

while

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

we can infer that

$$\underset{k\to \infty }{lim}{\displaystyle \frac{c_k}{\sqrt{2-c_{k-1}}}}=\infty .$$

Therefore, we can write

$$\underset{k\to \infty }{lim}{\displaystyle \frac{{\gamma }_k}{\sqrt{2-c_{k-1}}/c_k}}=1$$

from which it follows that

$$\underset{k\to \infty }{lim}arctan\left({\displaystyle \frac{\sqrt{2-c_{k-1}}}{c_k}}\right)=\underset{k\to \infty }{lim}arctan\left({\displaystyle \frac{1}{{\gamma }_k}}\right)$$

and proof of this lemma follows since [22]

$$\underset{k\to \infty }{lim}{2}^{k-1}arctan\left({\displaystyle \frac{\sqrt{2-c_{k-1}}}{c_k}}\right)={\displaystyle \frac{\pi }{4}}.$$

(23)

   □

Lemma 2. 

There is a limit

$${\displaystyle \frac{\pi }{4}}=\underset{k\to \infty }{lim}\left({\displaystyle \frac{2^{k-1}}{{\gamma }_k}}+{\displaystyle \frac{1}{{\theta }_{1,k}}}\right)=\underset{k\to \infty }{lim}\left({\displaystyle \frac{2^{k-1}}{v_1}}+{\displaystyle \frac{v_k-1}{v_k+1}}\right).$$

Proof. 

The proof follows immediately from the lemma 1 and due to equation (21).    □

The importance of the lemma 2 can be seen by comparing the accuracy of two approximations

$${\displaystyle \frac{\pi }{4}}\approx {\displaystyle \frac{2^{k-1}}{{\gamma }_k}}={\displaystyle \frac{2^{k-1}}{v_1}}$$

(24)

and

$${\displaystyle \frac{\pi }{4}}\approx {\displaystyle \frac{2^{k-1}}{{\gamma }_k}}+{\displaystyle \frac{1}{{\theta }_{1,k}}}={\displaystyle \frac{2^{k-1}}{v_1}}+{\displaystyle \frac{v_k-1}{v_k+1}}.$$

(25)

Since $v_k$ is very close to unity, the double-term rational approximation (25) can be simplified as

$${\displaystyle \frac{\pi }{4}}\approx {\displaystyle \frac{2^{k-1}}{v_1}}+{\displaystyle \frac{v_k-1}{2}}.$$

(26)

Mathematica code below shows how number of digits of $\pi$ can be computed by using a single- and double-term rational approximations (24) and (26), respectively.

Clear[k, c, flNum, accNum, v1, vk];

(* Integer k *)

k = 3000;

(* Use if needed: $RecursionLimit = 100000; *)

$RecursionLimit = 50000;

c[0] := c[0] = 0;

c[n_] := c[n] = SetAccuracy[Sqrt[2 + c[n - 1]], k];

(* Floor number *)

flNum = Floor[c[k]/Sqrt[2 - c[k - 1]]];

(* Accuracy number *)

accNum = Length[RealDigits[flNum][[1]]];

(* Coefficient v_1 *)

v1 = SetAccuracy[flNum, 2*accNum];

Print["At k = ", k, " number of digits of \[Pi] with single term: ",

  MantissaExponent[\[Pi] - 4 (2^(k - 1)/v1)][[2]] // Abs];

(* Compute coefficient v_k *)

vk = v1;

Do[vk = 1/2 (vk - 1/vk), k - 1];

Print["At k = ", k, " number of digits of \[Pi] with two terms: ",

    MantissaExponent[\[Pi] - 4 (2^(k - 1)/v1

    + (vk - 1)/2)][[2]] //Abs];

This code generates the following output:

At k = 3000 number of digits of À with single term: 902

At k = 3000 number of digits of À with two terms: 1805

As we can see from this example, once we know the value of the constant $v_k$, at $k=3000$ the number of digits of $\pi$ is doubled from 902 to 1805.

4.2. An Approximation of $\pi$ with Cubic Convergence

Now we show how to obtain a formula for $\pi$ with cubic convergence. Consider the following theorem.

Theorem 5. 

The two-term Machin-like formula (13) for π can be represented in trigonometric form as

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

Proof. 

Applying de Moivre’s theorem we can write

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

Substituting equation above into equation (10) and taking into consideration that ${\theta }_{1,k}=\beta$ we get

$${\theta }_{1,k}={\displaystyle \frac{cos\left(2^{k-1}\mathrm{Arg}({\kappa }_1+i{\lambda }_1)\right)}{1-sin\left(2^{k-1}\mathrm{Arg}({\kappa }_1+i{\lambda }_1)\right)}}.$$

Since

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

it follows that

$$\mathrm{Arg}({\kappa }_1+i{\lambda }_1)=arctan\left({\displaystyle \frac{2{\gamma }_k}{{\gamma }_k^2-1}}\right),x>0.$$

Consequently, the equation (10) can be represented as

$${\theta }_{1,k}={\displaystyle \frac{cos\left(2^{k-1}arctan\left(\frac{2{\gamma }_k}{{\gamma }_k^2-1}\right)\right)}{1-sin\left(2^{k-1}arctan\left(\frac{2{\gamma }_k}{{\gamma }_k^2-1}\right)\right)}}$$

(27)

and this completes the proof.    □

We may attempt to approximate the coefficient ${\theta }_{1,k}$. Since

$${\displaystyle \frac{2{\gamma }_k}{{\gamma }_k^2-1}}\to {\displaystyle \frac{2}{{\gamma }_k}}$$

with increasing k, the equation (27) can be approximated as

$${\theta }_{1,k}\approx {\displaystyle \frac{cos\left(2^{k-1}arctan\left(\frac{2}{{\gamma }_k}\right)\right)}{1-sin\left(2^{k-1}arctan\left(\frac{2}{{\gamma }_k}\right)\right)}}.$$

Since

$$arctan\left(x\right)\approx x,\left|x\right|\ll 1,$$

the approximation above can be further simplified as

$${\theta }_{1,k}\approx {\displaystyle \frac{cos\left(\frac{2^k}{{\gamma }_k}\right)}{1-sin\left(\frac{2^k}{{\gamma }_k}\right)}}.$$

(28)

Using the following trigonometric identities

$$cos\left(x\right)={\displaystyle \frac{2tan(x/2)}{1+{tan}^2(x/2)}}$$

and

$$sin\left(x\right)={\displaystyle \frac{1-{tan}^2(x/2)}{1+{tan}^2(x/2)}}$$

the equation (28) can be expressed as

$${\theta }_{1,k}\approx {\displaystyle \frac{1+{tan}^2\left(2^{k-1}/{\gamma }_k\right)}{1-tan\left(2^{k-1}/{\gamma }_k\right)}}.$$

This leads to

$${\displaystyle \frac{1}{{\theta }_{1,k}}}\approx {\displaystyle \frac{1-tan\left(2^{k-1}/{\gamma }_k\right)}{1+{tan}^2\left(2^{k-1}/{\gamma }_k\right)}}.$$

(29)

Since

$${tan}^2\left(2^{k-1}/{\gamma }_k\right)$$

converges faster to 1 than

$$tan\left(2^{k-1}/{\gamma }_k\right),$$

we can simplify approximation (29) as

$${\displaystyle \frac{1}{{\theta }_{1,k}}}\approx {\displaystyle \frac{1}{2}}\left(1-tan\left(2^{k-1}/{\gamma }_k\right)\right).$$

(30)

Substituting approximation (30) into equation (13), we can approximate the two-term Machin-like formula (1) for $\pi$ as follows

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

Consequently, from the lemma 2 we have

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

(31)

Consider now the following limit

$$\underset{x\to \pi /4}{lim}{\displaystyle \frac{2{sin}^2\left(x\right)}{tan\left(x\right)}}=1.$$

Since numerator of this limit is very close to its denominator near vicinity of $\pi /4$, we may replace tangent function in approximation (31) with twice of square of sine function. This leads to

$${\displaystyle \frac{\pi }{4}}\approx {\displaystyle \frac{2^{k-1}}{{\gamma }_k}}+{\displaystyle \frac{1}{2}}\left(1-2{sin}^2\left({\displaystyle \frac{2^{k-1}}{{\gamma }_k}}\right)\right)$$

or

$${\displaystyle \frac{\pi }{4}}\approx {\displaystyle \frac{2^k}{2{\gamma }_k}}+{\displaystyle \frac{1}{2}}cos\left({\displaystyle \frac{2^k}{{\gamma }_k}}\right)$$

or

$${\displaystyle \frac{\pi }{2}}\approx {\displaystyle \frac{2^k}{{\gamma }_k}}+cos\left({\displaystyle \frac{2^k}{{\gamma }_k}}\right).$$

(32)

Equation (32) provides a hint for computation of $\pi$. Using an heuristic assumption by change of the variable such that

$${\displaystyle \frac{2^k}{{\gamma }_k}}\to a_1,k\ge 1,$$

we may attempt to compute $\pi$ by using the following iteration

$$a_n=a_{n-1}+cos\left(a_{n-1}\right)$$

(33)

resulting in

$$\underset{n\to \infty }{lim}{a}_n={\displaystyle \frac{\pi }{2}}.$$

The following is a Mathematica code for computing digits of $\pi$ with help of the iterative formula (33).

Clear[a, accNum, n];

(* Initial accuracy number *)

accNum = 10;

(* Initial value of À/2 *)

a = SetAccuracy[3.145926/2, accNum];

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

Print["Iteration n", "  |  ", "Number of digits of À"];

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

(* Iteration *)

For[n = 1, n <= 8, a = SetAccuracy[a + Cos[a], accNum];

    Print[n, "            |  ", MantissaExponent[Pi

        - 2*a][[2]]//Abs]; accNum = 5*accNum; n++];

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

The output of this code:

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

Iteration n | Number of digits of À

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

1           | 8

2           | 26

3           | 76

4           | 232

5           | 698

6           | 2095

7           | 6288

8           | 18868

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

shows that the iterative formula (33) provides cubic convergence since number of digits of $\pi$ increases by a factor of 3 after each iteration.

4.3. A Numerical Solution for the Nested Radicals with Roots of 2

As we can see, approximation (26) doubles the number of digits of $\pi$. However, this approach can be used not only for computing digits of $\pi$. It can also be used to compute nested radical consisting of square roots of 2. In particular, by using iterative formula (20) we can compute the nested radical with roots of 2 as follows

$${\displaystyle \frac{v_k}{v_{k-n}}}\approx {\displaystyle \frac{\sqrt{2-c_n}}{c_{n+1}}}.$$

Mathematica code below shows how to generate digits of the nested radicals with roots of 2 at $n=3$

$${\displaystyle \frac{v_k}{v_{k-3}}}\approx {\displaystyle \frac{\sqrt{2-c_3}}{c_4}}={\displaystyle \frac{\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2}}}}}{\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2}}}}}}.$$

Clear[k, c, flNum, accNum, v]

(*Assign value of k*)

k = 5000;

(* Increase if needed: $RecursionLimit = 100000; *)

$RecursionLimit = 50000;

(* Define nested radicals *)

c[0] := c[0] = 0;

c[n_] := c[n] = SetAccuracy[Sqrt[2 + c[n - 1]], k];

(* Compute floor number *)

flNum = Floor[c[k]/Sqrt[2 - c[k - 1]]];

(*Set accuracy with accuracy number*)

accNum = Length[RealDigits[flNum][[1]]];

(*Compute v_k*)

v[1] := v[1] = SetAccuracy[Floor[flNum], accNum];

v[n_] := v[n] = 1/2 (v[n - 1] - 1/v[n - 1]);

Print[MantissaExponent[v[k]/v[k - 3] - Sqrt[2 - c[3]]/c[4]][[2]] //

   Abs, " computed digits of nested radical"];

This code produces the following output:

1506 computed digits of nested radical

It is interesting to note that due to relation

$${\displaystyle \frac{v_k}{v_{k-1}}}\approx {\displaystyle \frac{\sqrt{2-\sqrt{2}}}{\sqrt{2+\sqrt{2}}}}=\sqrt{2}-1$$

we can compute $\sqrt{2}$ and its total numbers of correct digits by using the following command lines:

Print["Computed square root of 2 is ", N[1 + v[k]/v[k - 1], 20],"..."];

Print[MantissaExponent[(v[k]/v[k - 1] + 1) - Sqrt[2]][[2]] //

    Abs, " computed digits of square root of 2"];

,

This Mathematica command line generates the following output:

Computed square root of 2 is 1.4142135623730950488...

1506 computed digits of square root of 2

There are several different methods [23,24,25] that can be used for computing digits of $\sqrt{2}$. However, the method of computation that we developed here can be implemented beyond $\sqrt{2}$ for efficient computation of the nested radicals with roots of 2 of kind

$${\displaystyle \frac{\sqrt{2-c_{n-1}}}{c_n}}.$$

These nested radicals are utilized in formulas like (23) and

$${\displaystyle \frac{\pi }{4}}=\sum _{k\ge 2}{2}^{k-1}arctan\left({\displaystyle \frac{\sqrt{2-c_{k-1}}}{c_k}}\right).$$

Furthermore, due to relation

$${\displaystyle \frac{\sqrt{2-c_{n-1}}}{c_n}}\approx {\displaystyle \frac{1}{2}}\sqrt{2-c_{n-1}},n\gg 1,$$

this technique can also be applied for computation of the nested radicals with roots of 2 of kind

$$\sqrt{2\pm \sqrt{2+\sqrt{2+\sqrt{2+\dots }}}}.$$

as an alternative to other known methods [10,16,26,27].

4.4. An Approximation to $\pi$ with Arbitrary Convergence

Consider how equation (21) can be used for computing digits of $\pi$. Rather than to apply it directly for computing digits of $\pi$, we can apply equation (12) to transform quotient

$${\displaystyle \frac{1}{{\theta }_{1,k}}}={\displaystyle \frac{v_k-1}{v_k+1}}$$

(34)

into reciprocal integers (compare equations (15) and (16) above before and after transformation). Another approach to compute $\pi$ is to approximate ${\theta }_k$ directly during process of iteration.

The following Mathematica code shows how to compute digits of $\pi$ by approximating equation (34) in iteration process.

Clear[k, accNum, c, v1, vk, n]

(* Integer k *)

k = 50;

(* Define array of accuracy numbers *)

accNum = {100000, 200000, 300000, 400000, 500000};

(* Define nested radicals *)

c[0] := c[0] = 0;

c[n_] := c[n] = SetAccuracy[Sqrt[2 + c[n - 1]], k];

n = 1;

Do[

  (* Setting accuracy *)

  v1 = SetAccuracy[Floor[c[k]/Sqrt[2 - c[k - 1]]], accNum[[n]]];

  (* Computing v_k *)

  vk = v1; Do[vk = 1/2 (vk - 1/vk), k - 1];

  Print["n = ", n, ", ", MantissaExponent[\[Pi] -

      4*(2^(k - 1)*ArcTan[1/v1] + ArcTan[(vk - 1)/(vk + 1)])][[2]] //

          Abs, " digits of À"];

  n++, 5];

The output of this code is

n = 1, 100014 digits of À

n = 2, 200014 digits of À

n = 3, 300014 digits of À

n = 4, 400014 digits of À

n = 5, 500014 digits of À

The number of digits in the Mathematica code above is determined by list variable accNum consisting of a sequence of 5 numbers as given by

accNum = {100000, 200000, 300000, 400000, 500000};

If we increase these numbers, then the number of correct digits of $\pi$ increases by the same factors. From this example we can see that by increasing parameters in the list variable accNum, we can archive an arbitrary convergence rate.

5. Conclusions

An iterative formula (20) that can be used for computing $\pi$ and nested radicals with roots of 2 is derived. It is shown how this iterative formula can be implemented to generate and approximate the two-term Machin-like formulas for $\pi$. Developing this approach, a simple trigonometric formula for $\pi$ with cubic convergence can be obtained. Some examples with Mathematica codes are provided.