Calculating Pi

The ratio of the circumference of a circle to its diameter has been known as a constant for several thousand years, although it was only named π by Leonard Euler in 1737. The Ancient Babylonians thought that π was 3 1/8, and they probably measured it experimentally by use of a thin rope. You can do this by driving a stake into the ground, stretching a rope from that stake to inscribe a circle as a groove, and then marking off the radius multiple times by placing the rope in the groove. It’ll go round six and a quarter times, roughly. (Note that it’s easy to measure off powers of a half by folding the rope at the midpoint, so it’s not surprising that they settled on 3 1/8.)

This value, and sometimes 3 1/7, was the accepted value of π through much of antiquity. Another value was 3: in the Old Testament, I Kings 7, verse 23 states: “And he made a molten sea, 10 cubits from the one brim to the other: it was round all about, and his height was five cubits: and a line of 30 cubits did compass it round about.” This indicates that the editors of the Bible thought π was 3, despite the experimental findings of the Ancient Egyptians.

Fast forward a few hundred years to about 250 BCE, and we meet up with Archimedes, one of the early giants of mathematics and physics. He applied a very modern view of calculating the circumference of a circle: rather than worrying about trying to find an exact value, he worked out a way of calculating an upper bound and a lower bound to the value of the circumference. That method could then be refined to give better and better estimations of the true value by reducing the upper bound and increasing the lower bound.

The method he used was to inscribe a regular polygon inside the circle. The perimeter was smaller than the circumference, forming the lower bound. Circumscribing a regular polygon outside the circle would provide the upper bound, since its perimeter would be larger than the circumference of the circle. Figure 1, below, shows the situation that Archimedes considered for a hexagon.

We can do this in our heads using basic trigonometry. The inner hexagon’s sides are all equal to the radius, and so the perimeter of the inner hexagon is six times the radius. The outer hexagon is slightly more complicated, but in essence the radius of the circle is the length of the perpendicular to the middle of the side. Using Pythagoras’ Theorem, we can work out that the length of the side is 2 ≈ √(1/3) times the radius, and therefore that the perimeter is 6.93 times the radius. From this we can deduce that π is between 3 and 3.47. Not very accurate, perhaps, but Archimedes went further than just a hexagon. He used 96-sided polygons to work out that π is between 3 10/71 and 3 1/7.

The most impressive thing about all of this is that he managed to do these calculations without the use of trigonometry (it had not yet been invented), and, even more impressively, without the use of our modern number notation (that is, the 10 digits, including 0, and positional notation). Things get even more astonishing when you remember that Archimedes lived under Roman rule, and they used the most awkward number notation possible for doing arithmetic: II + III = V, indeed. Now try and calculate the square root of 3 using this notation (as it happened, Archimedes used the ratio 265/153 to do this).

Despite all of these difficulties, Archimedes’ results remained the basis of calculating π for pretty much two thousand years. Fibonacci (he of the sequence 1, 1, 2, 3, 5, 8, 13... where each term is the sum of the previous two) had a crack at improving it in 1220 using the decimal notation newly imported from the Middle East, but couldn’t do much better than Archimedes: his approximation was 864/275, or 3.141818.

François Viète was possibly the last mathematician to use the Archimedean method for calculating π, and a pretty heroic attempt it must have been too. In 1593, he used a polygon of 393,216 sides to calculate that π was between 3.1415926535 and 3.1415926537. I can’t imagine the number of days it must have taken to work this out without being able to use any calculating devices.

After that time, new analytic methods were starting to appear that finally drove the nail into the coffin of Archimedes’ polygon method for calculating the value of π. In 1706, John Machin, then a professor of astronomy at Gresham College in London, devised a formula that would be used for hundreds of years afterwards to calculate ever more digits of π (he himself used it to calculate 100 decimals). The formula is trigonometric:

π/4 = 4 ≈ arctan(1/5) – arctan(1/239)

Although very concise, this formula doesn’t initially look too promising: how the heck do you calculate the arctan of 1/5, let alone 1/239? If you have a calculator, it’s pretty easy. Try it yourself and use your scientific calculator to evaluate this expression (making sure to set your calculator into radians mode, as these analytical methods use radians for angle measurements, not degrees). If everything works correctly, it will give an answer of 3.14159265359, which, sure enough, is π to 11 decimal places. But it can’t give it to any better approximation than that. If we’d like to calculate, say, 1,000 digits of the expansion of π, my trusty Hewlett-Packard calculator can’t get there from here.

So, apart from looking it up in a book of tables or using a calculator, how can you calculate the arctan of anything? Enter a contemporary of Machin called Brook Taylor. He invented the Taylor series, which was an important result in analytical methods. A Taylor series is an infinite sum of terms, and Taylor had shown the very nice result that arctan(x) = x – x³/3 + x⁵/5 – x⁷/7 + ... where each term in the series is the argument raised to the next odd number, divided by that odd number, and the term is either added or subtracted, alternately.

