**This haunts me …**

There is beauty here. It’s art, no less valid as art than a poem, or a sculpture, or a fine painting. It is simple and clean, yet infinitely precise. It is so simple that one may memorize it in seconds, yet it defines an infinite number of digits of an irrational mathematical constant. *What poem can match that?*

As simple as that formula is, one should be able to sit down with nothing more than pen and paper and derive it. It *has* to be possible: Draw some circles and triangles, think, rationalize, doodle, and come out with that formula. Thus far, I’ve not been able to do that. But I refuse to look up a solution, so I return to it from time to time.

I’ll get back to that beautifully simple formula in a bit, but first I would like to highlight the sloppy alternatives. There are any number of approaches that one may consider for approximating pi. The most straight forward, in my view, is the method if dissecting a circle into increasingly slender isosceles triangles and allowing the bases of these triangles to form the sides of a polygon, approximating a circle. Start with a hexagon embedded in a circle as shown in the picture to the right and divide the hexagon into six identical triangles. We’ll let the radius of the circle equal one, and since in this case the six triangles are equilateral, all sides are of of this length. This gives a first approximation of π of 6/2, or 3.

Next, divide each triangle in half. This results in a 12-sided polygon comprised of isosceles triangles of the type ABC in the diagram. The two sides AC and BC are still equal to 1 in length. Length AD is one half of the original triangles base, so it is 0.5. However, the length we need in order to calculate π is AB, and that will be slightly longer than 0.5 since it is the hypotenuse of the right triangle ADB. To find AB, we need BD, which we don’t have. However, we know that BC is equal to 1, and we can calculate CD since it is a side of the right triangle ADC, with hypotenuse AC equal to 1 and side AD equal to 0.5. Using the Pythagorean theorem, we find that CD is 0.866, meaning that BD must be equal to 0.134. So then, to calculate hypotenuse AB, we use the Pythagorean theorem again for right triangle ADB, with sides of 0.5 and 0.134. This gives us a length of AB as 0.5176, leading to our new approximation: π = (12 X 0.5176)/2 = 3.1056.

So yeah — by now, you aren’t reading any of this headache-inducing shit. But if you dick around with this stuff on paper a little while, you’ll come up with a formula for calculating the new side length whenever you split the triangles and double the number of sides of the polygon. After a couple of Advil, your formula will simplify down to look something like this:

And with each doubling of the number of sides of the polygon and recalculation of the side length, you’ll find that you approach π asymptotically as follows:

Polygon Sides | Side Length | Pi |

6 | 1.00000 | 3.000000 |

12 | 0.51764 | 3.105829 |

24 | 0.26105 | 3.132629 |

48 | 0.13081 | 3.139350 |

96 | 0.06544 | 3.141032 |

192 | 0.03272 | 3.141452 |

384 | 0.01636 | 3.141558 |

768 | 0.00818 | 3.141584 |

1536 | 0.00409 | 3.141590 |

3072 | 0.00205 | 3.141592 |

That was the painful way of approaching π, and it is clearly a lot messier and migraine catalyzing than the formula I showed at the beginning of this post. Let’s get back to that …

This is altogether a different kind of animal. Rather than approaching π asymptotically, this formula is an oscillator and decays towards π the further out one carries the calculation, as shown below.

It’s a painfully slow oscillation. Even after one hundred calculations, it still oscillates between about 3.13 and 3.15. However, because it oscillates evenly above and below the value of π, one may drastically boost the precision by applying a simple 2-point moving average as shown in the graph above. Furthermore, one may repeat this process, applying a 2-point moving average to the previous average and so on.

After three such applications of a 2-point moving average, the error is now only 0.000000008. Not too shabby! Look again at the simplicity of the formula from which this precise calculation of π was generated. *Tell me this isn’t art!* 😛

Clearly, this entire post has been one horrific exercise in overwhelming geekdom. But that’s okay! Mental exercises such as this keep you sharp and help you see the world in different ways. Some people like crossword or jigsaw puzzles; some like to calculate π in new ways and from different geometric constructs. Whether you are solving a word puzzle or calculating π, the process is one of creation, and the result is discovery and satisfaction. (Unless of course you fuck it up. Then … not so satisfying.) 🙂

