An iterative method for computing π by argument reduction of the tangent function

In this work, we develop a new iterative method for computing the digits of π by argument reduction of the tangent function. This method combines a modiﬁed version of the iterative formula for π with squared convergence that we proposed in a previous work and a leading arctangent term from the Machin-like formula. The computational test we performed shows that algorithmic implementation can provide more than 17 digits of π per increment. Mathematica codes, showing the convergence rate for computing the digits of π , are presented


Introduction
In 1706, English astronomer and mathematician John Machin discovered a formula for π, expressed in terms of two arctangents: By using this formula, he was the first to compute 100 digits of π [1][2][3][4].Nowadays, identities of the kind where the coefficients α j and β j are either integers or rational numbers, are named after him as the Machin-like formulas for π.
The Maclaurin series expansion of the arctangent function is given by arctan from which it follows that arctan(x) = x + O(x 3 ).( 4) Consequently, we can conclude that the convergence of the Machin-like Formula (2) for π is better when the coefficients β j are larger by their absolute values.
In order to estimate the efficiency of a given Machin-like formula, Lehmer introduced the measure [5][6][7][8][9]: A smaller value of µ indicates a higher efficiency of a given Machin-like formula.According to this formula, the measure µ is smaller when the total number of the terms J is lower and the coefficients β j are higher by their absolute values.A more detailed description and the significance of the measure (5) can be found in the literature [6,9].Although Lehmer's paper [5] assumes that β j ∈ Q, it seems that the measure (5) is applicable only when all values of β j are integers.Otherwise, even if a single value from the set {β j } represents a quotient, its presence causes complexities in the computation [10], which appear due to the rapidly increasing number of digits in the numerators in summation terms when a series expansion of the arctangent function, like Equations ( 3), (30) or (31), is used in the computation.That is why it is very desirable to generate the Machin-like formulas for π where β j ∈ Z.
All coefficients in any Machin-like formula for π satisfy the relation implying that the real and imaginary parts of this product must be equal to each other.This relation can be used for validation.For example, there is a simple and elegant proof of the original Machin Formula (1) for π [11] (5 + i) 4 (259 + i) −1 = 2(1 + i).
Although several iterative formulas with quadratic, cubic, quartic, quintic and nonic convergences have been discovered [12][13][14][15][16], they require undesirable surd numbers that appear over and over again at each consecutive step of the iteration.
Historically, algorithms based on the Chudnovsky and Machin-like formulas successively broke records in computing the digits of π [4].Currently, the Chudnovsky formula, providing linear convergence with 14 to 16 digits of π per increment of the summation terms [1,2,4], appears to be the most efficient.
In 2002, Kanada broke a record by computing more than one trillion decimal digits of π for the first time by using the following self-checking pair (also known as the Störmer-Takano pair) of Machin-like formulas [4,17]  with Lehmer's measures of 1.58604 and 1.7799, respectively.Although the current record, achieved by using the Chudnovsky formula, exceeds one hundred trillion digits of π [10], application of the Machin-like formulas may be promising to calculate a comparable number of digits due to availability of more advanced and powerful supercomputers then those used by Kanada more than 20 years ago.Furthermore, more Machin-like formulas with smaller Lehmer's measures have been discovered [18][19][20] may be more efficiently used as a self-checking pair since their Lehmer's measures of 1.34085 and 1.39524, respectively, are considerably smaller.Therefore, application of Machin-like formulas with small Lehmer's measure have colossal potential and can be competitive for computing the digits of the constant π.Equation ( 7) was reported by Wetherfield [20].Equation ( .
by using the identity [21] arctan In our earlier publication [22], we derived the following two-term Machinlike formula where the constant γ may be conveniently chosen according to a relation [23] The complete proof of the identity ( 9) is quite lengthy and, therefore, beyond the scope of the present work.However, the detailed derivation of identity (9), which is available in [22], can be briefly outlined as a determination of the value η at a given k and γ in the two-term Machin-like formula of kind and a subsequent reformulation of Equation (11) into Equation ( 9) with the help of de Moivre's formula (cos(x) + i sin(x)) n = cos(nx) + i sin(nx).
It is not difficult to prove that if the first constant γ of Equation ( 11) is an integer or a rational number, then the second constant η must also be a rational number.Specifically, when γ is either an integer or a rational number, then from the relation [22] η = 2 it follows that both the real and the imaginary parts of η must be rational numbers.However, from the identity it follows that η ∈ R since γ ∈ R. Consequently, we prove that η is a rational number at γ ∈ Q.
It should be noted that the two-term Machin-like formula of kind (11) that we considered in our paper [22] represents a practical interest.Recently, a group of independent researchers has developed a different method for determination of the value of η at given values of k and γ (see Table 2 in [10]).Unlike our iterative method described in [22], their algorithm is built on the basis of the rational functions R j (n, x) with a remarkable property (see [10] a for detailed description).
Following the same procedure that is described in our work [22], one can also obtain a generalization of identity (9) in the form Consequently, the two-term Machin-like Formula (11) can be generalized as where (compare with Equation (12) above) It should be noted that using the identity [21] arctan(x + y) = arctan(x) + arctan y 1 + (x + y)x the equation ( 14) can be represented as from which it follows that Both constants ϕ and γ in Equation ( 13) may be chosen conveniently.For example, consider a ratio of 22/7, representing a rough approximation of π.Therefore, we can write We can see now that Equation ( 9) is just a specific case ϕ = 2 k−1 of more general form of the two-term Machin-like Formula (13) for π.Therefore, the method described in our paper [22] can also be generalized to generate all two-term Machin-like formulas of kind (14), shown in Table 1 from the work [10] (see first row showing f 28  32 ).In our recent publication [24], using the two-term Machin-like Formula (9) for π, we found the following iterative formula where This equation can be employed to compute digits of the constant π with quadratic convergence (see Mathematica code provided in [24]).Motivated by recent publications [10,23,25,26] in connection to our works [21,22,24,27,28], we developed a new algorithm based on a modified version of the iterative Formula (15).
Although Equation ( 15) provides squared convergence in computing digits of π, its direct application results in a slow convergence rate in the intermediate steps of the calculations of the tangent function.This occurs because the argument of the tangent function in this equation tends to the relatively large value of π/4 as n increases.In this work, we propose a new method that resolves this problem.In particular, we show how Equation ( 15) can be rearranged and used in combination with arctangent terms of the Machin-like formula for π.Such an approach may be promising for efficient computation of π with more than 17 digits per increment of n in Equation ( 29) that will be discussed below.To the best of our knowledge, algorithmic implementation based on the combination of iterative and Machin-like formulas for computing digits of π has never been reported.

Machin-Like Formulas
In our recent work, we derived the following identity [21] where such that a 0 = 0 and a k = √ 2 + a k−1 are nested radicals of 2 and with an initial value B 1,k that can be computed by substituting A k into Equation ( 12) Equation ( 17) implies two important rules.First, since the integer B 0,k is not defined, it follows that at M = 0, the sum of arctangent functions then no further iteration is required, as the fractional part of the number B M +1,k does not exist.
We may compute the coefficient B 1,k by using Equation ( 20) at smaller values of the integer k.However, the computation slows down as k increases.
To resolve this problem, we proposed a more efficient method of computation based on a two-step iterative formula [22] with initial values Equation ( 22), based on the two-step iteration (21), is more efficient for computation of the constant B 1,k than Equation ( 20), since at larger values of the integer k, the rapidly growing exponent 2 k−1 in Equation ( 20) drastically decelerates the computation.
The case k = 4 requires more computations.In particular, we can see that and at M = 0, we have Consequently, we can write π 4 = 2 Repeating the same procedure over and over again up to M = 5, we can finally obtain the following seven-term Machin-like formula where is also an integer.Therefore, all iterations are completed.
From these examples, we can see that Hermann's ( 23), Machin's (1) and the derived (24) formulas for π belong to the same generic group, since all of them can be constructed from their generalized form (17) at different integers k, equal to 2, 3 and 4, respectively.Further, we will use Equation (24) as an example for computing digits of π.
The following Mathematica code: validates Equation ( 24) by returning True.This code applies the product relation (6) for verification.

Tangent Function
Since the tangent function can be represented as a series expansion where B 2n are the Bernoulli numbers, defined by the contour integral Consequently, the accuracy of the tangent function improves with decreasing the argument x.
Unfortunately, the series expansion of the tangent function ( 25) requires the determination of the Bernoulli numbers, which is itself a big challenge [29][30][31][32].For example, one of the most known formulas for computation of the Bernoulli numbers is based on double summation with the binomial coefficients that decelerate the computation.
There are other equations for computation of the tangent function [31][32][33] and one of the efficient techniques to perform computation of the tangent function is the Newton-Raphson iteration (see [27] for more details).In particular, the following iteration formula with an initial value that can be taken as can be employed.This iteration leads to tan(x) = lim n→∞ s n .
Iteration (27) leads to the quadratic convergence of the tangent function.In practice, however, quadratic convergence can be achieved only if a sufficiently large number of the summation terms are applied in the series expansions like (3), (30) or (31) to approximate the arctangent function.We have already applied the Newton-Raphson iteration (27) to compute the digits of π [27].
As an option, we can also apply the most common equation , where sine and cosine functions in the numerator and denominator are represented as the Maclaurin expansions.Unfortunately, this representation is not optimal for practical application since it requires separate computations of the expansion terms for the sine and cosine functions.However, if we rewrite the tangent function as and take into account that sin(2x) = 2 sin(x) cos(x), we can obtain the following identity Although this this representation of the tangent function is not common, its application is significantly advantageous, since each nth term in the expansion can be utilized again just by multiplying 2 2n+1 to obtain the expansion for the denominator Such a technique may accelerate computation and reduce memory usage.Thus, we can write the following iterative procedure according to Equation (28).
In this work, we used truncated Equation ( 29) as an alternative to Equations (25) and ( 27) since one of our objectives in this work is to develop an algorithm for computing the digits of π as simply as possible and without undesirable surd numbers.

Arctangent Function
Since in this work we utilize the terms from the Machin-like Formula ( 17), a series expansion with rapid convergence of the arctangent function should be applied.Based on our empirical results, we can consider two equations with rapid convergence that can be used for implementation.The first equation is Euler's series expansion [6,34] The second equation is given by the series expansion [35] arctan where the expansion coefficients are computed by a two-step iteration such that κ 1 (x, t) = 1/(xt), The derivation of the series expansion (31), shown in our work [35], is based on the Enhanced Midpoint Integration (EMI) formula (see also [22]) where we imply that The series expansion (31) appears to be rapid in convergence.In particular, even at M = 1, its convergence rate is faster by many orders of the magnitude than that of Euler's Formula (30) (see [35] for details).Thus, by taking M = 1, Equation ( 31) can be conveniently rearranged as where the expansion coefficient can be computed by the following two-step iteration Series expansions (30) and ( 32) are significantly faster in convergence than the Maclaurin series expansion (3).In this work, the truncated series expansion (32) is used.Although a value of integer M that is greater than 1 further improves the convergence rate, it may be preferable to apply Equation (32) rather than its generalization (31), since an increment of M just by 1 increases the number of the terms in series expansion (31) by a factor of 2.

Modified Iteration
Changing the variable θ n → 1/σ n in Equation ( 15) leads to a more convenient form where Consequently, the constant π can be found by iteration in accordance with Equation ( 16) such that π Comparing the limit ( 34) with (see [28] for derivation) we can obtain an important relation Consequently, the argument of the tangent function in Equation ( 33) tends to 1 with increasing n.Since, in Equation ( 33), the argument of the tangent function tends to 1, the convergence per iteration of n in applied Equation ( 29) is expected to be extremely slow.Indeed, if the argument of the tangent function is not small enough, then, in accordance with relation (26), we cannot gain a reasonable convergence of the tangent function.
This problem can be effectively resolved by introducing a constant such that c ≈ The value of this constant is close to σ n when n is sufficiently large (see Equation ( 35)).Therefore, with the help of this constant, we can modify the iteration Formula (33) by using the following procedure The expression for σ n can be simplified as follows or or where α = tan 2 k−1 c .Since |δ n | is supposed to be a small value, we may expect a reasonable convergence rate in iteration (36) at each increment n in Equation ( 29).The tangent function in Equation ( 36) is represented twice.Therefore, we can define τ n = tan 2 k−1 δ n and rewrite this equation in a more simplified form

Methodology
In our algorithm, we utilize a modified Equation (37).This yields the iteration based on a set of equations We can take a few initial terms in the Machin-like Formula ( 17) as a value to compute the constant c.
Return to Equation ( 24), corresponding to the case k = 4.If we take only the first term, then the constant first term of Equation ( 24) = arctan 1 10 (39) is not close enough to π/(4 × 2 k−1 ) to compute the digits of π.As a consequence, we cannot achieve a rapid convergence.Specifically, our empirical results show that with constant (39), the set (38) of iteration formulas provides convergence of four to five digits of π per increment of n in Equation (29).However, if we take the first two terms from Equation ( 24) such that first term of Equation ( 24) − arctan 1 84 second term of Equation ( 24) then the convergence rate significantly improves, providing ten digits of π per increment.Consider the following identity [36] tan This identity can be used for computation of the constant α.Thus, according to identity (41), we have Once the exact value of the constant α is calculated according to Equation (43), we can use the set of iteration Formula (38) for computing the digits of π.The Mathematica codes and their description are provided in the next section.
Although the convergence rate of 10 digits of π per increment of n in Equation ( 29) is relatively high, application of the Machin-like Formula (17) may not be efficient at a smaller value of the integer k.In particular, the computation of both arctangents arctan(1/10) and arctan(1/84) is expected to be slow, since the two integers in the argument denominators (10 and 84) are not big enough for rapid convergence in accordance with Equation (4).
Consider another case where k is sufficiently large.We can take, for example, k = 27.Applying identity (17) together with Equation ( 22), we can find that [22] Recently, the same equation was derived by Gasull et al. by a different method of computation (see Table 2 in [10] showing the row with f 522,185,807 522,185,816 ).At first glance, it may appear problematic to find the corresponding coefficient α, since at k = 27, the substitution 2 27−1 = 67108864 into identity (41) results in an expression that is impossible to compute due to the extremely large value of the exponent.However, application of the two-term Machin-like formula for π of kind (11) gives a big advantage, since the multiplier 2 k−1 is continuously divisible by 2. Thus, using the elementary trigonometric identity in the iteration process, we can compute the required constant α.Specifically, in accordance with this identity and due to relation x = tan (arctan(x)), at k = 27, the following iteration formula with initial value in which the integer 85445659 = ⌊A 27 ⌋ is calculated according to Equation (18).It is interesting to note that the number of digits (522,185,816) in this equation is the same as the number of digits in the denominator of Equation (44).At a larger value of the integer k, it is sufficient to take only the first term (leading term) of the Machin-like formula for π to achieve rapid convergence.In particular, using only the leading term from Equation (44) 67108864 arctan 1 85445659 , the set of iteration Formula (38) provides 17 to 18 digits of π per increment of n in Equation ( 29).The corresponding Mathematica code and its description are discussed in the next section.

Mathematica Codes and Description
The Mathematica codes consist of seven cells that can be copied and pasted directly to the Mathematica notebook.The first cell is given by the following code: that defines the two-step iterative method for the arctangent function, based on Equation (32).This function can be invoked by running the command atanF.
The second cell provides the code: defines the output format that includes the heading and ending parts for intermediate computed data.The header and the ending parts can be invoked by running the commands heading and ending, respectively.The code in the fourth cell:  computes the coefficient α according to Equation (41) and coefficient c at k = 4 by taking only the leading (first) term of the Machin-like Formula (24) for π.As we do not require the highest precision at each cycle of computation within the while loop, the precision increases with increasing n in Equation (29).The parameter of the precision is given as 5 + 5n, where n is the increment in Equation ( 29), since at n = 1, the number of correct digits of π is 5 and multiplier 5 is the largest number of digits of π per increment n.
Suppose we know the first 100 digits of π.These digits of π can be extracted from Mathematica by using the following division ⌊10 100 π⌋ 10 100 .
Once we know 100 digits of π, we can use the described iteration to double its number.Running this cell generates the following output:  ------------------------------200 digits of π after iteration that shows the first five and last fifteen cycles.As we can see, each increment of n by one gives four to five digits of π.After 39 cycles, the convergence slows down due to saturation.In particular, after 40th cycle, the number of digits doubles to 200.This saturation occurs since we have reached the limit of squared convergence in the determination of the digits of π.
The code in the fifth cell:  computes the coefficients c and α at k = 4 with the help of Equations ( 40), ( 45) and (46).Again, the precision increases with increasing n in Equation (29).The parameter of the precision is given as 12 + 10n, where n is the increment in Equation ( 29), since at n = 1, the number of correct digits of π is 12 and multiplier 10 is the largest number of digits of π per increment n.We can use the 200 already obtained digits of π to double it by iteration.
Running this cell produces the output:  -------------------------------402 digits of π after iteration that shows first five and last fifteen cycles.The convergence rate always remains the same, 10 digits per increment of n.After 39 cycles, the convergence slows down.In particular, after the 40th cycle, the number of digits of π doubles to 402.Further cycles do not contribute to increasing the number of digits since again we reached the limit in the determination of the digits of π.
The code in the sixth cell is the most interesting:  as it provides an excellent convergence rate at 17 to 18 digits per increment n in Equation ( 29) at k = 27.Although the exact value of the rational number α can be used (see Equation ( 46)), this code does not utilize it since a regular desktop or laptop computer requires a few hours for computation by using Equation (46).However, we can observe a rapid convergence computing this coefficient with tangent and arctangent approximations based on Equations (46) and (32), respectively.The code in this cell utilizes only a leading term of the Machin-like Formula (44) for π.This is possible to achieve since at larger value of k = 27, the argument 1/85445659 in the leading arctangent term becomes small.Consequently, this reduces the value of δ n and improves the convergence rate according to Equation (29).
The precision of this code increases with increasing n in Equation (29).The parameter of the precision is given as 25 + 18n, where n is the increment in Equation ( 29), since at n = 1, the number of correct digits of π is 25 and the multiplier 18 is the largest number of digits of π per increment n.Again, we can use the 402 already obtained digits of π to double it by iteration.
Since the proposed method requires at least the leading term of the Machin-like Formula (44) for π, we have to verify its convergence rate to ensure efficient computation.Thus, at k = 27, the corresponding value of the leading term in Equation ( 44 showing the convergence rate of the arctangent function.This code generates following output: Increment of n | Correct digits --------------------------------- As we can see, the code provides 16 to 17 correct digits of the arctangent function value (47) per increment of n in Equation (32).Furthermore, computing just a single arctangent term greatly minimizes Lehmer's measure (5).In particular, the value of the measure of the arctangent term (47) is only µ = 1 log 10 (85445659) ≈ 0.126077.

Conclusions
A new iterative method for computing the digits of π by argument reduction of the tangent function is developed.This method combines a modified iterative formula for π with squared convergence and a leading arctangent term from the Machin-like Formula (17).The computational test shows that algorithmic implementation can provide more than 17 digits of π per increment n in Equation ( 29).This method requires no surd numbers, and with an arbitrarily large k, there is no upper limit for the convergence rate.