Unconditional Applicability of Lehmer’s Measure to the Two-Term Machin-like Formula for π

Sanjar M. Abrarov
Rehan Siddiqui
Rajinder K. Jagpal
Brendan M. Quine

Lehmer defined a measure

μ=(∑_( j=1)^( J))(1/((Subscript[log, 10])(|(Subscript[β, j])|))),

where the (Subscript[β, j]) may be either integers or rational numbers in a Machin-like formula for π. When the (Subscript[β, j]) are integers, Lehmer’s measure can be used to determine the computational efficiency of the given Machin-like formula for π. However, because the computations are complicated, it is unclear if Lehmer’s measure applies when one or more of the (Subscript[β, j]) are rational. In this article, we develop a new algorithm for a two-term Machin-like formula for π as an example of the unconditional applicability of Lehmer’s measure. This approach does not involve any irrational numbers and may allow calculating π rapidly by the Newton–Raphson iteration method for the tangent function.

1. Introduction

In 1706, the English astronomer and mathematician John Machin discovered a two-term formula for π [1-3]

(π/4)=4arctan((1/5))-arctan((1/239)),

that was later named in his honor. This formula for π appeared to be more efficient than any others known by that time. In particular, due to the relatively rapid convergence of the right-hand side of (1), he was able to calculate 100 decimal digits of π [1]. Nowadays, identities of the form

(π/4)=(∑_(j=1)^J)(Subscript[α, j])arctan((1/(Subscript[β, j]))),

where (Subscript[α, j]) and (Subscript[β, j]) are either integers or rational numbers, are called the Machin-like formulas for π. Consequently, a two-term Machin-like formula for π is given by

(π/4)=(Subscript[α, 1])arctan((1/(Subscript[β, 1])))+(Subscript[α, 2])arctan((1/(Subscript[β, 2]))).

If in equation (3) the constants (Subscript[α, 1]) and (Subscript[β, 1]) are some positive integers and (Subscript[α, 2])=1, then the unknown value (Subscript[β, 2]) can be found as [4, 5]

(1/(Subscript[β, 2]))=(2/((((((Subscript[β, 1])+i)/((Subscript[β, 1])-i)))^(Subscript[α, 1]))+i))+i⇔(Subscript[β, 2])=(2/((((((Subscript[β, 1])+i)/((Subscript[β, 1])-i)))^(Subscript[α, 1]))-i))-i.

Furthermore, since we assumed that (Subscript[α, 1]) and (Subscript[β, 1]) are positive integers, from equation (4) it immediately follows that (Subscript[β, 2]) must be either an integer or a rational number.

In 1938, Lehmer [6] introduced a measure (see also [7])

μ=(∑_(j=1)^J)(1/((Subscript[log, 10])(|(Subscript[β, j])|))),

showing how much computational effort is required for a specific Machin-like formula for π. In particular, when μ is small, then less computational effort is required and, consequently, the computational efficiency of this formula is higher. Lehmer’s measure is smaller if there are fewer summation terms J and the constants (Subscript[β, j]) are larger in magnitude. For more efficient computation, the constants (Subscript[β, j]) should be larger by absolute value, since it is easier to approximate the arctangent function as its argument tends to zero (see [7] for more details).

It is also important to emphasize that in the same paper [6] Lehmer presented a few Machin-like formulas where some of the (Subscript[β, j]) are not integers but rational numbers. This signifies that Lehmer assumed that his measure (5) remains valid whether the (Subscript[β, j]) are integers or rational numbers.

In 2002 Kanada [8], using the following self-checking pair of Machin-like formulas for π,

(π/4)=44arctan((1/57))+7arctan((1/239))-12arctan((1/682))+24arctan((1/12943))

and

(π/4)=12arctan((1/49))+32arctan((1/57))-5arctan((1/239))+12arctan((1/110443)),

computed more than 1 trillion digits of π. These two examples show that the Machin-like formulas have colossal potential in the computation of decimal digits of π.

In 1997, Chien-Lih showed a remarkable formula [9]

(π/4)=183arctan((1/239))+32arctan((1/1023))-68arctan((1/5832))+12arctan((1/110443))-12arctan((1/4841182))-100arctan((1/6826318))

with μ≈1.51244. According to Weisstein [10], this Lehmer’s measure is the smallest known value for the (Subscript[β, j]) consisting of integers only. Later Chien-Lih [11] showed how Lehmer’s measure can be reduced even further by using an Euler type of identity in an iteration for generating the two-term Machin-like formulas like (3) such that (Subscript[β, 1]) and (Subscript[β, 2]) are rational numbers.

In [12], we derived the following simple identity (see also [13])

(π/4)=(2^(k-1))arctan(((Sqrt[ 2-(Subscript[c, k-1])])/(Subscript[c, k]))),	k=1,2,3,⋯,

where (Subscript[c, 0])=0 and (Subscript[c, k])=(Sqrt[2+(Subscript[c, k-1])]), and described how using this identity, another efficient method for generating the two-term Machin-like formula for π

