Whinings on BSD, virtualization, and other stuff.

Approximate the value of pi! (real slowly)

Even though an exact value of pi may never be found, there are many ways you can try and approximate it. One way is to simply draw a circle as perfectly as you can, then divide the circumference by the diameter, pi = C/d, while using the best precision available to you. There is a neat visualization of this on the pi wikipedia entry. The better the precision you use in your measurements, the better approximation of pi you get. Out of the many formulas mathmaticians have used, one stands out as using "simple" mathematics (multiplication, division, addition, and subtraction): Gregory-Leibniz Series.

The approximation can be simplified to:

pi = 4/1 - 4/3 + 4/5 - 4/7 + 4/9 - 4/11 + 4/13 - 4/15 ...

This turns out to be easily translated into code. Back when I first started learning shell scripting on UNIX, I decided to dust off my TI-83 and turn my old TI-BASIC program (not an advised platform) into a neat shell program utilizing bc(1) and some other script-fu to do just this. It was a nice exercise in learning shell scripting, but not very useful for pretty much anything, even testing CPU speeds, as this kind of thing is usually capped by the operating system. I recently wanted to create a program to do exactly that, though, so I decided to revisit the old leibniz shell script. I ended up making a quick and dirty shell script for FreeBSD that abuses the service servicename start rc.d process and writes all of it's variables to a file. There isn't much documentation for this script (because of the CPU limitations of /bin/sh) but the heart of the script is:

# Set default leibniz folder location

i=$1            # Take from input
scale=64        # How far out to take decimal place
seed=$(cat $dir/seed)   # Seed to know where it left off

# Start calculation
while [ $i -gt 0 ];
    # Run new work unit
    pi1=$(echo "scale=$scale; (4/$seed)-(4/($seed+2)) \
    +(4/($seed+4))-(4/($seed+6))" | bc)

    # Set seed for next run
    seed=$(expr $seed + 8)
    echo $seed > $dir/seed

    # Find old value for pi and add to new workunit
    pi2=$(cat $dir/pi)
    pi=$(echo "scale=$scale; $pi1 + $pi2" | bc)
    echo $pi > $dir/pi

    # Iteration increment
    i=$(expr $i - 1)

As you can see from the code, the scale is hardset to 64, meaning this calculation will tell bc(1) to chop everything after the 64th decimal place off. The higher the number, the higher the precision. You pass the iterations at the command line like leibniz 64 which tells the script how many times you want to run a given workunit. Each workunit is basically calculating (4/s)-(4/(s+2))+(4/(s+4))-(4/(s+6) where s is the seed of the workunit, often 1 when starting from scratch. The shell script then goes on to find the last approximation of pi, and adds it to the latest workunit. A new seed is written for the next workunit, a new approximation of pi is written, and the iteration function is ran to prevent infinite looping.

That's a fine and dandy proof of concept, but it does not work very well if you want to tax the CPU, like I wanted to do when testing the clustering capabilities of the data structure store redis. I then picked up the python scripting language again and hacked out a python version of the series that is still a work in progress. Leibniz-py can't distribute workunits to more nodes just yet, but that will be the end goal, having the ability to cluster the datastore nodes and worker nodes. That's for another time, for now, take a look at the python script:

# Import stuff
import redis
from decimal import *
# Create connection to localhost redis server
r = redis.StrictRedis(host='localhost', port=6379, db=0)
# Grab beginning variables
i = r.get('iterations')
scale = r.get('scale')
s = Decimal(r.get('seed'))
# Set the correct scale
getcontext().prec = int(scale)
# Start calculation
while (i > 0):
    pi1 = (Decimal(4) / Decimal(s)) - (Decimal(4) / Decimal(s + 2))
    pi2 = (Decimal(4) / Decimal(s + 4)) - (Decimal(4) / Decimal(s + 6))
    pi3 = Decimal(pi1) + Decimal(pi2)
    pi4 = Decimal(r.get('pi'))
    pi5 = Decimal(pi3) + Decimal(pi4)
    r.set('pi', pi5)
    # Get seed ready for next workunit
    s = s + 8
    r.set('seed', s)
    # Iteration increment
    i = int(i) - 1
# Print the latest approximation when done
lpi = Decimal(r.get('pi'))
print lpi

This version is not optimized to use redis in a very efficient matter, as we can clean up some cycles using pipelining. The focus of this is to break down the script to a basic work unit that can be used later on in a distributed fashion. As you can see, this script makes use of the Decimal() function where you can use numbers with large amounts of data after the decimal point. You can see the scale is set with getcontext().prec = int(scale) which is 64 by default (set in leibniz-reset). The iterations are set very low by default, you can actually bump that up to the millions on descent hardware and not have to wait forever for the run to complete.

This proof of concept, at the time of this writing, is fairly portable. It depends on python, redis, and a way for python to talk to redis, redis-py. Installing redis will take some minimal set up on most operating systems, and installing redis-py is as easy as running pip install redis if you have pip installed. Once everything is set up, you must run leibniz-reset to setup the database for leibniz. You can also edit variables in leibniz-reset in a text editor if you choose. For instance, I edit my leibniz-reset so that it runs the workunit ten million times so it will output pi correctly to 7 places after the decimal. I simply change the line r.set('iterations', '16') to r.set('iterations', '10000000') in lebniz-reset. Note this will take some time to run even on modest hardware. You can also run redis-cli in another terminal to check the current status of the approximation, or set a new iteration variable directly like:

> get pi

>set iterations 10000000

Receive Updates