Pi - Archimedes (2024)

It is clear from Part 1 that inorder to calculate π we are going to need a better method thanevaluating Gregory'sSeries.

Here is one which was originally discovered by Archimedes in 200BC. Mostof the proofs for series for calculating π require calculus or otheradvanced math, but since Archimedes didn't know any of those you canprove it without! In fact I thought I had discovered it for the firsttime when I was 13 or so, but alas I was 2000 years too late!

Archimedes knew that the ratio of the circumference of a circle to itsdiameter was π. He had the idea that he could draw a regular polygoninscribed within the circle to approximate π, and the more sides he drewon the polygon, the better approximation to π he would get. In modernmaths speak, we would say as the number of sides of the polygon tends toinfinity, the circumference tends to 2π.

Instead of drawing the polygon though, Archimedes calculated the lengthusing a geometric argument, like this.

Draw a circle of radius 1 unit centered at A. Inscribe an N sidedpolygon within it. Our estimate for π is half the circumference of thepolygon (circumference of a circle is 2πr, r = 1, giving 2π). As thesides of the polygon get smaller and smaller the circumference getscloser and closer to 2π.

The diagram shows one segment of the polygon ACE. The side of thepolygon CE has length d_n. Assuming we know d for an N sided polygon,If we can find an expression for the length CD = d_2n, the edge lengthof a polygon with 2N sides, in terms of d_n only, then we can improveour estimate for π. Let's try to do that.

We bisect the triangle CAE to make CAD and DAE the two new equalsegments of the 2N sided polygon.

Using Pythagoras's theorem on the right angled triangle ABC, we cansee:

$$\begin{align}AB^2 + BC^2 &= 1^2 = 1 \\AB &= \sqrt{1 - BC^2} \\\end{align}$$

Given that

$$BC = \frac{d_n}{2}$$

Substituting gives

$$\begin{align}AB &= \sqrt{1-\left(\frac{d_n}{2}\right)^2} \\BD &= 1 - AB \\ &= 1 - \sqrt{1-\frac{d_n^2}{4}} \\\end{align}$$

Using Pythagoras's theorem on the right angled triangle CBD

$$\begin{align}CD^2 &= BC^2 + BD^2 \\ &= \left(\frac{d_n}{2}\right)^2 + \left(1 - \sqrt{1-\left(\frac{d_n}{2}\right)^2}\right)^2 \\ &= \frac{d_n^2}{4} + \left(1 - 2\sqrt{1-\frac{d_n^2}{4}} + \left(1 - \frac{d_n^2}{4}\right)\right) \\ &= 2 - 2\sqrt{1-\frac{d_n^2}{4}} \\CD &= d_{2n} \\ &= \sqrt{2 - 2\sqrt{1-\frac{d_n^2}{4}}} \\\end{align}$$

d_2n is the length of one side of the polygon with 2N sides.

This means that if we know the length of the sides of a polygon with Nsides, then we can calculate the length of the sides of a polygon with2N sides.

What does this mean? Lets start with a square. Inscribing a square in acircle looks like this, with the side length √2. This gives an estimatefor π as 2√2, which is poor (at 2.828) but it is only the start of theprocess.

We can the calculate the side length of an octagon, from the side lengthof a square, and the side length of a 16-gon from an octagon etc.

It's almost time to break out python, but before we do so, we'll justnote that to calculate d_2n from d_n, the first thing you do iscalculate d² and at the very end you take the square root. This kind ofsuggests that it would be better to keep track of d² rather than d,which we do in the program below(pi_archimedes.py):

import mathdef pi_archimedes(n): """ Calculate n iterations of Archimedes PI recurrence relation """ polygon_edge_length_squared = 2.0 polygon_sides = 4 for i in range(n): polygon_edge_length_squared = 2 - 2 * math.sqrt(1 - polygon_edge_length_squared / 4) polygon_sides *= 2 return polygon_sides * math.sqrt(polygon_edge_length_squared) / 2def main(): """ Try the series """ for n in range(16): result = pi_archimedes(n) error = result - math.pi print("%8d iterations %.10f error %.10f" % (n, result, error))if __name__ == "__main__": main()

If you run this, then it produces:

Hooray! We calculated π to 8 decimal places in only 13 iterations.Iteration 0 for the square shows up the expected 2.828 estimate for π.You can see after iteration 13 that the estimate of π starts gettingworse. That is because we only calculated all our calculations to thelimit of precision of python's floating pointnumbers (about 17digits), and all those errors start adding up.

We can easily calculate more digits of π using Pythons excellentdecimal module. Thisallows you to do arbitrary precision arithmetic on numbers. It isn'tparticularly quick, but it is built in and easy to use.

Let's calculate π to 100 decimal places now. That sounds like asignificant milestone!(pi_archimedes_decimal.py):

from decimal import Decimal, getcontextdef pi_archimedes(n): """ Calculate n iterations of Archimedes PI recurrence relation """ polygon_edge_length_squared = Decimal(2) polygon_sides = 2 for i in range(n): polygon_edge_length_squared = 2 - 2 * (1 - polygon_edge_length_squared / 4).sqrt() polygon_sides *= 2 return polygon_sides * polygon_edge_length_squared.sqrt()def main(): """ Try the series """ places = 100 old_result = None for n in range(10*places): # Do calculations with double precision getcontext().prec = 2*places result = pi_archimedes(n) # Print the result with single precision getcontext().prec = places result = +result # do the rounding on result print("%3d: %s" % (n, result)) if result == old_result: break old_result = result

You'll see if you look at the pi_archimedes function that not a lothas changed. Instead of using the math.sqrt function we use the sqrtmethod of the Decimal object. The edge gets initialised toDecimal(2) rather than 2.0 but otherwise the methods are the same.The main method has changed a bit. You'll see that we set theprecision of the decimal calculations using thegetcontext().prec = ... call. This sets the precision for allfollowing calculations. There are other ways to do this which you cansee in the decimal moduledocs. We do the Archimedes calculations with 200 decimal placesprecision, then print the result out with 100 decimal places precisionby changing the precision and using the result = +result trick. Whenthe result stops changing we end, because that must be π!

If you run this, you'll find at iteration 168 it produces 100accurately rounded decimal places of π. So far so good!

There are two downsides to this function though. One is that the decimalarithmetic is quite slow. On my computer it takes about 6 seconds tocalculate 100 digits of π which sounds fast, but if you were to try for1,000,000 places you would be waiting a very very long time! The otherproblem is the square root. Square roots are expensive operations, theytake lots of multiplications and divisions and we need to do away withthat to go faster.

In Part 3 we'll be getting back to theinfinite series for arctan and on to much larger numbers of digits ofπ.