Algorithmic determination of a large integer in the two-term Machin-like formula for pi

In our earlier publication we have shown how to compute by iteration a rational number ${u_{2,k}}$ in the two-term Machin-like formula for pi of kind $$\frac{\pi}{4}=2^{k-1}\arctan\left(\frac{1}{u_{1,k}}\right)+\arctan\left(\frac{1}{u_{2,k}}\right),\qquad k\in \mathbb{Z},\quad k\ge 1,$$ where ${u_{1,k}}$ can be chosen as an integer ${u_{1,k}} = \left\lfloor{{a_k}/\sqrt{2-a_{k-1}}}\right\rfloor$ with nested radicals defined as ${a_k}=\sqrt{2+a_{k-1}}$ and $a_0 = 0$. In this work we report an alternative method for determination of the integer $u_{1,k}$. This approach is based on a simple iteration and does not require any irrational (surd) numbers from the set $\left\{a_k\right\}$ in computation of the integer $u_{1,k}$. Mathematica programs validating these results are presented.


Introduction
Historically, a computation of decimal digits of π was a big challenge until 1706, when the English astronomer and mathematician John Machin discovered a two-term formula for π as given by π 4 = 4 arctan that is named in his honor now. Using this remarkable formula he first was able to calculate 100 decimal digits of π [1][2][3]. Nowadays the identities of kind where A j and B j are rational numbers, are regarded as the Machin-like formulas for π. Interestingly that some of them, including the original Equation (1), can be proved geometrically [4,5]. It is very often when in the Machin-like formulas for π the constants A j and B j are both integers [6][7][8]. The more complete lists of formulas of kind (2) can be found in the references [9,10] and weblinks provided therein. The significance of the Machin-like formulas cannot be overestimated as their application may be one of the most efficient ways in computing π. Historically, only these formulas were able to compete with Chudnovsky formula [2,11] to beat the records in computation of the decimal digits of π. In particular, in 2002, Kanada first computed more than one trillion digits of π by using the following self-checking pair of the Machin-like formulas [12]  Such an achievement made by Kanada shows a colossal potential of the Machin-like formulas for computation of the decimal digits of π.
As a simplest case one can apply the Maclaurin expansion series for computation of the arctangent functions in Equation (2) arctan (x) = x − x 3 3 + x 5 5 − x 7 7 + · · · = ∞ n=0 (−1) n x 2n+1 2n + 1 and since this equation implies we can conclude that it would be very desirable to have the coefficient B j as large as possible by absolute value in order to improve the convergence rate.
Although the Maclaurin expansion series of the arctangent function can be simply implemented, its application is not optimal. The more efficient way to compute the arctangent functions in Equation (2) is to use the Euler's expansion formula [13] arctan Alternatively, the following expansion series can also be used for more rapid convergence. This formula can be obtained by a trivial rearrangement of the Equation (5) from our work [14] (see also [15]). In 1938 Lehmer introduced a measure defined as [6,16] This measure can be used to determine a computational efficiency of a given Machin-like formula for π. Specifically, when value of the constant µ is smaller, then less computational labour is required to compute π by a given Machin-like formula. Therefore, it is very desirable to reduce the number of the terms J and to increase B j by absolute value. The more detailed information about the Lehmer's measure µ and its significance for efficient computation of π can be found in literature [7].
In our previous publication using de Moivre's formula we derived the two-term Machin-like formula for π [15] where u 1,k can be chosen as an integer such that a set of nested radicals {a k } can be computed as a k = √ 2 + a k−1 starting from a 0 = 0. In the recent publication the researcher(s) from the Wolfram Mathematica [17] demonstrated an example for computation of π as a rational fraction. In particular, it was shown that with integer constant u 1,1000 the constant π can be computed such that the ratio results in more than 300 correct decimal digits. However, the method of computation shown in [17] is based on Equation (5) that involves the nested radicals consisting of multiple square roots of 2 [18][19][20][21][22][23]. Although Equation (5) helps generate the required integers u 1,k , it should not be generally used at larger values k since a function based on multiple square roots is not an elementary. Therefore, a simple method based on rational approximation would be preferable. In this work, we develop a new method of computation of the integer u 1,k that excludes application of the set of nested radicals {a k }. This approach is simple and does not require any irrational (surd) numbers in computation.

