Nick Craig-Wood's Articles

2019-12-25

Pi - Archimedes

2011-07-18 00:00

It is clear from Part 1 that in order to calculate π we are going to need a better method than evaluating Gregory's Series.

Here is one which was originally discovered by Archimedes in 200BC. Most of the proofs for series for calculating π require calculus or other advanced math, but since Archimedes didn't know any of those you can prove it without! In fact I thought I had discovered it for the first time when I was 13 or so, but alas I was 2000 years too late!

Archimedes knew that the ratio of the circumference of a circle to its diameter was π. He had the idea that he could draw a regular polygon inscribed within the circle to approximate π, and the more sides he drew on the polygon, the better approximation to π he would get. In modern maths speak, we would say as the number of sides of the polygon tends to infinity, the circumference tends to 2π.

Instead of drawing the polygon though, Archimedes calculated the length using a geometric argument, like this.

Diagram for proof of Archimedes π formula

Draw a circle of radius 1 unit centered at A. Inscribe an N sided polygon within it. Our estimate for π is half the circumference of the polygon (circumference of a circle is 2πr, r = 1, giving 2π). As the sides of the polygon get smaller and smaller the circumference gets closer and closer to 2π.

The diagram shows one segment of the polygon ACE. The side of the polygon CE has length d n. Assuming we know d for an N sided polygon, If we can find an expression for the length CD = d 2n, the edge length of a polygon with 2N sides, in terms of d n only, then we can improve our estimate for π. Let's try to do that.

We bisect the triangle CAE to make CAD and DAE the two new equal segments of the 2N sided polygon.

Using Pythagoras's theorem on the right angled triangle ABC, we can see:

\begin{eqnarray*}
AB^2 + BC^2 &=& 1^2 = 1 \\
         AB &=& \sqrt{1 - BC^2} \\
\end{eqnarray*}

Given that

\begin{equation*}
   BC = \frac{d_n}{2}
\end{equation*}

Substituting gives

\begin{eqnarray*}
   AB &=& \sqrt{1-\left(\frac{d_n}{2}\right)^2} \\
   BD &=& 1 - AB \\
      &=& 1 - \sqrt{1-\frac{d_n^2}{4}} \\
\end{eqnarray*}

Using Pythagoras's theorem on the right angled triangle CBD

\begin{eqnarray*}
CD^2 &=& BC^2 + BD^2 \\
    &=& \left(\frac{d_n}{2}\right)^2 + \left(1 - \sqrt{1-\left(\frac{d_n}{2}\right)^2}\right)^2 \\
    &=& \frac{d_n^2}{4} + \left(1 - 2\sqrt{1-\frac{d_n^2}{4}} + \left(1 - \frac{d_n^2}{4}\right)\right) \\
    &=& 2 - 2\sqrt{1-\frac{d_n^2}{4}} \\
CD = d_{2n} &=& \sqrt{2 - 2\sqrt{1-\frac{d_n^2}{4}}} \\
\end{eqnarray*}

d 2n is the length of one side of the polygon with 2N sides.

This means that if we know the length of the sides of a polygon with N sides, then we can calculate the length of the sides of a polygon with 2N sides.

What does this mean? Lets start with a square. Inscribing a square in a circle looks like this, with the side length √2. This gives an estimate for π as 2√2, which is poor (at 2.828) but it is only the start of the process.

Square inscribed in a circle

We can the calculate the side length of an octagon, from the side length of a square, and the side length of a 16-gon from an octagon etc.

Polygons inscribed in a circle

It's almost time to break out python, but before we do so, we'll just note that to calculate d 2n from d n, the first thing you do is calculate d² and at the very end you take the square root. This kind of suggests that it would be better to keep track of d² rather than d, which we do in the program below (pi_archimedes.py):

import math

def pi_archimedes(n):
    """
    Calculate n iterations of Archimedes PI recurrence relation
    """
    polygon_edge_length_squared = 2.0
    polygon_sides = 4
    for i in range(n):
        polygon_edge_length_squared = 2 - 2 * math.sqrt(1 - polygon_edge_length_squared / 4)
        polygon_sides *= 2
    return polygon_sides * math.sqrt(polygon_edge_length_squared) / 2

def main():
    """
    Try the series
    """
    for n in range(16):
        result = pi_archimedes(n)
        error = result - math.pi
        print("%8d iterations %.10f error %.10f" % (n, result, error))

if __name__ == "__main__":
    main()

If you run this, then it produces:

