Next Article in Journal
Basic Circuit Model of Voltage Source Converters: Methodology and Modeling
Previous Article in Journal
From Algebro Geometric Solutions of the Toda Equation to Sato Formulas
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Sanjar M. Abrarov
1,2,3,*,
Rehan Siddiqui
2,3,4,
Rajinder Kumar Jagpal
2,4 and
Brendan M. Quine
1,3,4
1
Thoth Technology Inc., Algonquin Radio Observatory, Achray Rd., RR6, Pembroke, ON K8A 6W7, Canada
2
Epic College of Technology, 5670 McAdam Rd., Mississauga, ON L4Z 1T2, Canada
3
Department Earth and Space Science and Engineering, York University, 4700 Keele St., Toronto, ON M3J 1P3, Canada
4
Department Physics and Astronomy, York University, 4700 Keele St., Toronto, ON M3J 1P3, Canada
*
Author to whom correspondence should be addressed.
AppliedMath 2024, 4(3), 868-888; https://doi.org/10.3390/appliedmath4030047
Submission received: 7 June 2024 / Revised: 5 July 2024 / Accepted: 14 July 2024 / Published: 19 July 2024

Abstract

:
In this work, we consider the properties of the two-term Machin-like formula and develop an algorithm for computing digits of π by using its rational approximation. In this approximation, both terms are constructed by using a representation of 1 / π in the binary form. This approach provides the squared convergence in computing digits of π without any trigonometric functions and surd numbers. The Mathematica codes showing some examples are presented.

1. Preliminaries