(π/4)=(2^(k-1))arctan((1/(Subscript[β, 1])))+arctan((1/(Subscript[β, 2]))),

with small Lehmer’s measure can be developed. In this approach, the constant (Subscript[β, 1]) can be chosen as a positive integer such that

(Subscript[β, 1])=⌊((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))⌋,

and, in accordance with equation (4), the constant (Subscript[β, 2]) in equation (7) can be found from

(Subscript[β, 2])=(2/((([((Subscript[β, 1])+i)/((Subscript[β, 1])-i)])^(2^(k-1)))-i))-i.

It is not reasonable to solve equation (9) directly to determine the rational number (Subscript[β, 2]), as its solution becomes tremendously difficult with increasing integer k. However, this problem can be effectively resolved by using a very simple two-step iteration procedure discussed in the next section. Therefore, our approach in generating the two-term Machin-like formula (7) for π with small Lehmer’s measure is much easier than Chien-Lih’s method [11].

Wetherfield [7] provides a detailed explanation clarifying the significance of Lehmer’s measure that shows how much computation is required for a given Machin-like formula for π when all the constants (Subscript[β, j])∈ℤ. However, it is unclear if this paradigm is also applicable when at least one number (Subscript[β, j]) is rational. More specifically, the problem that occurs in computing the two-term Machin-like formula for π in Chien-Lih’s [11] and our [4] iteration methods is related to the rapidly growing number of digits in the numerators and denominators of (Subscript[β, 1]) or (Subscript[β, 2]). This occurs simultaneously with an attempt to decrease Lehmer’s measure. As a result, the subsequent exponentiation in conventional algorithms makes computing the decimal digits of π inefficient. Therefore, the applicability of Lehmer’s measure for a Machin-like formula for π for the case (Subscript[β, j])∈ℚ is questionable. For example, Lehmer’s measure may be small, say less than 1. This means that less computational work is needed to calculate π. However, due to the large number of digits in the numerators and denominators in (Subscript[β, 1]) or (Subscript[β, 2]), a computer performs more intense arithmetic operations that make the runtime significantly longer. Consequently, we ask, Is Lehmer’s measure still applicable when at least one constant from the set (Subscript[β, j]) is not an integer but a rational number?

Motivated by an interesting paper in regard to equation (7) that was recently published [14], we further develop our previous work [4, 5]. In this article, we propose a new algorithm showing how unconditional applicability of Lehmer’s measure for the two-term Machin-like formula (7) for π can be achieved. We also describe how linear and quadratic convergence to π can be implemented.

2. Preliminaries

As mentioned, the number of summation terms J in equation (2) should be reduced in order to minimize Lehmer’s measure. Since at J=1 there is only one Machin-like formula for π,

(π/4)=arctan(1),

we consider the case J=2. We attempted to find a method that can be used to generate a two-term Machin-like formula for π with small Lehmer’s measure. Equation (7) provides an efficient way to do this.

In fact, the original Machin formula (1) for π appears quite naturally from the equations (7), (8) and (9) at k=3. Specifically, equation (8) provides

(Subscript[β, 1])=⌊((Subscript[c, 3])/(Sqrt[2-(Subscript[c, 3-1])]))⌋=⌊((Sqrt[2+(Sqrt[2+(Sqrt[2])])])/(Sqrt[2-(Sqrt[2+(Sqrt[2])])]))⌋=5.

Substituting into equation (9) results in

(1/(Subscript[β, 2]))=(2/((([(5+i)/(5-i)])^(2^(3-1)))+i))+i=-(1/239).

Consequently, at k=3 we get the constants (2^(k-1))=4, (Subscript[β, 1])=5 and (Subscript[β, 2])=-239 in equation (7). Since the arctangent function is odd (arctan(-x)=-arctan(x)), the constants for equation (3) can be rearranged as (Subscript[α, 1])=4, (Subscript[α, 2])=-1, (Subscript[β, 1])=5 and (Subscript[β, 2])=239. This corresponds to the original Machin formula (1) for π.

Theorem 2.1.

There are only four possible cases for a two-term Machin-like formula (3) for π when all four constants (Subscript[α, 1]), (Subscript[α, 2]), (Subscript[β, 1]) and (Subscript[β, 2]) are integers:

Text@Grid[{
{"",(Subscript[α, 1]),(Subscript[α, 2]),(Subscript[β, 1]),(Subscript[β, 2])},
{"Machin",4,-1,5,239},
{"Euler",1,1,2,3},
{"Hermann",2,-1,2,7},
{"Hutton",2,1,3,7}
},Frame->All,Alignment->Right]

The proof for Theorem 2.1 can be found in [10].

Lemma 2.2

If in equation (7) (Subscript[β, 1])∈ℤ, then (Subscript[β, 2])∉ℤ at any integer k>3.

Proof

As we can see from the four cases given in Theorem 2.1, the largest possible value for (Subscript[α, 1])=(2^(k-1)) is 4 and it occurs at k=3 (see the example above). Therefore, for any integer (Subscript[β, 1]) at integer k>3, the constant (Subscript[β, 2]) in the equation (7) cannot be an integer. □