Iterations Sides Result Error
0 4 2.8284271247 -0.3131655288
1 8 3.0614674589 -0.0801251947
2 16 3.1214451523 -0.0201475013
3 32 3.1365484905 -0.0050441630
4 64 3.1403311570 -0.0012614966
5 128 3.1412772509 -0.0003154027
6 256 3.1415138011 -0.0000788524
7 512 3.1415729404 -0.0000197132
8 1024 3.1415877253 -0.0000049283
9 2048 3.1415914215 -0.0000012321
10 4096 3.1415923456 -0.0000003080
11 8192 3.1415925765 -0.0000000770
12 16384 3.1415926335 -0.0000000201
13 32768 3.1415926548 0.0000000012
14 65536 3.1415926453 -0.0000000083
15 131072 3.1415926074 -0.0000000462

Hooray! We calculated π to 8 decimal places in only 13 iterations. Iteration 0 for the square shows up the expected 2.828 estimate for π. You can see after iteration 13 that the estimate of π starts getting worse. That is because we only calculated all our calculations to the limit of precision of python's floating point numbers (about 17 digits), and all those errors start adding up.

We can easily calculate more digits of π using Pythons excellent decimal module. This allows you to do arbitrary precision arithmetic on numbers. It isn't particularly quick, but it is built in and easy to use.

Let's calculate π to 100 decimal places now. That sounds like a significant milestone! (pi_archimedes_decimal.py):

from decimal import Decimal, getcontext

def pi_archimedes(n):
    """
    Calculate n iterations of Archimedes PI recurrence relation
    """
    polygon_edge_length_squared = Decimal(2)
    polygon_sides = 2
    for i in range(n):
        polygon_edge_length_squared = 2 - 2 * (1 - polygon_edge_length_squared / 4).sqrt()
        polygon_sides *= 2
    return polygon_sides * polygon_edge_length_squared.sqrt()

def main():
    """
    Try the series
    """
    places = 100
    old_result = None
    for n in range(10*places):
        # Do calculations with double precision
        getcontext().prec = 2*places
        result = pi_archimedes(n)
        # Print the result with single precision
        getcontext().prec = places
        result = +result           # do the rounding on result
        print("%3d: %s" % (n, result))
        if result == old_result:
            break
        old_result = result

You'll see if you look at the pi_archimedes function that not a lot has changed. Instead of using the math.sqrt function we use the sqrt method of the Decimal object. The edge gets initialised to Decimal(2) rather than 2.0 but otherwise the methods are the same. The main method has changed a bit. You'll see that we set the precision of the decimal calculations using the getcontext().prec = ... call. This sets the precision for all following calculations. There are other ways to do this which you can see in the decimal module docs. We do the Archimedes calculations with 200 decimal places precision, then print the result out with 100 decimal places precision by changing the precision and using the result = +result trick. When the result stops changing we end, because that must be π!

If you run this, you'll find at iteration 168 it produces 100 accurately rounded decimal places of π. So far so good!

There are two downsides to this function though. One is that the decimal arithmetic is quite slow. On my computer it takes about 6 seconds to calculate 100 digits of π which sounds fast, but if you were to try for 1,000,000 places you would be waiting a very very long time! The other problem is the square root. Square roots are expensive operations, they take lots of multiplications and divisions and we need to do away with that to go faster.

In Part 3 we'll be getting back to the infinite series for arctan and on to much larger numbers of digits of π.

Categories: maths, pymath, python | View Comments
Read and Post Comments

Pi - Gregory's Series

2011-07-17 00:00

Lets calculate π (or Pi if you prefer)! π is an irrational number (amongst other things) which means that it isn't one whole number divided by another whole number. In fact the digits of π are extremely random - if you didn't know they were the digits of π they would be perfectly random.

π has been calculated to billions of decimal places, but in this series of articles, we're going to aim for the more modest target of 100,000,000 decimal places calculated by a python program in a reasonable length of time. We are going to try some different formulae and algorithms following a roughly historical path.

You were probably taught at school that 22/7 is a reasonable approximation to π. Maybe you were even taught that 355/113 was a better one. You were unlikely to have been taught how π is actually calculated. You could draw a circle and measure the circumference, but that is very hard to do accurately.

Luckily maths comes to the rescue. There are 100s of formulae for calculating π, but the simplest one I know of is this:

\begin{equation*}
\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - \frac{1}{11} + \frac{1}{13} - \frac{1}{15} + \frac{1}{17} - \cdots
\end{equation*}

This is known as Gregory's Series (or sometimes the Leibnitz formula) for π. That ... stands for keep going forever. To make the next term in the series, alternate the sign, and add 2 to the denominator of the fraction, so + 1/19 is next then - 1/21 etc.

To see why this is true, we need to look at the Leibnitz formula for arctan (or tan⁻¹). This is:

\begin{equation*}
arctan(x) = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + \frac{x^9}{9} - \cdots ( -1 <= x <= 1 )
\end{equation*}

