Pi and Pythagoras

1.0 Introduction

This post contains a derivation of a recurrence relation that can be used to calculate π. While it's not useful for the purpose of calculating to any decimal place imaginable, what's notable is that it only uses the Pythagorean thereom and should therefore be easy to follow. Unsurprisingly, the recurrence relation is neither original nor unique because it falls under the category of using polygons to compute π (attributed to Archimedes). At the end of the notebook, the derived recurrence relation is implemented in code and run.

There is more than one way to use polygons in the computation of π. One could follow Archimedes' proof, which uses a theorem from Euclid's Elements. There are other discussions that make use of trigonometric functions. Of course, one could just use the Pythagorean theorem by reducing the n-sided polygons to right triangles, as done here.

I went through a few pages of google search results and found similar posts. (I posted this anyway because it was fun and mainly an excuse to create an IPython notebook.) The first link's code can achieve a much larger decimal place than that shown here because the author used a specialized decimal type instead of the default float type available in Python.

A good historical account of π that I found while briefly searching is found here.

2.0 Problem set up

Consider one quarter of a circle with arbitrary radius, r, between 0 and 90 degrees. The below will mathematically describe the outer perimeter of many triangles drawn inside this region. The outer perimeter is composed of the triangles whose hypotneuses are closest to the circle's arc. The quantity of triangles drawn will be proportional to an integer index, n. All variables make reference to n via a subscript.

Variable definitions

  • r is the radius.

  • Tn refers to any right triangle newly established at the n-th iteration.

  • Ln refers to any right triangle's leg.

  • Hn refers to the hypotneuse of any small triangle (TnS) that composes a portion of the outer perimeter.

  • Pn is the outer perimeter.

Subscript meanings

  • Subscripts B and S refer to any big and small triangles, respectively.

  • LnB indicates a leg that is unique to TnB.

  • LnS indicates a leg that is unique to TnS.

  • LnC indicates a leg that is common to TnB and TnS.

Below are explanatory images that indicate most of the above variables. They are for n=1, and the procedure to arrive at the drawn line segments is outlined in the following subsections. The first figure indicates triangles, and the second figure indicates line segment lengths. Omitted from the second figure are the radius labels (the radii are at 0o, 45o, and 90o). The figures are an attempt to convey that there are two sets (big/small) of identical triangles at each iteration, and that by calculating the properties of one triangle in each set, all other triangles' properties are known.

Blue lines always correspond to refinements made at the previous iteration. Green lines are the new refinements made at the current iteration. This combination of line segments produces the set of TnB and TnS at iteration n.

3.0 Analysis

3.1 Analysis for n = 0

Inside the starting region a right triangle is created by connecting points (0,r) and (r,0). The hypotneuse can be easily calculated by the Pythagorean thereom, and for n=0 the hypotneuse happens to be equal to the outer perimeter.

H0=r2+r2P0=H0=2r

3.2 Analysis for n = 1

The previous right triangle's hypotneuse (H0) is bisected by a line from the origin to the circle's perimeter at an angle exactly halfway between the legs (or catheti) of T0. This action creates four new right triangles since the hypotneuse of T0 forms a right angle with the bisecting line. These new triangles are two of T1B and two of T1S.

Since the hypotneuse of T0 has been bisected, the leg that is common to any T1B, T1S pair is known (this is L1C). The hypotneuse of T1B is known because it is just T0's leg, r. Therefore, we can solve for the missing leg, L1B.

L1C=H02L1B=r2L21C

The line segment labeled L1S is just the radius minus L1B. This permits calculation of H1, which is the hypotneuse of the small triangle T1S.

L1S=rL1BH1=L21S+L21C

Finally, the outer perimeter defined by the two new small right triangles is twice the value of a single small triangle's hypotneuse, H1.

P1=2H1=2L21S+L21C

3.3 Analysis at n = 2 and above

As the index increases, the efforts in the above are repeated. 2n+1 right triangles are created, of which half are TnB and half are TnS. Any TnB will have a leg that's shared with a TnS called LnC. This common leg is exactly half the hypotneuse from the previous iteration. Using this information one can solve for LnB, LnS, and Hn.

LnC=Hn12LnB=r2L2nCLnS=rLnB

With the legs of TnS known, one can find the hypotneuse of TnS using the Pythagorean theorem.

Hn=L2nS+L2nC=(rr2(Hn12)2)2+(Hn12)2

The above is a recurrence relation for Hn for all n>0. The initial value derived above, H0, is just 2r. The perimeter at an arbitrary index is comprised of 2n right triangles' hypotneuses.

Pn=2nHn=2n(rr2H2n14)2+H2n14