In 1876, English astronomer and mathematician John Machin demonstrated an efficient method to compute the digits of π by using his famous discovery [1,2,3,4]
π 4 = 4 arctan 1 5 arctan 1 239 .
In particular, due to the relatively rapid convergence of this formula, he was the first to compute more than 100 digits of π . Nowadays, the equations of the kind
π 4 = j = 1 J A j arctan 1 B j ,
where A j and B j are either integers or rational numbers, are named after him as the Machin-like formulas for π [1,2,3,4].
Theorem 1 below shows the arctangent Formula (2) for π . We can use this equation as a starting point to generate the two-term Machin-like formula for π of the kind
π 4 = 2 k 1 arctan 1 x + arctan 1 y ,
where k , x N and y 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 1 [5]:
π 4 = 2 k 1 arctan 2 a k 1 a k ,
where a 0 = 0 and a k = 2 + a k 1 are nested radicals.
Proof of Theorem 1.
Since
cos π 2 2 = 1 2 2 = 1 2 a 1 ,
from the identity
cos 2 x = 2 cos 2 x 1 ,
it follows, by induction, that
cos π 2 3 = 1 2 2 + 2 = 1 2 a 2 ,
cos π 2 4 = 1 2 2 + 2 + 2 = 1 2 a 3 ,
cos π 2 5 = 1 2 2 + 2 + 2 + 2 = 1 2 a 4 ,
cos π 2 k + 1 = 1 2 2 + 2 + 2 + + 2 k square roots = 1 2 a k .
Therefore, using
sin π 2 k + 1 = 1 cos 2 π 2 k + 1 ,
we obtain
π 2 k + 1 = arctan 1 cos 2 π 2 k + 1 cos π 2 k + 1 = arctan 2 a k 1 a k
and Equation (2) follows.    □
Lemma 1 and its proof below show how Equation (2) is related to the well-known limit (3) for π .
Lemma 1.
There is a limit such that [6]
π = lim k 2 k 2 2 + 2 + 2 + + 2 k 1 square roots .
Proof of Lemma 1. 
Let
lim k a k = x .
Then, we can write
lim k a k = lim k 2 + a k 1 = 2 + lim k a k
or
x = 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 the conclusion that
lim k a k = 2 .
Consequently, this gives
lim k 2 a k 1 = 2 lim k a k 1 = 2 lim k a k = 0 .
Since the formula (2) is valid for any arbitrarily large k, we can write
π 4 = lim k 2 k 1 arctan 2 a k 1 a k .
As arctan x / x 1 when x 0 , we have
π 4 = lim k 2 k 1 2 a k 1 a k = lim k 2 k 1 2 a k 1 / lim k a k = lim k 2 k 2 2 a k 1 = lim k 2 k 1 2 a k
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 π [7].
Lemma 2.
In the following equation
π 4 = 2 k 1 arctan 1 a k / 2 a k 1 + arctan 1 β k ,
the value β k is always a rational number.
Proof of Lemma 2. 
It is convenient to define
α k = a k / 2 a k 1 .
Using the identity
arctan 1 x = 1 2 i ln x + i x i
and taking into consideration that
π 4 = arctan 1 ,
Equation (4) can be represented as
π 2 i = ln α k + i α k i 2 k 1 ln β k + i β k i
or
i = α k + i α k i 2 k 1 β k + i β k i .
It is not difficult to see by substitution that the following formula [7]
β k = 2 α k + i α k i 2 k 1 i i
is a solution of Equation (6). Since k is a positive integer greater than or equal to 1, we can see that the real and imaginary parts of the expression
α k + i α k i 2 k 1 = α k 2 1 α k 2 + 1 + i 2 α k α k 2 + 1 2 k 1
cannot be an irrational number if α k is an integer. This means that the real and imaginary parts of the value β 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 β k is equal to zero. This completes the proof.    □
This is not the only method to generate the two-term Machin-like formulas for π of kind (4). Recently, Gasull et al. proposed a different method to derive the same equation (see Table 2 in [8]). Lemma 3 and its proof below shows how Equation (4) can be represented in a trigonometric form [7].
Lemma 3.
Equation (4) can be expressed as
π 4 = 2 k 1 arctan 1 α k + arctan 1 sin 2 k 1 arctan 2 α k α k 2 1 cos 2 k 1 arctan 2 α k α k 2 1
Proof of Lemma 3. 
Define
κ 1 = α k 2 1 α k 2 + 1
and
λ 1 = 2 α k α k 2 + 1
such that
β k = 2 κ 1 + i λ 1 2 k 1 i i ,
in accordance with Equation (7). Then, using de Moivre’s formula, we can write the complex number in polar form as
κ 1 + i λ 1 2 k 1 = κ 1 2 + λ 1 2 2 k 2 cos 2 k 1 Arg κ 1 + i λ 1 + i sin 2 k 1 Arg κ 1 + i λ 1 .
Thus, substituting this equation into Equation (9), we obtain
β k = cos 2 k 1 Arg κ 1 + i λ 1 1 sin 2 k 1 Arg κ 1 + i λ 1 = cos 2 k 1 Arg α k + i α k i 1 sin 2 k 1 Arg α k + i α k i .
Using the relation
Arg x + i y = arctan y x , x > 0 ,
we can write
Arg κ 1 + i λ 1 = arctan 2 α k α k 2 1 .
Consequently, from the identities (10) and (11), we obtain
β k = cos 2 k 1 arctan 2 α k α k 2 1 1 sin 2 k 1 arctan 2 α k α k 2 1
and Equation (8) follows.    □
It should be noted that computation of the constant β k by using Equation (7) is not optimal. Specifically, at lager values of the integer k, its application slows down the computation. The application of Equation (12) for the computation of the constant β k is also not desirable due to the presence of the trigonometric functions. Theorem 2 and its proof show how this problem can be effectively resolved [7].
Theorem 2. 
The rational number β k is given by
β k = κ k 1 λ k ,
where the coefficients κ k and λ k can be found by a two-step iteration
κ n = κ n 1 2 λ n 1 2 λ n = 2 κ n 1 λ n 1 , n = 2 , 3 , 4 , , k
with initial values
κ 1 = α k 2 1 α k 2 + 1
and
λ 1 = 2 α k α k 2 + 1 .
Proof of Theorem 2. 
We notice that the following power reduction
κ 1 + i λ 1 2 k 1 = κ 1 + i λ 1 2 2 2 k 1 powers of 2 = κ 2 + i λ 2 2 2 2 k 2 powers of 2 = κ 3 + i λ 3 2 2 2 k 3 powers of 2 = κ n + i λ n 2 2 2 k n powers of 2 = κ k 2 + i λ k 2 2 2 = κ k 1 + i λ k 1 2 = κ k + i λ k ,
where the numbers κ n and λ n can be found by two-step iteration (14), leads to
β k = 2 κ k i λ k i i = 2 κ k κ k 2 + λ k 1 2 + i 2 1 λ k κ k 2 + λ k 1 2 1 ,
according to Equation (7).
From the Lemma 2, it follows that the imaginary part of the value β k is equal to zero. Consequently, the equation above can be simplified as
β k = 2 κ k κ k 2 + λ k 1 2 .
However, since
i 2 1 λ k κ k 2 + λ k 1 2 1 = 0
it follows that
2 1 λ k κ k 2 + λ k 1 2 = 1 κ k 2 + λ k 1 2 = 2 1 λ k .
Substituting this result into the denominator of Equation (15), we obtain Equation (13). This completes the proof.    □
There is an interesting relation between β k and nested radicals a k 1 and a k . Specifically, comparing Equation (4) with Equation (18) from Lemma 4 below, one can see that
arctan 1 β k = 2 k 1 arctan a k / 2 a k 1 1 + a k / 2 a k 1 a k / 2 a k 1 ,
where
a k / 2 a k 1 = a k / 2 a k 1 a k / 2 a k 1
denotes the fractional part.
Lemma 4.
lim k arctan 1 β k = 0 .
Proof of Lemma 4. 
It is not difficult to see that using a change in the variable
y y 1 + x + y x
in the trigonometric identity
arctan x + arctan y = arctan x + y 1 x y
leads to
arctan x + y = arctan x + arctan y 1 + x + y x .
We define, for convenience, the fractional part as
γ k = a k / 2 a k 1 .
Then, from Equation (2), it follows that
π 4 = 2 k 1 arctan 1 α k + γ k
Solving the equation
1 α k + γ k = 1 α k + z ,
we obtain
z = γ k α k α k + γ k .
Therefore, we have
π 4 = 2 k 1 arctan 1 α k γ k α k α k + γ k .
Using the identity (16), we can rewrite Equation (17) as
π 4 = 2 k 1 arctan 1 α k arctan γ k 1 + α k α k + γ
or
π 4 = 2 k 1 arctan 1 a k / 2 a k 1 arctan a k / 2 a k 1 1 + a k / 2 a k 1 a k / 2 a k 1 .
As the fractional part
a k / 2 a k 1 < 1
while
lim k 1 + a k / 2 a k 1 a k / 2 a k 1 = ,
we conclude that
lim k 2 k 1 arctan a k / 2 a k 1 1 + a k / 2 a k 1 a k / 2 a k 1 = lim k arctan 1 λ k κ k = lim k arctan 1 β k = 0 .
   □