Theorem 2.3.

(Underscript[lim , k->∞])(Subscript[c, k])=(Underscript[lim , k->∞])(Underscript[(Underscript[(Sqrt[2+(Sqrt[2+(Sqrt[2+⋯+(Sqrt[2])])])]), ⏟]), k square roots])=2.

Proof

This is the simplest kind of Ramanujan nested radical and its proof is straightforward. Let X=(Underscript[lim , k->∞])(Subscript[c, k]). Then

X=(Underscript[lim , k->∞]) (Sqrt[2+(Subscript[c, k-1])])=(Sqrt[2+(Underscript[lim , k->∞])(Subscript[c, k-1])])=(Sqrt[2+(Underscript[lim , k->∞])(Subscript[c, k])]).

Squaring both sides,

(X^2)=2+(Underscript[lim , k->∞])(Subscript[c, k])⇔(X^2)=2+X

and solving results in two possible solutions, (Subscript[X, 1])=-1 and (Subscript[X, 2])=2. Since for any positive index k the value (Subscript[c, k]) is always positive, we exclude the solution (Subscript[X, 1])=-1. □

Lemma 2.4.

(Underscript[lim , k->∞])(Sqrt[2-(Subscript[c, k-1])])=(Underscript[lim , k->∞])(Sqrt[2-(Underscript[(Underscript[(Sqrt[2+(Sqrt[2+⋯+(Sqrt[2])])]), ⏟]), k-1 square roots])])=0.

Proof

The proof follows immediately from Theorem 2.3 since

(Underscript[lim , k->∞])(Sqrt[2-(Subscript[c, k-1])])=(Sqrt[2-(Underscript[lim , k->∞])(Subscript[c, k-1])])=(Sqrt[2-(Underscript[lim , k->∞])(Subscript[c, k])])=(Sqrt[2-2])=0. □

Lemma 2.5.

(Underscript[lim , k->∞])(((((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))))/(Subscript[β, 1]))=1.

Proof

Using equation (8) that defines (Subscript[β, 1]) by the floor function, the limit can be rewritten as

(Underscript[lim , k->∞])(((((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))))/(⌊((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))⌋))=1.

By definition, the fractional part given by the difference

frac(((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])])))=(((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])])))-⌊((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))⌋

is positive and cannot be greater than 1. Therefore, the limit (11) can be rewritten in the form

(Underscript[lim , k->∞])(((((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))))/(⌊((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))⌋))=(Underscript[lim , k->∞])((⌊((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))⌋+frac(((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))))/(⌊((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))⌋))=(Underscript[lim , k->∞])((⌊((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))⌋)/(⌊((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))⌋))+0=1. □

Lemma 2.6

We have that (Underscript[lim , k->∞])|(Subscript[β, 2])|=∞.

Proof

Equation (9) is too hard to work with directly, so we start with the limit (11) from Lemma 2.5. Since the limit is 1, the limit of the ratio of the reciprocals of its numerator and denominator must also be 1,

(Underscript[lim , k->∞])((((((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))))^(-1))/(1/(Subscript[β, 1])))=(Underscript[lim , k->∞])(((((Sqrt[2-(Subscript[c, k-1])])/(Subscript[c, k]))))/(1/(Subscript[β, 1])))=1.

From Theorem 2.3 and Lemma 2.4,

(Underscript[lim , k->∞])((Sqrt[2-(Subscript[c, k-1])])/(Subscript[c, k]))=((Sqrt[2-(Underscript[lim , k->∞])(Subscript[c, k-1])])/((Underscript[lim , k->∞])(Subscript[c, k])))=(0/2)=0

and

(Underscript[lim , k->∞])(1/(Subscript[β, 1]))=(Underscript[lim , k->∞])(1/(⌊((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))⌋))=(1/(⌊(((Underscript[lim , k->∞])(Subscript[c, k]))/(Sqrt[2-(Underscript[lim , k->∞])(Subscript[c, k-1])]))⌋))=(1/∞)=0.

Since both the numerator and denominator in (13) tend to zero as k tends to infinity and since arctan(x)->x as x->0, we can rewrite (13) as

(Underscript[lim , k->∞])((arctan(((Sqrt[2-(Subscript[c, k-1])])/(Subscript[c, k]))))/(arctan((1/(Subscript[β, 1])))))=1

or

(Underscript[lim , k->∞])(((2^(k-1))arctan(((Sqrt[2-(Subscript[c, k-1])])/(Subscript[c, k]))))/((2^(k-1))arctan((1/(Subscript[β, 1])))))=1.

Since equation (6) is valid for an arbitrarily large integer k, (15) implies that

(Underscript[lim , k->∞])(2^(k-1))arctan(((Sqrt[2-(Subscript[c, k-1])])/(Subscript[c, k])))=(Underscript[lim , k->∞])(2^(k-1))arctan((1/(Subscript[β, 1])))=(π/4).

However, we also have

(Underscript[lim , k->∞])((2^(k-1))arctan((1/(Subscript[β, 1])))+arctan((1/(Subscript[β, 2]))))=(π/4),

