I recently noticed a web page devoted to the German mathematician Ludolph Van Ceulen and the number PI. The author of the page includes his own "modest contribution to the search for the Ludolphine", namely, the recursive formula x[n+1] = x[n] - tan(x[n]) , pi/2 < x[0] < 3pi/2 which rapidly converges to the value of PI. I'm not sure if this is an intentional joke, but in any case it's interesting to review why this technique is not really useful, in spite of the fact that you can confirm its "correctness" on your pocket calculator. Of course, if we're permitted the explicit use of transendental functions such as tan(x) in our "search for the Ludolphine", then we could simply write pi = 4 invtan(1) and the "search" is over. On the other hand, if we limit ourselves to rational operations and root extractions, then we'll need to convert our transcendental functions to power series. For example, the formula pi = 4invtan(1) can be expanded into the series pi = 4 (1 - 1/3 + 1/5 - 1/7 +...) and of course we can derive more rapidly converging series in various ways. (For more on this, see Machin's Merit.) However, the recursive operation x -> x - tan(x) doesn't immediately translate into rational operations because it needs to work on the branch of x values between pi/2 and 3pi/2, whereas the series expansion of tan(x) requires x to be in the fundamental branch from -pi/2 to pi/2. Thus, to evaluate tan(x) with pi/2 < x < 3pi/2 using rational functions we must actually evaluate tan(x - pi). This is exactly how an electronic calculator or computer proceeds when asked to evaluate tan(x) for a value of x in this range. But notice that the iteration x -> x - tan(x-pi) will converge on "pi" regardless of the value we assume for "PI". For example, if our calculator was built with the stored value of PI erroneously equal to exactly 3.0, but it correctly computes tan(z) for values of z between -3/2 and +3/2, then the iteration will converge on 3. In other words, using the operation x -> x - tan(x) to "compute" PI is essentially equivalent to simply pressing the "PI" key on your calculator. The iteration procedure doesn't generate any new information about the value of pi beyond whatever value was stored in the computer's "PI" memory at the start. This is most clearly seen by noting that tan(z) ~ z for small arguments z, so the method on the web page is essentially x_new = x_old - (x_old - PI) so it's nothing but a disguised and inefficient way of reading whatever value of PI is stored in your calculator's memory. After posting the above comments on the internet I received email saying that I had under-estimated the utility of the "x-tan(x)" method, because we can compute tan(x) with x in the range PI/2 to 3PI/2 using rational operations by computing the ratio sin(x)/cos(x), noting that the series expansions for sine and cosine converge for all x. I agree that this is possible. On the other hand, I'm doubtful that this is what the author had in mind, because if we're going to use the iteration x -> x - sin(x)/cos(x) as x goes to PI, it's clear that the cos(x) in the denominator is superfluous, because it goes to -1. Thus we are effectively using the iteration x -> x + sin(x). If this was the intent, it seems odd that the method was described as x -> x - tan(x). Of course, even if we use the form x + sin(x), the value of sin(x) for arguments near PI is still computed in most calculators as -sin(x-PI) in order to get rapid convergence of the series, so we're still effectively just reading the value of PI that was already stored in the computer's memory if we use any pre- programmed sine function. On the other hand, the sine series DOES eventually converge for all x, so in principle we could use the iteration x -> x + sin(x) to perform a legitimate computation of PI. Recall that the power series for sin(x) is x^3 x^5 x^7 sin(x) = x - --- + --- - --- + ... 3! 5! 7! and we can evaluate the terms t[0], t[1],.. of this series by setting t[0] = x and then recursively computing t[n] = (-1)^n (x^2) t[n-1] / ((2n)(2n+1)) , n=1,2,3,... For example, the successive terms with x=3 are t[k] sum ------------ ------------ 3.000000000 3.000000000 -4.500000000 -1.500000000 +2.025000000 0.525000000 -0.433928571 0.091071429 +0.054241071 0.145312500 -0.004437906 0.140874594 +0.000256033 0.141130627 -0.000010973 0.141119654 +0.000000363 0.141120017 -0.000000010 0.141120007 This would constitute our first iteration, based on x[0]=3, and it would yield x[1] = 3.141120007. Then we would repeat this with x=x[1], by setting t[0]=x[1] and generating the terms of the sine series t[k] sum ------------ ------------ 3.141120007 3.141120007 -5.165380714 -2.024260707 2.548246281 0.523985574 -0.598633707 -0.074648133 0.082034725 0.007386592 -0.007358243 0.000028349 0.000465392 0.000493741 -0.000021866 0.000471875 0.000000793 0.000472668 -0.000000023 0.000472645 This gives as the result of our second iteration the value x[2] = 3.14159265(2). We could also try to salvage the original method of x -> x - tan(x) by using the addition formula for the tangent function. For example, to evaluate tan(3) we could start by computing tan(3/8) using the power series, and then compute tan(3/4), tan(3/2) and tan(3) by repeated application of the doubling rule tan(2x) = 2 tan(x) / (1 - tan(x)^2) In fact, if we don't mind these addition formula computations, we could ensure extremely fast convergence of the basic tangent series by starting with, say, tan(3/1024), and then doubling the argument ten times using the addition formula. However, depending on the implementation of the arithmetic functions that we are using, the divisions involved in the tangent addition formula may be more work than just evaluating a few more terms of the power series. (Incidentally, the approach of using the tangent addition formula leads directly to another the well-known "low tech" method of computing pi; see Machin's Merit.) Anyway, it's clear that there are several ways of getting a legitimate (though not outstandingly efficient) computation of PI from iterations like x -> x + sin(x) and x -> x - tan(x). This isn't surprising, because these are basically just applications of Newton's method to numerically find the roots of sin(x)=0 and tan(x)=0 respectively (noting that the slopes of these functions are +-1 at their roots). Even better, we can iterate on x -> x + cos(x) to solve cos(x)=0 for the root at PI/2. This has the advantage of working with a smaller argument, so convergence of the power series would be faster. Starting with x[0] = 3/2 we would compute x[1] = 1.570737202 and then x[2] = 1.570796327, which equals PI/2 to about 11 decimal places, and each of these two iterations would require (slightly) fewer terms than in the case of the sine formula. Of these three formulas, the best seems to be x -> x + cos(x), whereas the worst seems to be x -> x - tan(x), because it either contains a superfluous evaluation of the cosine in the denominator (using the sin/cos approach), or it reduces to the sine formula x -> x + sin(x), or it requires the use of the tangent addition formula or some other means to allow us to evaluate (without using a pre-computed value of PI) the tangent function in the range from PI/2 to 3PI/2. Nevertheless, this got me to thinking about other possible ways of deducing the value of pi from rational operations based on the tangent (rather than the arctangent) function. Recall that the power series expansion of tan(x) is tan(x) = x + (1/3)x^3 + (2/15)x^5 + (17/315)x^7 + ... and it might seem that this expansion can't do us much good because the special points on this curve related to pi occur at x values that are multiuples of pi. Thus without knowing pi we can't evaluate this series at any of these special points. We can, however, exploit this series in another way. We know that the series "blows up" when x equals pi/2 but it's finite for every positive x less than pi/2. This suggests that the coefficients of the series must, in some sense, grow at a rate corresponding to (pi/2)^n. As an example, consider the term (1/3)(x^3). The cube roots of the inverse of the coefficient of x^3 in the expansion of tan(x) is 3^(1/3) = 1.4422, which gives a very rough approximation of pi/2. A slightly better approximation is given by the 5th root of the inverse coefficient of x^5, i.e., we have (15/2)^(1/5) = 1.4962. Similarly, if C_25 denotes the coefficient of x^25 in the expansion of tan(x), then we find that 2(1/C_25)^(1/25) = 3.1415926... This observations leads us to some interesting relations between Bernoulli numbers and powers of PI, as explained below. To make this a little more formal, let's begin by noting that the expansion of tan(x) can be generated by dividing the simple expansion of sin(x) by the expansion of cos(x), which gives tan(x) = c[0] x + c[1] x^3 + c[2] x^5 + ... (1) where the coefficients are given recursively by c[0]/0! = 1/1! c[1]/0! - c[0]/2! = 1/3! (2) c[2]/0! - c[1]/2! + c[0]/4! = 1/5! c[3]/0! - c[2]/2! + c[1]/4! - c[0]/6! = 1/7! As x goes to PI/2 the value of tan(x) goes to infinity, so for any positive constant w there is an index n such that c[n] x^(2n+1) ~ w, from which it follows that / w \ 1/(2n+1) x ~ ( ----- ) (3) \ c[n]/ As n goes to infinity x goes to PI/2, so we have / w \ 1/(2n+1) PI lim ( ----- ) = ---- (4) n->inf \ c[n]/ 2 Clearly as n increases, the (2n+1)th root of the constant w goes to 1, so we also have the asymptotic relation / 1 \ 1/(2n+1) PI lim ( ----- ) = ---- (5) n->inf \ c[n]/ 2 This enables us to compute PI/2 simply by generating coefficients of the tangent series using the recursive formulas defined above. Of course, the coefficients for the tangent series are known to be expressible in terms of the Bernoulli numbers, i.e., we have 1 2 17 tan(x) = x + - x^3 + -- x^5 + --- x^7 + ... (6) 3 15 315 and in general the coefficient of x^(2n+1) is 2^(2n+2) [2^(2n+2) - 1] |B_{2n+2}| c[n] = ---------------------------------- (7) (2n+2)! Inserting this into the previous limiting relation between c[n] and PI gives the result / (2n+2)! \1/(2n+1) lim PI = ( ---------------------------- ) (8) n->inf \ 2 [2^(2n+2) - 1] |B_{2n+2}| / Of course, we could also use this relation in reverse, to give an approximation for the Bernoulli numbers in terms of PI. However, if we solve for B_{2n+2} in terms of PI^(2n+1) we will no longer be taking the 1/(2n+1)th root, so the constant factor of "w" becomes important. Recall that w was defined as any fixed positive constant, but now we need to determine its optimum value in order to match the Bernoulli numbers. Solving (4) for w, and inserting from (7) for the value of c[n], we have the following expression for the value of w needed to make (4) an exact equality for a given index n 2^(2n+2) - 1 w[n] = 2 ------------ |B_{2n+2}| PI^(2n+1) (9) (2n+2)! The values of w[n] for the first few n are listed below: n w[n] w[n] PI --- --------------------------------- -------- 0 PI / 2 = 1.570796 4.934801 1 PI^3 / 24 = 1.291928 4.058711 2 PI^5 / 240 = 1.275082 4.005788 3 2 PI^7 / 4725 = 1.278430 4.016306 4 4 PI^9 / 93555 = 1.274505 4.003975 5 2764 PI^11/ 638512875 = 1.273552 4.000984 The right hand column lists the product w[n] PI, and clearly suggests that the optimum value of w[n] is 4/PI (although I don't see how to prove this). So, inserting this value into (9) and solving for B_{2n+2}, we have the approximate formula for the Bernoulli numbers 2 (2n+2)! |B_{2n+2}| ~ ------------------------ (10) [2^(2n+2) - 1] PI^(2n+2) If we define the index k = n+1 and neglect the -1 in the first term of the denominator, the expression becomes 2 (2k)! |B_{2k}| ~ ---------- (11) (2PI)^(2k) Sure enough, this corresponds to lower limit of the the well-known bounds 2(2k)! 2(2k)! ------------------------ > |B_{2k}| > ---------- (12) [1 - 2^(1-2k)](2PI)^(2k) (2PI)^(2k) as listed in Abramowitz and Stegun. These bounds are quite good. For example, the absolute value of B_12 is 691/2730, and the bounds given by (12) are 691.1674801 / 2730 and 690.829996 / 2730, so it already gives a fairly accurate estimate even at this low value of k. Moreover, if we go back to equation (10) and DON'T neglect the -1 term in the denominator, we get |B_12| ~ 690.9986969 / 2730 which is extremely precise, and very nearly the average of the upper and lower bounds. (Actually, the average has the factor 4095/4094 rather than 4096/4095, and gives the marginally more accurate estimate of 690.9987381 / 2730.)

Return to MathPages Main Menu