Preliminaries
Suppose that The simplest case when α = γ = 1. However, there may be infinitely many identities of kind (6). Let us show the infinitude of this kind of formulas.
Theorem 2.1. There are infinitely many numbers α and γ satisfying the relation (6).
The derivation of Equation (7) is very simple and can be found in [14].
It is interesting to note that using Equation (7) one can easily prove the well-known Equation (8) for π below.
From this limit it immediately follows that Consequently, the ratio The argument of the arctangent function in Equation (7) tends to zero as the integer k increases. Therefore, in accordance with relation (3) we can write from which it follows that As we can see, this proof of Equation (8) is as easy as the one shown in [18]. We also note that the following limit that we used in the proof is valid since the identity (7) remains valid at any arbitrarily large positive integer k. It is also easy to prove the infinitude of the Machin-like formulas (2) for π; substituting x = 1 into the expansion series [24] arctan we obtain the following identity [25] leading to π 4 = arctan (1) , N = 1, and so on. Although the identity (9) shows infinitude of the Machin-like formulas for π, its number of the terms increases with increasing N. We can also show a simple proof for infinitude of the two-term Machin-like formulas for π. Lemma 2.3. For real α and β 1 , there are infinitely many two-terms Machinlike formulas for π of kind Proof. The Lemma 2.3 follows directly from the Theorem 2.1 that implies infinitude of equations of kind (6). In order to show this relation, we assume that α and γ in Equation (6) are both positive numbers and represent γ as a sum γ = ζ + δ, where δ is any small number that can be chosen arbitrarily such that γ >> |δ|. Thus, we can rewrite the Equation (6) in form From the inequality γ >> |δ| it follows that ζ ≈ β. Therefore, we can approximate By introducing now an error term ε, we can infer that Consequently, we have where η is defined such that The Equation (11) is of the same kind as that of given by Equation (10). This completes the proof since for any equation of kind (6) we can always construct an equation of kind (10).
When α and β 1 in Equation (10) are known, then the unknown value β 2 is given by The derivation of Equation (12) can be shown from the following identity Thus, substituting this identity into Equation (10) results in Exponentiation on both sides leads to Solving this with respect to the constant β 2 leads to Equation (12). (10) ∀k ≥ 2 the multiplier α = 2 k−1 and β 1 is a rational number greater than 1, then β 2 is also a rational number.

Theorem 2.4. If in Equation
Proof. Define σ 1 and τ 1 such that Then, it is not difficult to see by induction that where by the following two-step iteration we have Consequently, Equation (12) can be rewritten in form and, therefore, the principal value argument can be replaced by the arctangent function as follows Arg Consequently, we can write As we can see from this equation, the constant β 2 ∈ R. This signifies that the imaginary part of Equation (15) must be equal to zero. Therefore, from Equation (15) we get The values σ k and τ k are rational since, according to two-step iteration (14) all values σ n and τ n at any intermediate steps of iterations are rational. Therefore, the constant β 2 must be a rational number.
Chien-Lih proposed a method showing how to reduce the Lehmer's measure by using the Euler's-type identity in an iteration for generating the two-term Machin-like formulas for π. However, our method of generating the two-term Machin-like formula for π based on the two-step iteration (14) is much easier than the method proposed by Chien-Lih in the work [26].

Derivation
Using Equation (16) we can rewrite the Equation (10) as In general, the constant β 1 may be either rational or irrational number. However, it is more convenient to apply notation u 1,k that is defined by Equation (5) instead of β 1 . Such a notation is to emphasize that the constant u 1,k is an integer dependent upon on k. Thus, with this notation the two-term Machin-like formula for π can be represented as where in accordance with Equation (12) we have now It is interesting to note that by taking k = 3, we get Substituting u 1,3 = 5 into Equation (17) we obtain u 2,3 = −239. Considering that 2 k−1 = 2 2 = 4 and substituting these two constants into Equation (10) we derive an original Machin-like formula (1) for π.
Since the value 2 k−1 rapidly increases with increasing k, application of the Equation (18) if not effective to compute the second constant u 2,k . However, the two-step iteration (14) perfectly resolves this issue. Specifically, implying that the initial values for the two-step iteration (14) are we can find the second constant as We can derive again the original Machin-like formula (1) for π by using the two-step iteration (14) at k = 3. This leads to the following Finally, using Equation (19) we can find the second constant to be It should be noted that the second constant u 2,k is an integer only at k = 2 and k = 3. At k > 3 it is not an integer but a rational number.
The next example is k = 6. The first constant is an integer given by The second constant u 2,6 is a rational number that can be computed either by Equation (18) or, more efficiently, by two-step iteration (14) u 2,6 = − 2634699316100146880926635665506082395762836079845121 38035138859000075702655846657186322249216830232319 .
The following Mathematica code validates the two-term Machin-like formula for π at k = 6: Alternatively, the second constant can also be found by using the following identity in trigonometric form that follows from Equation (16). It should be noted that the constant u 2,k must be a rational number as it has been shown by Theorem 2.4. We can see that from Equation (5) it follows that u 1,k << u 2 1,k , k >> 1.
Consequently, the following ratio can be simplified as given by Replacing the arguments of the sine and cosine functions in Equation (20), we can approximate the two-term Machin-like formula (17) for π as Using the identities for the double angle sin (x) = 2 tan (x/2) 1 + tan 2 (x/2) , after some trivial rearrangements we obtain Recently, it has been noticed in publication [17] that the ratio 2 k−1 /u 1,k approximates π/4 reasonably well when integer k = 1000. The following theorem shows why accuracy of this ratio improves with increasing k. Theorem 3.1. There is a limit Proof. By definition of the floor function we have where by definition the fractional part cannot be smaller than zero and greater than or equal to unity 0 ≥ frac a k √ 2 − a k−1 < 1.
Therefore, we can write Since the fractional part cannot be smaller than 0 and greater than 1 we can conclude that Consequently, we can infer that From Theorem 2.2 we know that Consequently, from the relation (3) we get and the limit (22) follows.
Lemma 3.2. There is a limit such that Proof. We know that the limit is valid since the identity (17) remains valid at any arbitrarily large integer k. We also know that from the Theorem 3.1 and relation (3) it follows that This signifies that lim k→∞ 2 k−1 arctan 1 However, according to relation (3) this equation can be simplified as and the proof follows.
Applying Theorem 3.1 to the left side of approximation (21) yields .
We note that both arguments of the arctangent function tend to zero with increasing k. Therefore, referring to the relation (3) again we can simplify the approximation above as . Using the Theorem 3.1 it follows that at k → ∞ Consequently, the value with increasing k. This leads to Comparing Equations (17) and (21) one can see that and due to relation (23) this approximation can be further simplified to Although this equation only approximates the second constant u 2,k , its accuracy, nevertheless, improves with increasing k. Perhaps, the approximation (25) can also be used at larger values of the integer k as an alternative to the exact formula (19) based on the two-step iteration (14).
We can see consistency of the approximations (23) and (25) with Theorem 3.1 and Lemma 3.2. In particular, when k tends to infinity the constants u 1,k and u 2,k also tend to infinity. Therefore, in order to enhance a convergence rate, it is important to obtain the integer u 1,k in the two-term Machin-like formula (17) for π as large as possible. Once the value of the first constant u 1,k is determined, the second constant u 2,k can be computed by using Equation (19) based on two-step iteration formula (14). For example, at k = 27 the value u 1,27 = 85,445,659. The corresponding Lehmer's measure is µ ≈ 0.245319 only. Such a small Lehmer's measure implies a rapid convergence rate. In particular, we can observe 16 correct decimal digits of π per term increment. This can be confirmed by running a Mathematica program provided in [27].
At k = 27 Equation (19) yields a rational number u 2,27 consisting of 522,185,807 digits in numerator and 522,185,816 digits in denominator. Such a quotient with huge numbers in numerator and denominator is not unusual and can also be observed in Borwein integrals. Specifically, Bäsel and Baillie in their work [28] showed that a formula for π can be generated with a quotient consisting of 453,130,145 and 453,237,170 digits in its numerator and denominator, respectively. The interested readers can download the exact number u 2,27 with all digits from [29].