**The Monte Carlo Approach …**

In the comments section, KAB posted a link to the Monte Carlo approach to π and I started playing around with that. I’ll give the method an “A” for simplicity but an “F” for precision and ease of execution. It is a beautiful approach and the math is simple. But it’s the sort of thing you would want to program as a loop and leave running overnight (or over a few days). If all you have to work with are simple spreadsheets, as is the case with me right now, this approach isn’t so hot.

The premise is that, if you inscribe a circle in a square, the ratio of the areas will be equal to π/4. You can determine the ratio of the areas by randomly picking points within the square and determining whether or not they are also inside the circle. The easiest way to go about this is to use Cartesian coordinates.

Use a value of “1” for the radius of the circle and use a random number generator to pick two random numbers — one for X and one for Y. Then, sum the squares of these numbers to determine the square of the radial distance this point is from the origin. If the result is a number greater than 1, then the point is outside of the circle. If, on the other hand, the point is equal to or less than 1, then the point is inside both the circle and the square. Note that I’m only calculating an r-squared here: You don’t need to bother taking the square root of the number since it won’t effect whether the result is greater than or less than one. Let’s say you do this 1,000 times, choosing random X,Y pairs and summing the squares. Out of those 1,000 points, you may get 800 values that are less than or equal to one and are therefore inside both the circle and the square. So, to calculate π, you multiply 800/1000 by 4 to get 3.2.

The method is not too terribly efficient, obviously. But it does have one advantage: It does not rely on previous calculations. You can calculate 1,000 X,Y pairs today, jot down the number of results that are less than or equal to one, and do the same thing the next day with fresh calculations. Furthermore, while the methods I described previously in this post eventually result in calculation errors due to computer limits, this method avoids those problems. Finally, while the example result I gave above of 3.2 is pretty far from the actual value of π, that was from only 1,000 calculations. A trillion calculations would give a much better answer! 🙂

This method neither approaches π asymptotically nor as a decaying oscillation. It just kind of ** meanders** towards π. If ever an “LOL” is appropriate in a dry mathematical discussion, this would be it. LOL! 😛

It is worth noting that the Monte Carlo approach is going to be very susceptible to any biases in rounding or in the methods of random number generation. I used Google Sheets for my calculations and to generate the graph shown above. In my first attempt, using 160,000 XY pairs, I found that my approximation seemed to drift a little to the high side and my final result was 3.146, which isn’t really that great. I wasn’t sure if I had a real bias problem or if it was just a coincidental artifact that would be smoothed out with more data, so I decided to redo the project. I repeated the approximation using one million XY pairs and obtained a much better result of 3.1415. I’m confident now that the “drift” problem I had previously was just a localized artifact.

It may sound like one million XY pairs is a lot of data, but actually it isn’t. Each of those calculations, regardless of the number of digits of the generated random numbers, is ultimately converted to a simple “yes” or “no,” depending on whether the sum of the squares is greater than or less than 1. My calculation, therefore, breaks down to simple count of “no” values: 4 x 785,375 / 1,000,000 = 3.1415. The closest I could have come to π would have been 785,398 “no” values, or 3.141592.

My favorite method is using Monte Carlo simulations to approximate pi:

http://mathfaculty.fullerton.edu/mathews/n2003/MonteCarloPiMod.html

Essentially throwing darts at a circle drawn in a square

Never tried that, but I might have to dick around with it at some point.

The first method I posted — calculating from polygons if increasing number of sides, is actually Hell on a spreadsheet. Google sheets went totally wonky on me after about 25,000 sides, since I was a number too small for Google sheets to handle properly. It kind of goes crazy for a little while, giving various goofy results, and then it fully goes to shit.

That’s one of the nice features of the oscillating equation …. it will handle much, much smaller numbers before crapping out.

Post updated to include Monte Carlo. 🙂

“It is worth noting that the Monte Carlo approach is going to be very susceptible to any biases in rounding or in the methods of random number generation.”

Yes, I agree, and it is a good test for random number generators. Of course you have to run it to infinity to get Pi if it is a perfect generator.

But if you don’t have an analytic solution to a problem, and you can simulate using Monte Carlo methods, it is fantastically approximating a useful result in the real world.