since equation (7) is also valid at an arbitrarily large integer k. Comparing the limits (16) and (17), we get

(Underscript[lim , k->∞])(2^(k-1))arctan((1/(Subscript[β, 1])))=(Underscript[lim , k->∞])[(2^(k-1))arctan((1/(Subscript[β, 1])))+arctan((1/(Subscript[β, 2])))]=(π/4).

Equation (18) is valid if and only if

(Underscript[lim , k->∞])arctan((1/(Subscript[β, 2])))=0,

so (Underscript[lim , k->∞])|(Subscript[β, 2])|=∞, which is (12). □

Lemma 2.7.

Lehmer’s measure (5) may be vanishingly small.

Proof

The limit (14) implies that

(Underscript[lim , k->∞])(Subscript[β, 1])=∞.

From (12) and (19) we conclude that as k->∞, both (Subscript[β, 1]) and (Subscript[β, 2]) tend to infinity. Thus, according to equation (5), Lehmer’s measure μ for the two-term Machin-like formula (7) for π tends to zero with increasing k. □

Theorem 2.8

If equation (8) holds, then the constant (Subscript[β, 2]) is always negative.

Proof

There is only one single-term Machin-like formula (10) for π such that the constants (Subscript[α, 1]) and (Subscript[β, 1]) are both integers. Therefore, in equation (6), at any integer k≥2 the argument of the arctangent function (Sqrt[2-(Subscript[c, k-1])])/(Subscript[c, k]) cannot be represented as the reciprocal of an integer; that is, (Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]) is not an integer. Therefore,

frac(((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])])))>0

and

((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))>⌊((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))⌋, k=2,3,4,….

This implies

arctan(((Sqrt[2-(Subscript[c, k-1])])/(Subscript[c, k])))<arctan(1/⌊((Subscript[c, k])/(Sqrt[2-(Subscript[c, k-1])]))⌋),

or

arctan(((Sqrt[2-(Subscript[c, k-1])])/(Subscript[c, k])))<arctan((1/(Subscript[β, 1]))),

or

(2^(k-1))arctan(((Sqrt[2-(Subscript[c, k-1])])/(Subscript[c, k])))<(2^(k-1))arctan((1/(Subscript[β, 1]))).

We make (20) into an equality by adding a negative error term ε such that

(2^(k-1))arctan(((Sqrt[2-(Subscript[c, k-1])])/(Subscript[c, k])))=(2^(k-1))arctan((1/(Subscript[β, 1])))+ε.

Defining the constant (Subscript[β, 2]) in accordance with equation (9), we find the error term to be

ε=arctan((1/(Subscript[β, 2]))).

Since ε is negative, the constant (Subscript[β, 2]) is also negative. □

3. Iteration Methods

3.1 Arctangent Function

Since in equation (7) the constant (Subscript[β, 1]) is an integer, the first arctangent function term can be computed by any existing method. For example, we can use Euler’s formula for the arctangent function

arctan(x)=(∑_(m=0)^∞)(((2^(2m))(((m!))^2))/((2m+1)!))((x^(2m+1))/(((1+(x^2)))^(m+1))).

Chien-Lih used this formula to develop his iteration method for generating the two-term Machin-like formula for π [11] and later he found an elegant derivation of this formula [15].

We derived another series expansion of the arctangent function [5, 12]:

arctan(x)=i(∑_(m=1)^∞)(1/(2m-1))((1/(((1+2i/x))^(2m-1)))-(1/(((1-2i/x))^(2m-1)))).

It interesting that generalizing the derivation method that was used to get equation (23), we can find by induction the identity

arctan(x)=(∑_(m=1)^M)arctan((( M x)/((M^2)+(m-1)m (x^2)))),

where M is the order of arctan function expansion; that yields simple approximations like

arctan(x)≈(∑_(m=1)^M)(( M x)/((M^2)+(m-1)m (x^2)))

and

arctan(x)≈(∑_(m=1)^M)((( M x)/((M^2)+(m-1)m (x^2)))-(1/3)((((( M x)/((M^2)+(m-1)m (x^2)))))^3)) ,

since arctan(x)=x-(x^3)/3+O ((x^5)).

The representation (23) of the arctangent function is not optimal for algorithmic implementation, since it deals with complex numbers. Fortunately, as we showed in [4], this series expansion can be significantly simplified to

arctan(x)=2(∑_(m=1)^∞)(1/(2m-1))(((Subscript[g, m])(x))/((Subscript[g, m]^2)(x)+(Subscript[h, m]^2)(x))),

where the expansion coefficients are computed by iteration:

(Subscript[g, 1])(x)=2/x, (Subscript[h, 1])(x)=1,
(Subscript[g, m])(x)=(Subscript[g, m-1])(x)(1-4/(x^2))+4(Subscript[h, m-1])(x)/x,
(Subscript[h, m])(x)=(Subscript[h, m-1])(x)(1-4/(x^2))-4(Subscript[g, m-1])(x)/x.