Implementation
At first glance, the approximation (24) does not look interesting as its both sides contain the constant u 1,k and it is unclear how to represent it in explicit form. However, sample computations we performed with this approximation show that it can be implemented effectively. In particular, we noticed that application of approximation (24) in iteration provides a result that tends to be more accurate with increasing the integer k.
Consider for example k = 10. In this case we have that With initial guess for u 1,10 , say 1000, after just 5 iterations (self-substitutions) we obtain u  We have ceased the iterative process after 5-th cycle since two successive numbers coincide with each other at fourth and fifth iterations with same output 651.899. Application of the Equation (24) is convenient since for each consecutive increment of the integer k the following inequality remains valid This inequality follows from the property of the floor function that is used in Equation (5); the constant u 1,k+1 should be either equal to 2u 1,k or larger it by unity (see [17] for some examples). Thus, based on Equation (24) and inequality (26) we can make the following assumption .
The computational tests we performed shows that this formula provides correct results for a large range of the integer k ≥ 2. However, its general applicability for any arbitrarily large k yet to be proved.
There are different methods to approximate the tangent function in Equation (27). One of the ways is to truncate the following expansion series where B n are the Bernoulli numbers, defined by a contour integral Although this series expansion is rapid in convergence, its application may not be optimal since it requires the determination of the Bernoulli numbers. One of the ways to compute them is given by the following identity We can see that this formula involves the double summation and, therefore, cannot be rapid in principle especially at larger orders of n. Although other methods of computation of the Bernoulli numbers are more efficient, their implementations require quite sophisticated algorithms [30][31][32].
Alternatively, the tangent function may also be computed by using continued fractions [33][34][35]. However, algorithmic implementation of the continued fractions may not be optimal for our particular task.
This problem can be resolved by noticing that at each consecutive step of iteration the integer u 1,k increases, and because of this the argument of the tangent function decreases. Consequently, it may be reasonable to utilize argument reduction method for computation of the tangent function [36]. As a simplest case we can use, for example, the double angle identity providing argument reduction by a factor of two tan (2x) = 2 tan (x) 1 − tan 2 (x) .
This approximation implies that the arctangent function can be calculated in a simple iteration over and over again by defining the following function where Tangent function can be computed more accurately by defining since according to expansion series (28) we can also infer that tan (x) = x + x 3 /3 + O (x 5 ).
The following Mathematica command lines execute the program for computation of the integer constant u 1,k by using Equations (27) and (29): posed iteration formula (30) provides quadratic convergence to π without any irrational (surd) numbers.

Conclusions
In this work we propose a method for determination of the integer u 1,k . In particular, the algorithmic implementation of the Formula (27) shows that it can be used as an alternative to Equation (5) requiring a set of the nested radicals {a k } defined as a k = √ 2 + a k−1 and a 1 = 0. This method is based on a simple iteration and can be implemented without any irrational (surd) numbers.