Tag: pi

Posted on in Mathematics, Programming

Liu Hui’s algorithm for calculating Pi

Liu Hui [刘徽] was a Chinese mathematician who lived in the Wei period. One of the things he was famous for was his method of calculating \pi iteratively.

This approach to calculate \pi (also used by Archimedes), is to calculate the area of a polygon inscribed in a circle. As the number of sides of the polygon increases, its area becomes closer and closer to the area of the circle.

Liu Hui found a simple formula to find the side length of an inscribed polygon of 2n sides, if the side length of a polygon with n sides is known:

\begin{array}{l} k_6 = 1 \\ k_{2n} = \sqrt{2+k_n} \\ S_n = \sqrt{2-k_n} \end{array}

Here k is a temporary variable, and S_n is the side length of an inscribed polygon with n sides. Something like this:

We start with a hexagon. The radius of the circle is 1, the area \pi. The side length of the hexagon is 1.

To calculate the next k value, all we need to do is do an addition and a square root:

\begin{array}{ll} k_6 = 1 & S_6 = \sqrt{2-1} = 1 \\ k_{12} = \sqrt{2+1} & S_{12} = \sqrt{2-\sqrt{2+1}} \approx 0.518 \\ k_{24} = \sqrt{2+\sqrt{2+1}} & S_{24} = \sqrt{2-\sqrt{2+\sqrt{2+1}}} \approx 0.261 \end{array}

The area of a regular polygon is A = \frac{1}{2}nsa where n is the number of sides, s is the side length, and a is the apothem. As the number of sides increases, the apothem becomes closer and closer to the radius, so here we’ll just let a=1.

We then have the formula for the area of the polygon: P_n = \frac{1}{2} n S_n, where P_n is the area of a polygon with n sides.

Let’s do a few here:

\begin{array}{l} P_6 = 3 \sqrt{2-1} = 3 \\ P_{12} = 6 \sqrt{2-\sqrt{2+1}} \approx 3.106 \\ P_{24} = 12 \sqrt{2-\sqrt{2+\sqrt{2+1}}} \approx 3.133 \\ P_{48} = 24 \sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2+1}}}} \approx 3.139 \end{array}

To save us some tedious work, here’s a Haskell program to do the calculations for us:

k 6 = 1
k n = sqrt $ 2 + k (n/2)
s n = sqrt $ 2 - k n
p n = (n/2) * s n
main = mapM_ (putStrLn . (\n -> show (round n) ++ " " ++ show (p n))) $ iterate (*2) 6

Here’s part of the output:

6 3.0
12 3.1058285412302498
24 3.132628613281237
48 3.139350203046872
96 3.14103195089053
192 3.1414524722853443
384 3.141557607911622
768 3.141583892148936
1536 3.1415904632367617
3072 3.1415921060430483
6144 3.1415925165881546
12288 3.1415926186407894
24576 3.1415926453212157
49152 3.1415926453212157
98304 3.1415926453212157

This is about the maximum you could get without arbitrary precision decimals: the area of a polygon with 98304 sides. This gives us 3.14159265 which is 8 digits of \pi.

With a big-decimal library, it would be possible to calculate \pi to any precision.

Happy pi day!