Note that the result of arctan(x) is in radians. I'm not going to prove that formula here as it involves calculus, but you can see a proof of the arctan formula if you want.

arctan is the inverse of tan, so arctan(tan(x)) = x. tan(x) is the familiar (hopefully) ratio of the opposite and adjacent sides of a right angled triangle and how that relates to the angle.

If we set the angle to 45°, which is π/4 radians, then:

   tan(45°) = 1
⇒ tan(π/4) = 1
⇒ arctan(1) = π/4

which when substituted in the above becomes:

\begin{equation*}
\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - \cdots
\end{equation*}

This may be the simplest infinite series formula for π but it has the disadvantage that you'll need to add up a great number of terms to get an acceptably accurate answer.

Here is a bit of python to calculate it (pi_gregory.py):

import math

def pi_gregory(n):
    """
    Calculate n iterations of Gregory's series
    """
    result = 0
    divisor = 1
    sign = 1
    for i in range(n):
        result += sign / divisor
        divisor += 2
        sign = -sign
    return 4 * result

def main():
    """
    Try Gregory's series
    """
    for log_n in range(1,8):
        n = 10**log_n
        result = pi_gregory(n)
        error = result - math.pi
        print("%8d iterations %.10f error %.10f" % (n, result, error))

if __name__ == "__main__":
    main()

All the important stuff happens in pi_gregory and hopefully you'll see how that relates to the infinite series for arctan. Note the trick of keeping the sign of the term in a variable as 1 or -1 and then multiplying by this. This is a common trick, easier to write than an if statement and probably faster.

The program above prints:

Iterations Result Error
10 3.0418396189 -0.0997530347
100 3.1315929036 -0.0099997500
1000 3.1405926538 -0.0009999997
10000 3.1414926536 -0.0001000000
100000 3.1415826536 -0.0000100000
1000000 3.1415916536 -0.0000010000
10000000 3.1415925536 -0.0000001000

That looks like our old friend π doesn't it! However you can see that it takes 10 times as many iterations to add another decimal place of accuracy to the result. That isn't going to get us even 10 decimal places in a reasonable length of time... We'll have to get a bit cleverer. (If the program printed out something very different then you need to run it with python 3 not python 2 see the introduction.) You can also see something very odd in the error terms - the result is correct to 10 decimal places, all except for one digit. This is a known oddity and you can read up about it The Leibniz formula for pi .

Those of you paying attention will have noted that we used math.pi from the python standard library in the above program to calculate the difference of the result from π. So python already knows the value of π! However that is only a double precison value (17 digits or so) and we are aiming for much more. We are going to leave the world of double precision floating point behind, and calculate a lot more digits of π, much quicker, in the next exciting episode (Part 2)!

Categories: maths, pymath, python | View Comments
Read and Post Comments

This is going to be a series of articles which are all about having fun with maths and Python. I'm going to show you a bit of maths and then how to apply it with a python program, then work out how to improve it.

I decided to write these after looking in my doodles folder on my computer and finding 100s of programs to calculate this bit of maths, or explore that bit of maths, and I thought it would be interesting for me to share them. Hopefully you'll find it interesting too!

I'm going to use python 3 for all python code. This enables me to future proof these articles and encourage everyone who hasn't tried Python 3 to have a go. The programs will run with minor modifications on python 2 - just watch out for division which has changed from python 2 to 3. In fact if you run any of these programs under python 2 then put from __future__ import division at the top of them.

Python is a really excellent language for investigations into maths. It is easy to read and has excellent long number support with infinite precision integers and the decimal module. I'm going to attempt to demonstrate what I want to demonstrate not using any third party modules, but I may occasionally use the gmpy module.

I'm going to pitch the level of maths below calculus, no more than the maths that everyone should have learnt at school. If I want to use some bit of maths outside that range I'll explain it first. This does means that sometimes I won't be able to show a proof for things, but I'll add a link to an online proof if I can find one.

You are going to see maths written both in traditional style:

\begin{equation*}
x=\frac{-b \pm \sqrt {b^2-4ac}}{2a}
\end{equation*}

And in computer style which every programmer knows anyway. This involves a few more brackets, but it is much easier to translate to computer code. For example:

(-b +/- sqrt(b**2 - 4*a*c)) / (2*a)

In Python computer maths we write:

  • * to mean multiply
  • / to mean divide
  • ** to mean to the power
  • // to mean divide but give the result as an integer
  • % to mean the same as the modulo or mod operator, ie do the division and take the remainder.

I'm going to make free with unicode fonts so you'll see:

  • π - pi
  • x² - x squared
  • √x - square root of x
  • ± - plus or minus

If these aren't appearing properly, then you'll need to try a different browser, or install some more fonts.

As for the programming level, anyone who knows a little bit of python should be able to work out what is happening - I won't be doing anything too fancy!

