# 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!

## 6 thoughts on “Liu Hui’s algorithm for calculating Pi”

How about the error analysis? How do we know how many digits of our result are correct, given some precision on taking the square roots?

1. Hmm, I actually don’t know how you would calculate that. But I imagine it converges very slowly (compared to modern algorithms), and you need to store a lot of digits for intermediate values if you want to calculate pi digits.

2. Pingback: Who Discovered Pi?

1. bob says:

Liu Hiu lived in the common era.. Not before Christ

3. Anonymous says:

it would be cool if you would write this article in a way a layman would undestand

4. I have a method of calculating pi which I think is new. For N operations, it produces O(N^2) accurate digits. The methods I’ve read about in the literature seem to be:
Euler’s series, so no. of digits D is about log10(N)
power series, D is proportional to N
quadratically-convergent methods, so D is about 2^N or 3^N or (other small integer)^N

My method is in between: D is proportional to N^2.

It’s a two-stage method. It uses p steps of a first algorithm, then q steps of a second algorithm. The error is about (1/2)^(pq). The number of operations is about p+q times a constant. One step of the first algorithm uses two additions, a multiply and a square root. One step of the second algorithm uses a multiply, an add and a divide. The division is by an integer.

I have a (conceptually) related algorithm for sine and cosine. Again, it’s a two-part method with similar properties of computation time and size of error.
Anyone interested?