Below are plots for n=2 and n=3. Only Tn and Tn1 are drawn to reduce graphical clutter. The new iteration's triangles become more difficult to see as the index increases so plots of n4 are not shown.

4.0 Computing π

The circumference of the circle is approximated by 4Pn, where the factor of four is present since we only considered the region between 0 and 90 degrees (all other regions are identical). The quantity 4Pn will always be slightly smaller than the true circumference but would converge in the limit as n. Equating 4Pn with the true circumference, π2r, and solving for π yields:

π=limn2Pnr

Since computing something to infinity is not possible, let's accept that we're stuck with an approximation and recognize it using in the below expression for π. These equations are readily implemented in code with a for loop and can be evaluated for as large an n that our computing resources will permit.

H0=2rHn=(rr2H2n14)2+H2n14 , for n>0π2n+1rHn

The downside to this approach is that as n increases, Hn gets smaller and smaller. For example, consider the unit circle with r=1. As n increases, the above equation presents a precision challenge for any computer since it involves the difference of a constant (unity) with an ever-decreasing quantity (H2n1/4) that approaches zero as n.

4.1 In terms of incremental area

It is also possible to use the triangles' areas to compute π. The additional area contributed by each iteration can be computed by analyzing TnS, the area of which is written below by using the Pythagorean theorem for LnS.

ATnS=12LnCLnS=12Hn12H2n(Hn12)2=Hn14(rr2(Hn12)2)

There are 2n such TnS in the region between 0 and 90 degrees. Therefore, multiplying ATnS by 2n+2 gives the total incremental area at an arbitrary index, ΔAn.

ΔAn=2n+2ATnS=2nHn1(rr2(Hn12)2)

The total area at some index, An, is an accumulation of ΔAn plus the initial value, A0.

A0=2r2An=A0+nΔAn , for n>0

As for the perimeter case, if the area-based result is used to compute π then it is at best an approximation. When implemented in code it has the same precision limitations as the perimeter case.

πAnr2

5.0 Python implementation

In [8]:
from math import sqrt

# http://www.piday.org/million/
PI = '3.141592653589793238462643383279502884197169399375' # and many more ...

def compare_pi(value):
    """
    Compare the approximation to pi given by 'value' to the actual value.
    - 'value' argument must be a float
    """
    # first, convert 'value' into a super-long string
    value = '%.050g' % value
    # next, iterate over the characters in the PI string
    # if there is a mismatch between PI and value, stop iterating
    for i in xrange(len(PI)):
        if PI[i] != value[i]:
            print 'correct portion = %s' % value[0:i]
            # the conditional ignores each string's decimal symbol
            print 'mismatch at digit %d' % (i if i > 0 else 1)
            break

# assign the radius and number of iterations
# make r a float to avoid unintentional integer division
r = 1.0
# in its current form, the best this script can do is 16 decimal places
# this occurs for iterations >= 24
iterations = 24

# during every loop:
#   (1) compute the current hypotneuse using the previous hypotneuse
#   (2) compute pi using the current hypotneuse
for n in xrange(iterations):
    if n == 0:
        H_current = sqrt(2) * r
    else:
        H_current = sqrt(
            (r - sqrt(r ** 2 - H_previous ** 2 / 4)) ** 2 +
            H_previous ** 2 / 4
        )

    pi_approximated =  2 ** (n + 1) * H_current / r
    H_previous = H_current

    # conditional that avoids printing each iteration's details
    if n > 1 and n < iterations - 2:
        if n in [2, 3, 4]:
            print '...'
        continue

    print '--------------------------------------------------------'
    print 'iteration number %d' % (n + 1)
    print 'n = {:}, which gives {:,} triangles'.format(n, 2 ** (n + 1))
    print 'calculated pi   = %.030g' % pi_approximated
    compare_pi(pi_approximated)
--------------------------------------------------------
iteration number 1
n = 0, which gives 2 triangles
calculated pi   = 2.82842712474619029094924371748
correct portion = 
mismatch at digit 1
--------------------------------------------------------
iteration number 2
n = 1, which gives 4 triangles
calculated pi   = 3.061467458920718698323071294
correct portion = 3.
mismatch at digit 2
...
...
...
--------------------------------------------------------
iteration number 23
n = 22, which gives 8,388,608 triangles
calculated pi   = 3.14159265358977624060798916616
correct portion = 3.1415926535897
mismatch at digit 15
--------------------------------------------------------
iteration number 24
n = 23, which gives 16,777,216 triangles
calculated pi   = 3.14159265358979000737349451811
correct portion = 3.14159265358979
mismatch at digit 16

Download this post's IPython notebook here.

Comments !

social