Lemma 5 and its proof below show how to obtain a single-term rational approximation (34) for π by truncating Equation (19).
Lemma 5.
π 4 = lim k 2 k 1 a k / 2 a k 1 .
Proof of Lemma 5. 
From the Lemma 4, it immediately follows that
π 4 = lim k 2 k 1 arctan 1 a k / 2 a k 1 .
According to Lemma 1
lim k 2 a k 1 = 0
and
lim k a k = 2 .
Therefore, we can infer that
lim k a k / 2 a k 1 =
or
lim k 1 a k / 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], we propose a methodology for determination of the coefficients α k without computing nested square roots of 2 (see Equation (5)) and develop an algorithm providing the squared convergence per iteration in computing the digits of π . To the best of our knowledge, this approach is new and has 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 publication [15], we showed how the two-term Machin-like formula for π represented in form (8) is used to derive an iterative formula
θ n , k = 1 1 θ n 1 , k + 1 2 k 1 tan 2 k 1 θ n 1 , k , k 1 ,
where the initial value can be taken as
θ 1 , k = 2 k
such that
π 4 = 2 k 1 lim n 1 θ n , k .
This iterative formula provides a squared convergence in computing the digits of π (see the Mathematica code in [15]).
Taking the change in the variable θ n 1 / θ n in Equation (21) yields a more convenient form
θ n , k = θ n 1 , k + 2 k 1 tan 2 k 1 θ n 1 , k , k 1 .
Consequently, the limit (22) can be rearranged now as
π 4 = 2 k 1 lim n θ n , k .
Comparing this limit with Equation (2), we can see that this iteration procedure results in
lim n θ n , k = arctan 1 α k
and is used de facto for the determination of the arctangent function. The detailed procedure showing how to implement the computation with a high convergence rate is given in our recent publication [16].
Alternatively, we can transform the two-term into multi-term Machin-like formulas for π consisting of only the integer reciprocals. In order to do this, we can use the following equation [14]
π 4 = 2 k 1 arctan 1 α k + m = 1 M arctan 1 μ m , k + arctan 1 μ M + 1 , k ,
where
μ m , k = 1 + μ m 1 , k μ m 1 , k μ m 1 , k μ m 1 , k
with an initial value
μ 1 , k = β k .
For example, by taking k = 2 in Equation (24), we can construct
π 4 = 2 2 1 arctan 1 α 2 + arctan 1 μ 1 , 2 = 2 arctan 1 2 arctan 1 7 .
This equation is commonly known as the Hermann’s formula for π [17,18].
By taking k = 3 , Equation (24) gives
π 4 = 2 3 1 arctan 1 α 3 arctan 1 μ 1 , 3
representing the original Machin’s formula (1) for π .
The case k = 4 requires more calculations to obtain the multi-term formula for π consisting of only integer reciprocals. In particular, at  M = 1 , we obtain
π 4 = 2 4 1 arctan 1 α 4 + arctan 1 μ 1 , 4 = 8 arctan 1 10 arctan 1758719 147153121
At M = 2 and M = 3 , Equation (24) yields
π 4 = 2 4 1 arctan 1 α 4 + arctan 1 μ 1 , 4 + arctan 1 μ 2 , 4 = 8 arctan 1 10 arctan 1 84 arctan 579275 12362620883
and
π 4 = 2 4 1 arctan 1 α 4 + arctan 1 μ 1 , 4 + arctan 1 μ 2 , 4 + arctan 1 μ 3 , 4 = 8 arctan 1 10 arctan 1 84 arctan 1 21342 arctan 266167 263843055464261 ,
respectively.
Repeating this procedure over and over again, at  M = 5 , we end up with a 7-term Machin-like formula for π consisting of only integer reciprocals
π 4 = 2 4 1 arctan 1 α 4 + m = 1 5 arctan 1 μ m , 4 + arctan 1 μ 6 , 4 = 8 arctan 1 10 arctan 1 84 arctan 1 21342 arctan 1 991268848 arctan 1 193018008592515208050 arctan 1 197967899896401851763240424238758988350338 arctan 1 Ω ,
where
Ω = 11757386816817535293027775284419412676799191500853701 8836932014293678271636885792397 ,
is the largest integer (see the Mathematica code in [16] that validates this seven-term Machin-like formula for π ).
As we can see, Hermann’s (25), Machin’s (1) and the derived (26) formulas for π belong to the same generic group, as all of them can be constructed from the same equation-template (24).
It is not difficult to show that applying the following identity
arctan 1 x = arctan 1 x + 1 + arctan 1 x 2 + x + 1 , x 1 , 0 ,
to Equation (26), we obtain
π 4 = 8 arctan 1 10 arctan 1 84 arctan 1 21342 arctan 1 991268848 arctan 1 193018008592515208050 arctan 1 197967899896401851763240424238758988350338 n = 1 N arctan 1 ω n , N 2 ,
where
ω 1 = Ω + 1 , ω 2 = ω 1 1 2 + ω 1 + 1 , ω 3 = ω 2 1 2 + ω 2 + 1 , ω N 1 = ω N 2 1 2 + ω N 2 + 1 , ω N = ω N 1 1 2 + ω N 1 .
This expanded variation of Equation (26) together with iterative formula (23) can also be used for computing the digits of π (see Section 5 in [16] for more details).
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 [19]
arctan x = n = 0 2 2 n n ! 2 2 n + 1 ! x 2 n + 1 1 + x 2 n + 1 .
The second equation is [7,20]
arctan x = 2 n = 1 1 2 n 1 g n x g n 2 x + h n 2 x ,
where
g n x = g n 1 x 1 4 / x 2 + 4 h n 1 x / x
and
h n x = h n 1 x 1 4 / x 2 4 g n 1 x / x
with initial values
g 1 x = 2 / x
and
h 1 x = 1 .
The computational test we performed shows that Equation (28) is faster in convergence than Equation (27). Recently, we proposed the generalization of the arctangent series expansion (28) [20].

