Tangent to PI
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