From Euclid to Gregory |
Proposition 3 of Book VI of Euclid's Elements states that |
If an angle of a triangle be bisected, and the straight line cutting the angle cut the base also, then the segments of the base will have the same ratio as the remaining sides of the triangle... |
A typical example of this elementary proposition is shown below. |
Proposition 3 simply asserts that ae = cd. Euclid presents a characteristically clever proof of this fact, but for the modern reader the proposition is more easily seen in terms of areas. The two triangles on the bases d and e have the same altitudes, so the ratio of their areas is d/e. The same two triangles can also be seen as sharing a common base b, and if we drop perpendiculars from the outer vertices to this [extended] line, we have two similar right triangles with the hypotenuses a and c, so the heights of these triangles above the base b are proportional to a and c respectively. Hence the ratio of areas is a/c, so we have a/c = d/e. (The proposition is even more immediate if we recall that the area of a triangle is the dot product of the edge vectors.) |
Archimedes made use of this simple fact in computing bounds on the value of p. Suppose that for some integer N we know the length 2L0 of an edge of a regular N-gon circumscribing a unit circle, as illustrated below for a square (N = 4). |
Half the perimeter of the N-gon is therefore NL0. With N = 4 we have L0 = 1, and the half perimeter is 4, which is a very crude upper bound on the half-perimeter p of the circle. Archimedes noticed that we can improve this bound given by the perimeter of an N-gon by computing the perimeter of a 2N-gon, in which each edge subtends half the angle subtended by the edges of the N-gon. Thus we merely need to bisect the angle subtended by one of the (half)edges, and then determine the new (half)edge length by means of Euclid's Proposition VI;3. For example, in the figure above, the segment L1 is half the edge of an 8-gon circumscribing the unit circle, and by Euclid's proposition we have |
from which it follows that |
The half-perimeter of the octagon is therefore equal to 3.313..., a somewhat better upper bound on p. Of course, knowing L1, we can also compute the hypotenuse H1 of the right triangle formed by the segment L1 and the radius of the circle by Pythagoras' formula |
Repeating this process, we can compute the half-edge length L2 of a 4N-gon and the corresponding hypotenuse H2, and so on. At each stage an upper bound for p is given by 2kNLk, which clearly converges on p as k increases. (Historically, Archimedes began his calculation with a hexagon, i.e., with N = 6, for which and . Also, he chose to work with the reciprocals of the edge lengths.) |
Naturally we can give a more succinct description of the above facts in terms of trigonometric functions, and we can include both upper and lower bounds. For example, letting 2x denote the angle subtended by the original edge in the above figure, we see that sin(2x)=L0/H0, cos(2x)=1/H0, tan(2x)=L0, and sin(x)=L1/H1, cos(x)=1/H1, tan(x)=L1, so equation (1) is just an expression of the familiar identity |
The right-most expression shows why it was convenient for Archimedes to work with the reciprocal functions cot(x)=1/tan(x) and csc(x)=1/sin(x), because then this identity is simply |
We can also divide through equation (2) by H1 and replace the resulting expressions with their trigonometric equivalents to give |
We can, of course, combine equations (3) and (4) (for successive steps) into a single recursive equation involving the cotangent function |
This formula can also be derived directly from the tangent addition formula |
which can be solved to give tan(x) in terms of tan(2x) as |
This is essentially equivalent to (5), except that it is expressed in terms of the tangent rather than the cotangent. If we write the basic tangent addition formula as a quadratic in tan(x), and divide through by tan(2x) and tan(x)2, we get the equivalent equation |
Solving for cot(x) in terms of cos(2x) gives equation (5). |
Just as the half-perimeter of a regular N-gon circumscribing a unit circle is N tan(pi/N), it's easy to see that the half-perimeter of the inscribed N-gon is N sin(pi/N), so we can compute lower bounds on p simply by finding the values of sin(x) or, equivalently, csc(x), for each x. Equation (3) provides the means of doing this, given the values of cot(2x) and cot(x). |
Archimedes actually performed his entire calculation twice, once using a lower bound on , and once using an upper bound. The result for the inscribed polygons with the lower bound on gave him an overall lower bound on p, and the result for circumscribed polygons with the upper bound on gave the overall upper bound on p. Of course, each step of the algorithm involves the extraction of another square root, so in each case he selected rational approximations that we lower or upper bounds, as appropriate. As with his bounds on , he did not explain how he arrived at these rational approximations, presumably because the method was well-known. In most cases he was taking the square root of the sum of a large square and a small square, so he could easily begin with the root of the large square and then adjust using something like the binomial expansion. |
One disadvantage in the method used by Archimedes to compute bounds on p is that the values of the cotangent and cosecant roughly double on each step. Alternatively, if we use the tangent and sine functions, the terms are reduced by half on each step. These are both undesirable qualities, because the terms are added to a fixed quantity on each step, and there is a tendency to lose significant information when adding two quantities of very different magnitude. In other words, the algorithm is somewhat ill-conditioned from a computational standpoint. To avoid this problem, we would like to find a recursion that gives a sequence of values converging on p itself, rather than growing (or shrinking) geometrically. In other words, we would like to find a recursive relation for two sequences Tj and Sj such that |
Note the factors of 2 in Sj+1 and Tj+1. This means that these values are the half-perimeters of their respective polygons, and need not be multiplied or divided by 2n at the end of the computation. Now, we already have part of a recurrence on these values, because substituting into equation (3) gives |
which just expresses the fact that 2tan(a) is the harmonic mean of tan(2a) and sin(2a). The arbitrary factor k cancels out, because (3) is homogeneous. Of course, in terms of the reciprocal functions tj = 1/Tj and sj = 1/Sj this can be written as the arithmetic mean |
To find the companion equation for Sj+1 in terms of Tj and Sj we need another homogeneous trigonometric identity involving the sine and tangent, so that the factor k will again cancel out. Combining the two half-angle formulas |
we have |
Replacing the tangent and sine functions with T and S functions gives |
which can also be expressed in terms of the reciprocals as |
Equations (6) and (7) enable us to compute successively better approximations of p (or 1/p) simply by repeatedly evaluating harmonic (or arithmetic) and geometric means. For example, beginning with T0 = 4 and S0 = , we have the sequences shown in the attached Table. |
Interestingly, the recurrence given by (6b) and (7b) is essentially identical to the recurrence that James Gregory (1638-1675) found as a means of evaluating the areas or arc lengths of sections of a circle or a right hyperbola (the latter being based on the hyperbolic trigonometric functions, which have relations directly analogous to the circular functions). His algorithm is sometimes expressed in the form |
With the initial values a0 = 2, b0 = 4, these two sequences converge on p. Needless to say, this is the same as the S and T sequences defined above, with the indices of the S sequence shifted by -1, so that we have aj = Sj-1 and bj = Tj. |
Interestingly, the mixed indices on the right hand side of this recurrence has a very significant effect on the convergence properties of the sequence. If we make both of the right hand sides simply functions of the prior terms (and revert back to the reciprocal form) then we have the familiar arithmetic-geometric mean |
which converges quadratically. (Note that the geometric mean is really just the arithmetic mean in the logarithmic domain, so these two functions are, in a sense, natural duals of each other. Also, it's interesting that the geometric mean has the same form if the arguments are replaced by their reciprocals, whereas the arithmetic mean is transposed with the harmonic mean.) Letting M(a0,b0) denote the limiting value of this recurrence with the initial values a0,b0, Gauss discovered (in 1799) that |
as well as the useful definite elliptic integral |
Somewhat surprisingly, it wasn't until 1976 (Salamin and Brent) that the highly convergent properties of the arithmetic-geometric mean were applied to the computation of the number p. They deduced the relation |
which is the basis for algorithms that double the number of correct digits on each iteration. |