Both series expansions (22) and (24) converge rapidly and need no undesirable irrational numbers to compute π. However, the computational test we performed shows that the series expansion (24) converges more rapidly by many orders of magnitude than Euler’s formula (22) (see Figures 2 and 3 in [4]). Therefore, the series expansion (24) is more advantageous and can be taken to compute the first arctangent function term from the two-term Machin-like formula (7) for π.

The second arctangent function term in equation (7) should not be computed by straightforward substitution of the constant (Subscript[β, 2]) into equation (24). As mentioned, computing with a ratio of numbers with many digits should be avoided. Instead, the second arctangent function term in equation (7) can be computed by Newton–Raphson iteration.

3.2 Rational (Subscript[β, 2])

Once the value of the integer k is chosen, it is not difficult to determine the integer (Subscript[β, 1]) by using equation (8) with the help of Mathematica. However, determining the second constant (Subscript[β, 2]) with (9) is very hard with increasing k, as already mentioned. To overcome this, we proposed a different method [4]. We define the very simple two-step iteration for n=2,3,4,…,

{({{(Subscript[u, n])=(Subscript[u, n-1]^2)-(Subscript[v, n-1]^2)}, {(Subscript[v, n])=2(Subscript[u, n-1])(Subscript[v, n-1])}})

where

(Subscript[u, 1])=(((Subscript[β, 1]^2)-1)/((Subscript[β, 1]^2)+1)),

and

(Subscript[v, 1])=((2(Subscript[β, 1]))/((Subscript[β, 1]^2)+1)).

Then

(Subscript[β, 2])=((2(Subscript[u, k]))/((Subscript[u, k]^2)+((((Subscript[v, k])-1))^2))).

For k>3, the second arctangent term in (7) deals only with a rational number (Subscript[β, 2]). As k increases, the number of digits in the numerator and denominator of the constant (Subscript[β, 2]) increases. For example, at k=6, equation (8) yields

(Subscript[β, 1])=⌊((Sqrt[2+(Sqrt[2+(Sqrt[2+(Sqrt[2+(Sqrt[2+(Sqrt[2])])])])])])/(Sqrt[2-(Sqrt[2+(Sqrt[2+(Sqrt[2+(Sqrt[2+(Sqrt[2])])])])])]))⌋=40

and using (25) we find

(Subscript[β, 2])=-(2634699316100146880926635665506082395762836079845121/38035138859000075702655846657186322249216830232319)≈-69.27013796024857670135….

Consequently, the two-term Machin-like formula (7) for π is generated as

(π/4)=32arctan((1/40))-arctan((38035138859000075702655846657186322249216830232319/2634699316100146880926635665506082395762836079845121)).

Lehmer’s measure for this two-term Machin-like formula for π is μ≈1.16751.

However, if we take k=27, then

(Subscript[β, 1])=⌊((Overscript[(Overscript[(Sqrt[2+(Sqrt[2+(Sqrt[2+…+(Sqrt[2])])])]), ⏞]), 27square roots])/(Sqrt[2-(Underscript[(Underscript[(Sqrt[2+(Sqrt[2+…+(Sqrt[2])])]), ⏟]), 26 square roots])]))⌋=85445659

and using (25) we get

(Subscript[β, 2])=-((Overscript[(Overscript[2368557598 … 9903554561, ⏞]), 522,185,816 digits])/(Underscript[(Underscript[9732933578 … 4975692799, ⏟]), 522,185,807 digits]))≈-2.43354953523904089818…×(10^8).

The corresponding two-term Machin like formula for π is

(π/4)=67108864arctan((1/85445659))-arctan(((Overscript[(Overscript[9732933578 … 4975692799, ⏞]), 522,185,807 digits])/(Underscript[(Underscript[2368557598 … 9903554561, ⏟]), 522,185,816 digits]))),

for which Lehmer’s measure is only μ≈0.245319. Such a large number of digits in the numerator and denominator in the second arctangent function may look unusual. However, some formulas for π obtained from the Borwein integrals involving the sinc function can also result in ratios of integers with a large number of digits. For example, Bäsel and Baillie reported a formula for π that uses a quotient with 453,130,145 digits in the numerator and 453,237,170 digits in the denominator [16]; you can download a file with all the digits of the constant (Subscript[β, 2]) from [17].

As we can see from these examples, Lehmer’s measure decreases with increasing k. However, that occurs simultaneously with a rapid increase in the number of digits in the numerator and denominator of the constant (Subscript[β, 2]). As a result, taking powers of such fractions becomes very slow. This raises the doubt that Lehmer’s measure (5) is indeed relevant for a given Machin-like formula for π when at least one coefficient (Subscript[β, j]) is rational.

To resolve this problem, we considered applying the Newton–Raphson iteration [18]. Specifically, we showed that each consecutive iteration doubles the number of correct digits in the second term of the arctangent function in equation (7). This method is based on the iteration formula:

(Subscript[y, p+1])=(Subscript[y, p])-(1-(sin^2)((Subscript[y, p])))(tan((Subscript[y, p]))-(1/(Subscript[β, 2]))),

such that

(Underscript[lim , p->∞])(Subscript[y, p])=arctan((1/(Subscript[β, 2]))).

The most important advantage of (26) is that the rational number 1/(Subscript[β, 2]) is not involved in the computation of the trigonometric functions. As we can see from (26), 1/(Subscript[β, 2]) is no longer problematic because it is not taken to a power, nor is it the argument of sine or tangent, which consumes most of the runtime. Instead, it is only applied in a single subtraction. This single subtraction (that can be implemented by changing the precision) takes a negligibly small amount of time as compared to the time to compute sin((Subscript[y, p])) or tan((Subscript[y, p])).

To reduce the number of trigonometric functions from two to one, it is convenient to put equation (26) into the form

(Subscript[y, p+1])=(Subscript[y, p])-(1-(((((2tan(((Subscript[y, p])/2)))/(1+(tan^2)(((Subscript[y, p])/2))))))^2))(((2tan(((Subscript[y, p])/2)))/(1-(tan^2)(((Subscript[y, p])/2))))-(1/(Subscript[β, 2]))),

using the elementary trigonometric identities

sin(x)=((2tan(x/2))/(1+(tan^2)(x/2)))

and

tan(x)=((2tan(x/2))/(1-(tan^2)(x/2))).

The tangent function can be found, for example, by using the equation

tan(x)=((sin(x))/(cos(x)))=(((∑_(n=0)^∞)((((-1))^n)/((2n+1)!))(x^(2n+1)))/((∑_(n=0)^∞)((((-1))^n)/((2n)!))(x^(2n)))),

representing the ratio based on the Maclaurin series expansions for the sine and cosine functions. Alternatively, we can use the series expansion

tan(x)=(∑_(n=1)^∞)(((((-1))^(n-1))(2^(2n))((2^(2n))-1)(Subscript[B, 2n]))/((2n)!))(x^(2n-1)),

where (Subscript[B, 2n]) is a Bernoulli number. There are several other ways to compute the tangent, like continued fractions [19, 20] or Newton–Raphson iteration again [21]. Perhaps the argument reduction method for the tangent function can also be used to improve accuracy, but we did not implement that, to keep the algorithm as simple as possible.

3.3 Tangent Function

Equation (27) is based on the Newton–Raphson iteration method. However, this iteration formula contains the tangent function. We showed [18] that once the first arctangent term is computed, the number of correct digits in π can be doubled at each consecutive step of the iteration as with the Newton–Raphson iteration method. To approximate the tangent function with high accuracy, we can apply the Newton–Raphson iteration again [21]. The derivation of the iteration-based equation for the tangent function is not difficult. Let γ=arctan(z).

Using the Newton–Raphson iteration formula

(Subscript[z, q+1])=(Subscript[z, q])-((f((Subscript[z, q])))/(f'((Subscript[z, q])))),

where

f(z)=arctan(z)-γ and f'(z)=(d/(d z))(arctan(z)-γ)=(1/(1+(z^2))),

we get

(Subscript[z, q+1])=(Subscript[z, q])-(1+(Subscript[z, q]^2))(arctan((Subscript[z, q]))-γ)

such that

tan(γ)=(Underscript[lim , q->∞])(Subscript[z, q]).

Substituting γ=(Subscript[y, p])/2 into (28) yields

(Subscript[z, q+1])=(Subscript[z, q])-(1+(Subscript[z, q]^2))(arctan((Subscript[z, q]))-((Subscript[y, p])/2)).

The iteration-based expansion of the series (24) for the arctangent function converges very rapidly. Therefore, we can apply it to compute the arctangent function in equation (29).

4. Implementation

4.1 Linear Convergence

The series expansion (24) of the arctangent function can be computed as follows.

atan[M_,x_]:=atan[M,x]=2 Sum[1/(2 m-1)g[m,x]/(g[m,x]^2+h[m,x]^2),{m,M}]

g[1,x_]:=g[1,x]=2/x;
g[m_,x_]:=g[m,x]=g[m-1,x](1-4/x^2)+4/x h[m-1,x]

h[1,_]=1;
h[m_,x_]:=h[m,x]=h[m-1,x](1-4/x^2)-4/x g[m-1,x]

Next we define the nested radicals consisting of square roots of 2.

c[0]=0;
c[k_]:=Sqrt[2+c[k-1]]

This computes the constants (Subscript[β, 1]) and (Subscript[β, 2]) for the two-term Machin-like formula (7) for π at k=6.

β1=Floor[c[6]/Sqrt[2-c[6-1]]]

For the constant (Subscript[β, 2]), we use the iteration-based formula (25) instead of equation (9).

u[1]=(β1^2-1)/(β1^2+1)

u[n_]:=u[n-1]^2-v[n-1]^2

v[1]=2 β1/(β1^2+1)

v[n_]:=2 u[n-1] v[n-1]

β2=2 u[6]/(u[6]^2+(v[6]-1)^2)

Define a function for the Lehmer’s measure corresponding to a two-term Machin-like formula (7) for π.

LehmerMeasure[β1_,β2_]:=(1/(Log10[Abs[β1]]))+(1/(Log10[Abs[β2]]))

Here is Lehmer’s measure for the case k=6.

LehmerMeasure[β1,β2]//N

The accuracy improves with each iteration, so we do not need to use the highest accuracy at each step. At k=6, the Newton–Raphson iteration-based formula (29) gives 4 to 5 correct digits of the tangent function at each step. Therefore, it is reasonable to use the argument 5 q+2 in SetPrecision, where 2 is taken to minimize rounding or truncation errors.

Recall that (Underscript[lim , p->∞])(Subscript[y, p])=arctan((1/(Subscript[β, 2]))). As an initial guess for the Newton–Raphson iteration, we choose arctan((1/(Subscript[β, 2]))), since this is close to the actual value of the tangent function tan ((Subscript[y, p])/2).

This sets up the recurrence for (Subscript[z, q]) as defined in (29).

z[1,_]=SetPrecision[ArcTan[1/β2],5]

z[q_,x_]:=z[q,x]=SetPrecision[
z[q-1,x]-(1+z[q-1,x]^2)(atan[q,z[q-1,x]]-x),
5 q+2
]

This part of the program, which computes the tangent function, takes most of the runtime. It uses the Newton–Raphson iteration built on the basis of the series expansion (24) of the arctangent function (see (29)).

This sets up the Newton–Raphson iteration formula (27).

y[1]=1/β2

y[p_]:=y[p]=With[{
a=2 z[p,y[p-1]/2],
b=z[p,y[p-1]/2]^2
},
SetPrecision[y[p-1],5 p+2]-(1-(a/(1+b))^2) (a/(1-b)-1/β2)
]

The next part of the program invokes the value tan ((Subscript[y, p])/2) and performs just a few arithmetic operations. As mentioned, the numerator and denominator of 1/(Subscript[β, 2]) contain many digits, but they are not involved in computing the tangent function. Only a minor part of the time is needed to subtract this number (see equation (27)). Consequently, Lehmer’s measure applies unconditionally.

The table shows how the values of the arctangent function gain digits at each iteration step.

Table[{p-1,y[p]},{p,2,11}]//TableForm

This shows the linear rate of convergence to π in terms of the number of accurate digits.

Table[Abs[MantissaExponent[N[π-4 (2^(6-1)atan[1000,1/β1]+y[p]),75]][[2]]],{p,2,11}]

ListLinePlot@%

After the first two iterations, each of the following iterations adds five correct digits to the approximation of π.

Differences@%%

The convergence rate increases as Lehmer’s measure decreases, which can be readily confirmed by increasing k and readjusting the SetPrecision parameter in this algorithm. No undesirable irrational numbers are needed. Furthermore, since Lehmer’s measure can be made vanishingly small, there is no upper bound in the convergence rate per iteration.

4.2 Quadratic Convergence

Consider a variation of the algorithm based on the Newton–Raphson iterations for the tangent function that can be implemented to get quadratic convergence to π. Assume that only d=50 decimal digits of π are known at the beginning.

d=50;

Y=Floor[SetPrecision[Pi/4-32 atan[Ceiling[2 d/4],1/40],2 d]10^d]/10^d

We determined experimentally that at k=6, equation (24) provides four or five correct digits of π at each successive step. Therefore, it is sufficient to take d/4 terms. However, to exclude truncation and rounding errors, we use  d/2 terms of the arctangent function.

We define acc[q] to determine the accuracy of computation and for the number of summation terms in approximation (24) of the arctangent function.

acc[q_]:=If[4 2^q<2d,4 2^q,2d+10]

The multiplier 4×(2^q) was found experimentally. When 4×(2^q)>2d, we restrict its rapid growth to 2d+10, as we do not need extra accuracy at this stage. We add 10 to exclude truncation and rounding errors. We have taken the initial value z=−0.007.

We use equation (29).

Z[1]:=-0.007;
Z[q_]:=Z[q]=SetPrecision[
Z[q-1]-(1+Z[q-1]^2)(atan[Ceiling[acc[q-1]/4],Z[q-1]]-Y/2),
acc[q]
]

For our choice of d, here is the required index.

q=1;While[acc[q]<2 d,q++];q++;q

Table[Z[q],{q,7}]//Column

Once the tangent function is computed, we substitute it into equation (27). The value Y2 has double the accuracy of Y with only one iteration.

Y2=Y-(1-(2 Z[q]/(1+Z[q]^2))^2)(2 Z[q]/(1-Z[q]^2)-1/β2)

The approximate values of the arctangent function by iteration and by the built-in function ArcTan match to 100 places.

N[Y2,2 d]
N[ArcTan[1/β2],2d]

Since we got Y2 with double the accuracy, it can now be used to compute π with significantly improved accuracy.

pi2=4(32atan[Ceiling[4d/4],1/40]+Y2)

This shows that the number d2 of correct decimal digits of pi2 doubled from 50 to 101.

d2=Abs[MantissaExponent[N[Pi-pi2,10000]][[2]]]

More directly, this shows the complete match between the computed approximation of π and that provided by Mathematica.

SetPrecision[pi2,d2]

N[π,d2]

The first arctangent function in the two-term Machin-like formula (7) for π can also be found by using the same algorithm based on the Newton–Raphson iteration. Consequently, this method results in quadratic convergence to π. However, unlike the Brent–Salamin algorithm (also known as the Gauss–Brent–Salamin algorithm) with quadratic convergence to π [2], our approach does not involve any irrational numbers. The number of summation terms in equation (24) and the number of iteration cycles for computation of the tangent function (29) decrease with increasing k. This can be confirmed by using the code given here. To the best of our knowledge, this is the first algorithm showing the feasibility of quadratic convergence to π without using any irrational numbers.

5. Conclusion

In this article, we presented a new algorithm to compute the two-term Machin-like formula (7) for π and showed an example where the condition (Subscript[β, j])∈ℤ was not necessary in order to validate Lehmer’s measure (5). Since this algorithmic implementation lets us avoid subsequent exponentiation of the second constant (Subscript[β, 2]), this approach may be promising for computing π more rapidly without using irrational numbers.

Acknowledgments

This work is supported by National Research Council Canada, Thoth Technology, Inc., York University, Epic College of Technology and Epic Climate Green (ECG) Inc. The authors wish to thank the reviewer for useful comments and recommendations. Constructive suggestions from the editor that improved the content of this work are greatly appreciated.

References

P. Beckmann, A History of Pi, Boulder, CO: Golem Press, 1971.

L. Berggren, J. Borwein and P. Borwein, Pi: A Source Book, 3rd ed., New York: Springer, 2004.

J. Borwein and D. Bailey, Mathematics by Experiment. Plausible Reasoning in the 21st Century, 2nd ed., Wellesley, MA: AK Peters, 2008.

S. M. Abrarov and B. M. Quine, “An Iteration Procedure for a Two-Term Machin-like Formula for Pi with Small Lehmer’s Measure.” arxiv.org/abs/1706.08835.

S. M. Abrarov and B. M. Quine, “The Two-Term Machin-like Formula for Pi with Small Arguments of the Arctangent Function.” arxiv.org/abs/1704.02875.

D. H. Lehmer, “On Arccotangent Relations for π,” The American Mathematical Monthly, 45(10), 1938, pp. 657–664. doi.org/10.2307/2302434.

M. Wetherfield, “The Enhancement of Machin’s Formula by Todd’s Process,” The Mathematical Gazette, 80(488), 1996, pp. 333–344. doi.org/10.2307/3619567.

J. S. Calcut, “Gaussian Integers and Arctangent Identities for π,” The American Mathematical Monthly, 116(6), 2009, pp. 515–530. www.jstor.org/stable/40391144.

H. Chien-Lih, “More Machin-type Identities,” The Mathematical Gazette, 81(490), 1997, pp. 120–121. doi.org/10.2307/3618793.

E. W. Weisstein. “Machin-like Formulas” from Wolfram MathWorld—A Wolfram Web Resource. mathworld.wolfram.com/Machin-LikeFormulas.html.

H. Chien-Lih, “Some Observations on the Method of Arctangents for the Calculation of π,” The Mathematical Gazette, 88(512), 2004, pp. 270–278. www.jstor.org/stable/3620848.

S. M. Abrarov and B. M. Quine, “A Formula for Pi Involving Nested Radicals,” The Ramanujan Journal, 46, 2018, pp. 657–665. doi.org/10.1007/s11139-018-9996-8.

OEIS. “Decimal Expansion of Pi,” A000796. oeis.org/A000796.

Wolfram Cloud. “A Wolfram Notebook Playing with Machin-like Formulas.” develop.open.wolframcloud.com/objects/exploration/MachinLike.nb.

H. Chien-Lih, “An Elementary Derivation of Euler’s Series for the Arctangent Function,” 89(516), 2005, pp. 469–470. doi.org/10.1017/S0025557200178404.

U. Bäsel and R. Baillie, “Sinc Integrals and Tiny Numbers.” arxiv.org/abs/1510.03200.

S. M. Abrarov and B. M. Quine, “The Rational Number (Subscript[u, 2]) for the Two-Term Machin-like Formula for Pi Computed by Iteration.” yorkspace.library.yorku.ca/xmlui/handle/10315/33173.

S. M. Abrarov and B. M. Quine, “Efficient Computation of Pi by the Newton–Raphson Iteration and a Two-Term Machin-like Formula,” International Journal of Mathematics and Computer Science, 13(2), 2018, pp. 157–169. ijmcs.future-in-tech.net/13.2/R-Abrarov.pdf.

J. Havil, The Irrationals: A Story of the Numbers You Can’t Count On, Princeton, NJ: Princeton University Press, 2012.

M. Trott. “Continued Fraction Approximations of the Tangent Function” from the Wolfram Demonstrations Project—A Wolfram Web Resource. demonstrations.wolfram.com/ContinuedFractionApproximationsOfTheTangentFunction.

J.-M. Muller, Elementary Functions. Algorithms and Implementation, 3rd ed., Boston: Birkhäuser, 2016.