Furthermore, if the argument is less than 1, the terms decrease in value pretty quickly, and if the argument is close to 0 then they decrease very rapidly. Using the Taylor series, Machin was able to calculate the value of π to 100 decimal places, and others went even further in their quest for the number. Nowadays, π is known to some trillion digits (1,241,100,000,000, according to Wikipedia). For reasons of ease of calculation with binary computers, it’s usually calculated in hexadecimal and then converted to decimal.

But the question remains for us: how do we, with our mortal desktop PCs, do this calculation? After all, the floating point types we use have only 17 significant digits, which is just not enough to get 20 digits of π, let alone 1,000 or more. There’s nothing for it but to write a set of routines to deal with ‘big numbers’, numbers with thousands of digits. Big numbers are generally implemented as arrays of integers, and for ease of use, as fixed point numbers. That is, one part of the big number array is the part of the number before the decimal point and the other part the digits after it. In our example, the part before the decimal point is going to be small (after all, π is 3 point something), and the part after the point is going to be where all the action occurs.

To make the arithmetic conceptually easier to understand, I’ll make each element of the array an integer that contains a number from 0 to 99,999 – that is, be responsible for five decimal digits of the final number. This will also help in the final printing out of our 1,000 digits of π: after all, we’ll need to be able to read the result of our calculation. One way to look at this concept is to say that we’re going to be doing arithmetic in base 100,000, rather than base 10.

Let’s rearrange Machin’s formula a bit:

π = 16 ≈ arctan(1/5) – 4 ≈ arctan(1/239)

Now we can see that our big number routines must include multiplication by an integer (by 16, and by 4), subtraction of two big numbers and an arctan routine that acts on the reciprocal of an integer (so would seem to indicate that we should have division by an integer as well).

Looking at the Taylor series for arctan, it would confirm that we need division by an integer, and also addition of two big numbers. The series can be then calculated by having a big number sum and two temporary big numbers, the first to store the powers of the reciprocal, and the second as the store to calculate each term.

So, all in all, we have five calculation routines and a single routine to set up a big number variable in the first place. For 1,000 digits, we need 200 elements of the array after the fixed decimal point, and one before. It makes sense to have some extra elements on the end as well to control underflow, so we’ll make the big number array 210 elements long (by the way, I’m assuming that integers are 32-bit entities). The routine that sets up the big number must allocate the 210 element array and initialize all elements to zero.

Let’s tackle the easiest calculation routine first: addition. What we have are two arrays of integers, both the same size, with each element representing five decimal digits, and we’ll assume that the second operand gets added to the first (this means we don’t have to have a separate big number result). We proceed just as we learned in school when we first came across addition. We start from the right, and add up corresponding elements. If the result of adding two elements is greater than or equal to 100,000, we carry 1 to the next element addition to the left. Figure 2, below, shows this fairly simple concept in action.

By now you’ll have worked out that subtraction also works the same way as you learned at school. We go from right to left again, and we use the concept of a ‘borrow’ of 1 from the element to the left if the value being subtracted is less than the value being subtracted from. See Figure 2 again to see a demonstration of how this works. For both of these operations, we shall ignore overflow. That is, if the resulting big number on addition is greater than or equal to 100,000.0, we’ll lose the thousands part; if the result of a subtraction is less than 0.0, we’ll be wrapping. For this implementation, neither of these will happen.

Multiplication by an integer is the next step. It will come as no surprise that we use our trusty primary school knowledge again here. We proceed from right to left, multiplying each element by the integer. We have to take account of carry again, as in the addition example. But note something interesting: since we only multiply by a maximum of 16, the result of multiplying an element will not overflow an integer (in fact, we can multiply by a maximum of about 20,000 with safety). Figure 3, below, shows the operation.

Next on the list of basic operations is division by an integer. Here, we’re going to be treading on some very thin ice indeed. First off, we proceed from left to right. We divide each element (that is, big number ‘digit’) by the divisor to give us the resulting digit at that position. We then ‘mod’ the element value by the divisor to give us the carry for the next element. This carry, as we’ve indicated before, cannot exceed 20,000, because for the next element we have to take the current carry and multiply it by 100,000 and add it to the element value before doing the next division. This is, as they say, definitely thin ice.

For a start, we cannot divide by 239 ≈ 239 (as we’d expect to do) since that’s 57,121 and it would produce remainders greater than 20,000. Hence we’ll have to divide by 239 twice in our loop. How about the odd divisor? How big will that grow? My calculations show that, for 10,000 digits of π, it won’t get above 14,310, since 0.2^14,310 has more than 10,000 zeros after the decimal point (and would therefore be assumed to be zero).

Finally, having the four basic operations, we can easily combine them to calculate π from Machin’s formula and the Taylor series. In 1957, the year I was born, a Ferranti Pegasus 1 computer calculated 7,480 digits of π, a record at the time. It took 33 hours. Using this algorithm implemented in C# and a modern quad-core desktop to run it on, I just now worked out that the 7480th digit after the decimal point was 4. It took less than a second. That’s some pretty impressive development.

Julian Bucknall