2.2. Tangent Function

Generally, transformation of the two-terms to multi-terms formulas for π with integer reciprocals is not required. In particular, we can use the Newton–Raphson iteration method. For example, both arctangent terms in the two-term Machin-like formula (4) for π can be computed directly by using the following iterative formulas
σ n = σ n 1 1 2 tan σ n 1 / 2 1 + tan 2 σ n 1 / 2 2 2 tan σ n 1 / 2 1 tan 2 σ n 1 / 2 1 α k
and
τ n = τ n 1 1 2 tan τ n 1 / 2 1 + tan 2 τ n 1 / 2 2 2 tan τ n 1 / 2 1 tan 2 τ n 1 / 2 1 β k
with initial values
σ 1 = 1 α k
and
τ 1 = 1 β k
such that
arctan 1 α k = lim n σ n
and
arctan 1 β k = lim n τ n .
Since the convergence of the Newton–Raphson iteration is quadratic [21], with proper implementation of the tangent function, we may achieve an efficient computation.
The tangent function can be expanded as
tan x = n = 1 1 n 1 2 2 n 2 2 n 1 B 2 n 2 n ! x 2 n 1 = x + x 3 3 + 2 x 5 15 + 17 x 7 315 + 62 x 9 2835 + ,
where B 2 n are the Bernoulli numbers that can be defined by the exponential generating function
x e x 1 = n 0 B n x n n ! .
However, application of Equation (29) is not desirable since the computation of the Bernoulli numbers B 2 n itself is a big challenge [22,23,24,25].
In order to resolve this problem, we propose the following limit [16]:
tan x = lim n 2 p n 2 x q n x ,
where
p n x = p n 1 x + r n 1 x ,
q n x = q n 1 x + 2 2 n 1 r n 1 x ,
such that
r n = 1 n 2 n + 1 ! x 2 n + 1 ,
with initial values
p 0 x = 0 ,
q 0 x = 0 .
Specifically, it has been shown that at k = 27 , in the two-term Machin-like formula (4) for π , the application of Equation (30) results in more than 17 digits of π per increment n (see the Mathematica codes in [16]).
Since the multiplier x 2 n + 1 in Equation (30) increases exponentially with increasing integer n, the iteration process rapidly reduces the number of digits in the mantissa when the tangent function with some fixed accuracy is computed. Decreasing the number of digits in the mantissa may be advantageous as the computation requires less memory usage.