I'll tag all the articles with pymath so you can find them or subscribe to the rss.

Happy Calculating!

Nick Craig-Wood

nick@craig-wood.com

Categories: pymath | View Comments
Read and Post Comments

Technology changes in 2010

2010-12-22 14:56

Technology that has come and gone from my life in 2010, and new stuff I've learnt. More than I expected!

These notes are really for me to look back at over the years, but you might find them interesting.

Old

Things gone from my life. Gone but not forgotten :-

  • Mutt - finally replaced with Thunderbird. Mutt is a terminal based email client, just in case you haven't heard of it. It works really well for power users, is really configurable, but in these days of HTML emails and attachments it was no longer keeping up. I loved him - he was a faithful mutt, but he had to go. (Actually he is still there at the command line ready to go if needed - but keep it under your hat!)
  • Courier V34 modem - years overdue, I finally disconnected it and consigned it to the pit of oblivion. It used to receive faxes, but I don't think anyone has sent one for several years. It also used to run a BBS and an Internet node, which was useful at times, but now GPRS/3G/Wifi is everywhere I want to go it is no longer needed.
  • ISDN - 128k of digital goodness just wasn't good enough any more and became too expensive. Great in its hey-day with 2 channels, 10 numbers and a whopping 64k digitial Internet connection. Obsoleted by ADSL and VOIP.
  • Nokia E90 & Symbian - sorry... we've had good times, but I've found someone else
  • KDE - you worked really well until version 4 where the bugs and the slowness really kicked in.
  • Usenet - This was my first hint of online communities on the Internet in 1994. Unfortunately the cool kids have moved on, and so, reluctantly, have I!

New

New technology :-

  • Thunderbird - for people with complicated email lives the plugin system makes life bearable! Here are my top plugins
    • Nostalgy - adds keyboard shortcuts to save, copy, goto folders
    • Quickfolders - Put your favourite folders along the top
    • Copy Sent to Current - copy sent mail to the current folder not the Sent folder - makes a new and better way of doing email
    • Display Mail User Agent - snoop on what other people use for their mail client
    • Enigmail - Use PGP to sign/encrypt your emails seamlessly
    • External Editor - edit emails in your favourite editor (eg emacs)
    • Identity Chooser - If you use multiple identities you need this!
    • Mail Redirect - proper bounce of email (not forward)
  • IMAP - makes it easy to migrate your email setup. You can still use it with procmail so some things haven't changed ;-) I can use it with my Android phone and the excellent K9Mail.
  • Blogofile - which makes this (fully static) website. It replaces a pile of crufty hand built python scripts which did the same thing, just not nearly as well!
  • Ubuntu on the desktop - I finally gave in and replaced Debian testing with Ubuntu on my main laptop so I can spend less time fixing stuff and more time doing real work! Debian still rules on the server though.
  • VOIP - I've signed up with a 100% Voip provider and migrated my 10 telephone numbers to them. It went reasonably smoothly! I took the ISDN card out of my Asterisk PBX. We've been doing this at work for ages - it seemed like time to migrate at home too.
  • Android - is Google's new phone operating system. It runs on Linux and is mostly open source which are big positives in my view. I've been enjoying my Android phone greatly and have written a number of applications - Oxo3d is one.
  • Gnome - well not new actually, I last gave Gnome up in about 2005 when I got irritated at no longer being able to edit my menus. Unfortunately KDE 4 isn't shaping up for me, so it is hello again!

Books

Tech books I've read this year. No python books this year - I think I've read them all now!

  • Java in a Nutshell (again) programming to write Android programs. I learnt original Java ages ago, but is is quite a different experience with generic types. Still not my favourite language - reminds me too much of C++, but bearable.
  • Programming In Scala - just for fun - maybe I'll use it with Android one day. Scala is one of the new breed of languages which runs on the Java VM. Scala is a (mostly) functional language created by one of the creators of Java itself.
  • Version Control with Git - distributed revision control system. Very clever idea with lots and lots of command lines to explore. If it scales to the Linux kernel then surely it is good enough for my projects?
  • BGP: Building Reliable Networks with the Border Gateway Protocol - how the Internet really works!
  • Unlocking Android: A Developer's Guide - how to make Android applications
  • Running IPv6 - yes it is coming - IPv4 addresses are running out very soon.

These are affiliate links, so if you click then buy I'll earn a few pence towards my next tech book!

Categories: tech | View Comments
Read and Post Comments

Oxo3d v0.94 released

2010-11-21 18:48

Oxo3d v0.94 released.

Major changes are

  • Fix missing textures on milestone/droid
Categories: oxo3d, android | View Comments
Read and Post Comments
« Previous Page -- Next Page »