This post looks at computing by hand a result for pi using a process first described in the 3rd century. This same algorithm yielded the world’s most accurate valuation for the ratio of circle circumference to diameter until the 16th century.

A Short History

Many early civilisations, approximated the ratio of circle circumference to diameter to be 3. Mathematicians in China before the 3rd century had produced various other esimates to close in on a more accurate valuation.

Liu Hui (3rd century), unsatisfied with these approximations, describes a method to produce successively more accurate valuations for the ratio of a circle’s circumference to its diameter. His process looks at the perimeter of an n-sided regular polygon inscribed in a circle and provides a method of doubling the number of sides to produce better approximations with minimum computation[1].

There is some uncertainty as to which approximations can be attributed to Liu, although it is known that he obtained pi to be 3.14 (2 decimal places) by applying the procedure up until a 96-sided shape.

Zu Chongzhi (429-500 AD) later produced an approximation yielding pi between 3.1415926 and 3.1415927 – correct to 6 decimal places. This level of precision is unmatched until the 16th century. It is believed that this was achieved using Liu’s algorithm[1].

The method can be compared to that used by Archimedes who used inscribed and circumscribed polygons with different numbers of sides to produce bounds for the area of the circle and hence pi.

Calculation Strategy

The mathematics and a modern derivation of Liu Hui’s algorithm is explained in a separate post. Using the slightly reformulated expression, we’ll aim to recreate Liu’s result (approximation derived from the perimeter of an inscribed 96-sided shape) using a similar means of doubling the number of shape sides. Our circle will be a unit-circle (r = 1); and each iteration requiring only an addition and a square root:

\[k_{6} = 1\] \[k_{2n} = \sqrt{2 + k_n}\] \[A_{n} \approx \pi \approx \frac{n}{2} \times \sqrt{2 - k_{n}}\]

This gives the final algorithm \(A_n\):

PROCEDURE approximate(iterations) IS
    n := 6
    kn := 1
    FOR I := 1 TO iterations LOOP
        kn := sqrt(2 + kn)
        n := n * 2

    sn := sqrt(2 - kn) 
    An := sn * n / 2

    RETURN An

Intermediate results will be given as decimal values with 17-digits of precision (as many as this human calculator can fit on a sheet of paper: much more than required, but useful if we want to continue with shapes with many more sides in the future – more on this later).

Note that while it’s a similar logic to what Liu Hui followed, the execution is quite different. For a start, the calculation (including square root extraction) would have been calculated using rod calculus – using physical counting rods to represent digits in calculations. Maybe one day this could be interesting to learn more about and try. A result would have been expressed as a pair of integers[2] corresponding to the circumference and diameter of the circle rather than a final decimal ratio.

For now, I’ll stick with paper and a ‘modern’ 16th century method for computing square roots on paper. I describe this method in more detail in a separate post.

Working

Iterations

n	Working
6	The side-length of a hexagon in a unit circle is 1 so \(s_n = 1\). \(k_n = (2 - {s_n}^2)\) \(\Rightarrow k_6 = 1\)
\(k_{12} = \sqrt{2+k_6}\)
\(k_{24} = \sqrt{2+k_{12}}\)
\(k_{48} = \sqrt{2+k_{24}}\)
\(k_{96} = \sqrt{2+k_{48}}\)
\(k_{192} = \sqrt{2+k_{96}}\)

Approximation (A₉₆)

\(A_{96} \approx \pi \approx \sqrt{2 - k_{96}}\)

Approximation (A₁₉₂)

\(A_{192} \approx \pi \approx \sqrt{2 - k_{192}}\)

Final Result

Calculating the difference between the results for the approximations \(A_{192}\) and \(A_{96}\) gives \(D_{192}\). Using the bounds from Liu Hui’s inequality \(A_n < \pi < A_n + D_n\) we can verify the precision of the result.

Hence \(\pi = 3.141\), with a precision of 4 decimal digits (3 decimal places).

Limits and Intermediate Precision

After all that manual computation, the result is admittedly a little underwhelming. The intermediate calculations are calculated to a precision (17 digits) far beyond that required for this result.

But this raises the question – how far could we take this? And how does the precision of the intermediate results affect the final number of correct digits we could achieve with this algorithm?

A small program using arbitrary precision decimals is written to compare the final result with the known value of pi. With intermediate calculations limited to 17 digits of precision as done here, we can get a best final result with 9 correct digits of pi (3.141592653) after 13 iterations of the algorithm (+ a few more hours of computing square roots by hand).