3. Algorithmic Implementation

In our recent publication, we showed that [15]
β k 2 1 tan 2 k 1 α k .
Consequently, in accordance with Equation (4), we obtain the following approximation:
π 4 2 k 1 arctan 1 α k + 1 2 1 tan 2 k 1 α k .
At sufficiently large k, the value 1 / α k < < 1 . Therefore, according to Lemma 5, in this equation, we can replace the first arctangent term by a rational number 2 k 1 / α k . This gives
π 4 2 k 1 α k + 1 2 1 tan 2 k 1 α k .
Unfortunately, we cannot compute efficiently the tangent function in this approximation since its argument 2 k 1 / α k is not a small number as it tends to π / 4 with increasing k. However, again by taking into consideration that 1 / α k < < 1 from which it follows that
1 α k arctan 1 α k ,
we can write
π 4 2 k 1 α k + 1 2 1 tan 2 k 1 arctan 1 α k .
Now we can take advantage of the fact that the multiplier 2 k 1 is continuously divisible by 2. Therefore, we can use the trigonometric identity
tan 2 arctan x = 2 x 1 x 2
k 1 times over and over again. Thus, this leads to the following iterative formula:
η n x = 2 η n 1 x 1 η n 1 2 x
with an initial value
η 1 x = 2 x 1 x 2
such that
η k 1 1 α k = tan 2 k 1 arctan 1 α k .
Since the left side of the equation above provides an exact value without (tangent and arctangent) trigonometric functions, we can regard this equation
π 4 2 k 1 α k + 1 2 1 η k 1 1 α k
as a rational approximation. This rational approximation of the two-term Machin-like formula for π can be used in an algorithm providing a quadratic convergence. This can be achieved with help of Theorem 3.
Theorem 3. 
There is a conversion formula
lim k n = 1 k 1 10 n + 1 α n mod 2 = 0.01010001011111001100 ,
such that
0.01010001011111001100 2 = 1 π .
Proof. 
The proof is related to the parity of the integer α k . According to Lemma 5, we can write
π 4 2 k 1 α k
or
1 π α k 4 × 2 k 1 .
Consequently, if the integer α k is even, then
α k 4 × 2 k 1 α k 1 4 × 2 k 2 = 0 2 .
However, if the integer α k is odd, then
α k 4 × 2 k 1 α k 1 4 × 2 k 2 = 0.00 00 k + 1 zeros 1 2 .
This means that α k contributes a binary digit 0.00 00 k + 1 zeros 1 2 to the previous value α k 1 if and only if α k is odd. This completes the proof.    □
Consider an example. There are four consecutive values α 4 = 10 , α 5 = 20 and α 6 = 40 and α 7 = 81 . Since the first three values are even, we have
α 4 4 × 2 4 1 = α 5 4 × 2 5 1 = α 6 4 × 2 6 1 = 0 . 0101 2 .
However, since α 7 = 81 is odd, we obtain
α 7 4 × 2 7 1 = 0 . 0101 2 + 0 . 00000001 2 = 0 . 01010001 2 .
Consider how the number of digits of π can be doubled without computing square roots for the nested radicals a k . We can take, for example, k = 7 . This yields
α 7 = a 7 2 a 7 1 = 2 + 2 + 2 + 2 + 2 + 2 + 2 2 2 + 2 + 2 + 2 + 2 + 2 = 81 .
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 1 / π in the binary form according to Theorem 3. Thus, ignoring the first two initial zeros in the binary output of Equation (33), we have a corresponding sequence
( s ) 1 = 1 , 0 , 1 , 0 , 0 , 0 , 1 , 0 , 1 , 1 , 1 , 1 , 1 , 0 , 0 , .
This sequence can be obtained by using the built-in Mathematica directly using 1 / π . 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 (36):
{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)
( s ) 1 7 = 1 , 0 , 1 , 0 , 0 , 0 , 1
and apply it accordingly as
α k = 2 α k 1 , if k th binary digit is 0 2 α k 1 + 1 , if k th binary digit is 1 .
Explicitly, defining α 0 = 0 , this step-by-step procedure results in
α 1 = 2 α 0 + 1 = 2 × 0 + 1 = 1 , α 2 = 2 α 1 + 0 = 2 × 1 + 0 = 2 , α 3 = 2 α 2 + 1 = 2 × 2 + 1 = 5 , α 4 = 2 α 3 + 0 = 2 × 5 + 0 = 10 , α 5 = 2 α 4 + 0 = 2 × 10 + 0 = 20 , α 6 = 2 α 5 + 0 = 2 × 20 + 0 = 40 , α 7 = 2 α 6 + 1 = 2 × 40 + 1 = 81 .
Thus, we can see how a very simple procedure can be used to determine the value of the rational number α 7 = 81 without using a sophisticated Equation (35) consisting of 14 nested square roots of 2.
At α 7 = 81 , the corresponding Machin-like formula is
π 4 = 2 7 1 arctan 1 α 7 + arctan 1 β 7 = 64 arctan 1 81 arctan 2154947582 4298183679 111 digits 4599489202 6981324801 113 digits ,
where the constant
β 7 = 4599489202 6981324801 113 digits 2154947582 4298183679 111 digits
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 Equation (37) by returning True.
Suppose that we do not know the sequence other than s 1 7 . However, with the help of Equation (32), we can find other digits of π in the iterative process. In particular, using Equation (31), we have
η 7 1 1 81 = 2310519339 5639754240 113 digits 2288969863 1341570561 113 digits 1.00941448647564092749
Substituting this value into Equation (32), we can find a significantly better approximation of π .
The following is the Mathematica code:
Print["Equation (34) at k = 7: ",
  MantissaExponent[N[\[Pi]-4∗(64/81),20]][[2]] // Abs,
    " digits of \[Pi]"];
Print["Equation (32) at k = 7: ",
  MantissaExponent[
    N[\[Pi]-4∗(64/81+1/2∗(1-1.00941448647564092749)),
      20]][[2]]//Abs," digits of \[Pi]"];
which produces the following output:
Equation (34) at k = 7: 1 digits of π
Equation (32) at k = 7: 4 digits of π
The initial sequence s 1 7 helped us to find the value α 7 = 81 . Now, due to the higher accuracy of Equation (32), we can generate the sequence in which its upper index is doubled
s 1 14 = 1 , 0 , 1 , 0 , 0 , 0 , 1 , 0 , 1 , 1 , 1 , 1 , 1 , 0
and with the help of this sequence, we can find the corresponding value α 14 = 10,430.
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 s 1 8 , then we obtain α 16  = 41,722 instead of the correct value α 16 = 41,721. Therefore, the upper index of the sequence should be slightly less than two.
The two-term approximation (32) doubles the number of digits of π as compared to the single-term approximation (34). This means that using the sequence s 1 k , we can obtain all sequences s 1 k + 1 , s 1 k + 2 , s 1 k + 3 , etc., up to s 1 k 0 , where k 0 is an integer slightly smaller than 2 k . Our numerical results show that doubling the value of k does not always provide the correct sequence, as a few binary digits at the end of the sequence s 1 2 k occasionally may not be correct. However, when we use the empirical equation
k 0 = 2 1 32 k ,
then the corresponding sequence s 1 k 0 is a sub-sequence of the infinite sequence (36) 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 π 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 π doubles at each iteration.
It should be noted that the tangent function tan 2 k 1 / α k can also be computed by combining together Equations (30) and (31). In particular, taking σ k 1 , we can apply the identity
tan 2 k 1 α k = tan 2 σ 2 k 1 σ α k = η σ tan 2 k 1 σ α k .
The following is the Mathematica code where the identity (38) is implemented by using Equations (30) and (31) at σ = 100 and k = 7323 :
Clear[tanF]
(∗ Computing tangent ∗)
tanF[n_,x_,k_]:=Module[{p=0,q=0,m=1,x0=N[x,k]},xSq=x0^2;
  While[m<=n,
    r=(-1)^(m-1)/(2∗m-1)!∗x0;
    p=p+r;
    q=q+2^(2∗m-1)∗r;
    m++;x0=x0∗xSq];
    (2∗p^2)/q];
(∗ Case k = k0 = 7323 ∗)
n=1;
\[Sigma]=100;
While[n<=10,numb=tanF[n,2^(k0-1-\[Sigma])/\[Alpha][k0],k0];
  Print["n = ",n,": ",MantissaExponent[\[Pi]-
  4∗(2^(k0-1)/\[Alpha][k0]+
  1/2∗(1-\[Eta][\[Sigma],numb,k0]))][[2]]//Abs,
  " digits of \[Pi]"];n++];
This code returns the following output:
n = 1: 60 digits of π
n = 2: 121 digits of π
n = 3: 182 digits of π
n = 4: 244 digits of π
n = 5: 306 digits of π
n = 6: 368 digits of π
n = 7: 430 digits of π
n = 8: 492 digits of π
n = 9: 554 digits of π
n = 10: 617 digits of π
This output shows a high convergence rate. Specifically, we can see that each increment n in Equation (30) contributes more than 60 digits of π . Generally, this convergence rate has no upper limitation. However, in order to increase σ in Equation (38), we have to increase the value of the integer k.
Although this convergence rate is significantly higher than those in some modern algorithms for computing digits of π [4], even without nested square roots of 2, the proposed method at large values k still requires a powerful computer to compute the coefficients α k by using the binary form of 1 / π .
The presence of the multiplier 2 k 1 σ in Equation (38) is continuously divisible by 2. Therefore, its application may be advantageous for more rapid computation since each multiplication by 2 implies just a binary shift of the mantissa.
The Mathematica code below shows the values of α 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 α k are the same as those reported in [9].
Since the constants α k have been computed already, we can use them to validate the formula (33) for 1 / π 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.010100010112
20  | 0.0101000101111100110002
30  | 0.01010001011111001100000110110112
40  | 0.010100010111110011000001101101110010011102
50  | 0.0101000101111100110000011011011100100111001000100002
————————————————————-
Built-in Mathematica:
1 / π  = [0.010100010111110011000001101101110010011100100010000...]2
according to Equation (33). The original binary representation of the number 1 / π generated by the built-in Mathematica is also shown for comparison (see also [26] for binary sequence for 1 / π ).

4. Conclusions

We consider the properties of the two-term Machin-like formula for π and propose its two-term rational approximation (32). Using this approach, we develop an efficient algorithm for computing digits of π with squared convergence. The constants α k in this approximation are computed without nested square roots of 2.

Author Contributions

S.M.A. developed the methodology, wrote the codes and prepared a draft version of the manuscript. R.S., R.K.J. and B.M.Q. verified, reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article.

Acknowledgments

This work was supported by National Research Council Canada, Thoth Technology Inc., York University and Epic College of Technology. The authors wish to thank the Reviewers for their constructive comments and recommendations.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Beckmann, P. A History of Pi; Golem Press: New York, NY, USA, 1971. [Google Scholar]
  2. Berggren, L.; Borwein, J.; Borwein, P. Pi: A Source Book, 3rd ed.; Springer-Verlag: New York, NY, USA, 2004. [Google Scholar]
  3. Borwein, J.; Bailey, D. Mathematics by Experiment—Plausible Reasoning in the 21st Century, 2nd ed.; Taylor & Francis Group: Abingdon, UK, 2008. [Google Scholar]
  4. Agarwal, R.P.; Agarwal, H.; Sen, S.K. Birth, growth and computation of pi to ten trillion digits. Adv. Differ. Equ. 2013, 2023, 100. [Google Scholar] [CrossRef]
  5. Abrarov, S.M.; Quine, B.M. A formula for pi involving nested radicals. Ramanujan J. 2018, 46, 657–665. [Google Scholar] [CrossRef]
  6. Servi, L.D. Nested square roots of 2. Am. Math. Mon. 2003, 110, 326–330. [Google Scholar] [CrossRef]
  7. Abrarov, S.M.; Quine, B.M. An iteration procedure for a two-term Machin-like formula for pi with small Lehmer’s measure. arXiv 2017, arXiv:1706.08835. [Google Scholar]
  8. Gasull, A.; Luca, F.; Varona, J.L. Three essays on Machin’s type formulas. Indag. Math. 2023, 34, 1373–1396. [Google Scholar] [CrossRef]
  9. Wolfram Cloud. A Wolfram Notebook Playing with Machin-like Formulas. Available online: https://www.wolframcloud.com/obj/exploration/MachinLike.nb (accessed on 5 June 2024).
  10. Campbell, J. Nested radicals obtained via the Wilf–Zeilberger method and related results. Maple Trans. 2023, 3, 16011. [Google Scholar] [CrossRef]
  11. Maritz, M.F. Extracting pi from chaos. Coll. Math. J. 2024, 55, 86–99. [Google Scholar] [CrossRef]
  12. Spíchal, L. Using the golden section to approximate π. Math. Mag. 2022, 97, 315–320. [Google Scholar] [CrossRef]
  13. Alferov, O. A rapidly converging Machin-like formula for π. arXiv 2023, arXiv:arxiv:2403.09654. [Google Scholar]
  14. Abrarov, S.M.; Jagpal, R.K.; Siddiqui, R.; Quine, B.M. A new form of the Machin-like formula for π by iteration with increasing integers. J. Integer Seq. 2022, 25, 22.4.5. [Google Scholar]
  15. Abrarov, S.M.; Jagpal, R.K.; Siddiqui, R.; Quine, B.M. Algorithmic determination of a large integer in the two-term Machin-like formula for pi. Mathematics 2021, 9, 2162. [Google Scholar] [CrossRef]
  16. Abrarov, S.M.; Jagpal, R.K.; Siddiqui, R.; Quine, B.M. An iterative method for computing π by argument reduction of the tangent function. Math. Comput. Appl. 2024, 29, 17. [Google Scholar] [CrossRef]
  17. Borwein, J.M.; Borwein, P.B. Pi and the AGM—A Study in Analytic Number Theory and Computational Complexity; Wiley & Sons Inc.: Hoboken, NJ, USA, 1987. [Google Scholar]
  18. Chien-Lih, H. More Machin-type identities. Math. Gaz. 1997, 81, 120–121. [Google Scholar] [CrossRef]
  19. Chien-Lih, H. An elementary derivation of Euler’s series for the arctangent function. Math. Gaz. 2005, 89, 469–470. [Google Scholar] [CrossRef]
  20. Abrarov, S.M.; Siddiqui, R.; Jagpal, R.K.; Quine, B.M. A generalized series expansion of the arctangent function based on the enhanced midpoint integration. AppliedMath 2023, 3, 395–405. [Google Scholar] [CrossRef]
  21. Ypma, T.J. Historical development of the Newton—Raphson method. SIAM Rev. 1995, 37, 531–551. [Google Scholar] [CrossRef]
  22. Knuth, D.E.; Buckholtz, T.J. Computation of tangent, Euler, and Bernoulli numbers. Math. Comp. 1967, 21, 663–688. [Google Scholar] [CrossRef]
  23. Harvey, D. A multimodular algorithm for computing Bernoulli numbers. Math. Comput. 2010, 79, 2361–2370. [Google Scholar] [CrossRef]
  24. Bailey, D.H.; Bauschke, H.H.; Borwein, P.; Garvan, F.; Vanderwerff, M.T.J.D.; Wolkowicz, H. Computational and Analytical Mathematics; Springer: New York, NY, USA, 2013. [Google Scholar]
  25. Beebe, N.H.F. The Mathematical Function Computation Handbook; Springer International Publishing AG: New York, NY, USA, 2017. [Google Scholar]
  26. The Online Encyclopedia of Integer Sequences. Expansion of 1/Pi in Base 2. OEIS: A127266. Available online: https://oeis.org/A127266 (accessed on 5 June 2024).
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Abrarov, S.M.; Siddiqui, R.; Jagpal, R.K.; Quine, B.M. A Rational Approximation of the Two-Term Machin-like Formula for π. AppliedMath 2024, 4, 868-888. https://doi.org/10.3390/appliedmath4030047

AMA Style

Abrarov SM, Siddiqui R, Jagpal RK, Quine BM. A Rational Approximation of the Two-Term Machin-like Formula for π. AppliedMath. 2024; 4(3):868-888. https://doi.org/10.3390/appliedmath4030047

Chicago/Turabian Style

Abrarov, Sanjar M., Rehan Siddiqui, Rajinder Kumar Jagpal, and Brendan M. Quine. 2024. "A Rational Approximation of the Two-Term Machin-like Formula for π" AppliedMath 4, no. 3: 868-888. https://doi.org/10.3390/appliedmath4030047

APA Style

Abrarov, S. M., Siddiqui, R., Jagpal, R. K., & Quine, B. M. (2024). A Rational Approximation of the Two-Term Machin-like Formula for π. AppliedMath, 4(3), 868-888. https://doi.org/10.3390/appliedmath4030047

Article Metrics

Back to TopTop