Code Blog: Ramblings about random bits of code from different languages
Contents: Code Blog Home | Tags | Featured Posts
Thursday, November 26, 2009
Computing The Decimal Places of π
One way of calculating π is continued fractions. One continued fraction which gives π is
So basically a quarter of π is less than 1. The amount that it is less than one is between 1/3 and 1/5 and is found by iterating the denominator of these fractions by two and negating them each time (minus a third, add a fifth, minus a seventh, add a ninth, etc.) - the more iterations of this you do, the more accurate your approximation of 1/4 π becomes, then just multiply this number by 4 to get an approximation of π.
This continued fraction can also be written as an infinite series, known as the
Leibniz formula for π:
I started by writing down the basis of the continued fraction on paper and thinking about how it would work. The easiest way was obvious - a
for loop. I needed to run a loop a LOT of times and needed to think about how I was going to do the alternating negation. Since I was only using the odd numbers from 3 onwards, I set the loop to start at 3 and do steps of 2. I thought about what the relationship between the denominators of the fractions that were added compared to that of those that were subtracted and realised that this solution lied in the remainder of the number when divided by four:{D-} Denominators of subtracted fractions: {3,7,11,...}
{D+} Denominators of added fractions: {5,9,13,...}
All elements of {D-} have a remainder of 3 when divided by 4.
All elements of {D+} have a remainder of 1 when divided by 4.
Finding the remainder of an integer when divided by a certain number is a very useful programming method known as modulo, and is usually denoted by %, e.g. 4%2 = 0 (the remainder of 4/2 is 0); 5%2 = 1 (the remainder of 5/2 is 1). In MATLAB, the function syntax is actually
rem(x,y) (remainder when x/y) as opposed to the % used in other languages such as Java and PHP.So I could use this in my π program. The loop iteration variable would always be the denominator and I could tell it whether to add the fraction or subtract it based on whether the remainder when divided by 4 was 1 or 3.
So now I had the program worked out in my head, it seemed rather simple. I typed it out:
p = 1;
for i = 3:2:10000000000000
if rem(i,4) == 3
p = p - 1/i;
else
p = p + 1/i;
end
fprintf('%10i Iterations: Pi is approximately equal to %1.16f...\n',i,4*p)
end
p as the variable for π as pi is a pre-defined function (equal to π. oddly enough) in MATLAB. The fprintf statement prints the number of iterations (for indication purposes) and the current calculated π value to 16 decimal places. This is printed at every single iteration so you can see live, in real time, how each decimal place is converging, and how many iterations it takes to get each one. The if statement checked to see if the iteration variable i gave a remainder of . The equation started with p=1 and then with each iteration, added (or subtracted) a fraction of 1 over the iteration number (denoted by i) which goes from 3, in steps of 2, up to 10,000,000,000,000 (ten trillion - a very large number just so it would keep running constantly for a long period of time, plenty long enough to see the first few decimal places converge.The output screen would show:
1 Iterations: Pi is approximately equal to 2.6666666666666670...2 Iterations: Pi is approximately equal to 3.4666666666666668...
3 Iterations: Pi is approximately equal to 2.8952380952380956...
4 Iterations: Pi is approximately equal to 3.3396825396825403......and continue until the iteration got to 10 trillion (which would take an awful long time with a home PC!) or I stopped it running!
The value for π would converge one decimal place at a time, so it would fluctuate between 3.0 and 3.2 (with each fraction being added and subtracted - and each fraction getting smaller and less significant) until the fractions being added and subtracted were so insignificant that they did not affect the first decimal place. Then some time later it would converge to 3.14, then 3.141, etc.
It took 626 iterations to get it to converge to 2 D.P. (3.14):
620 Iterations: Pi is approximately equal to 3.1432029585040153...621 Iterations: Pi is approximately equal to 3.1399849375868794...
622 Iterations: Pi is approximately equal to 3.1431977889925018...
623 Iterations: Pi is approximately equal to 3.1399900905161586...
624 Iterations: Pi is approximately equal to 3.1431926525657983...
625 Iterations: Pi is approximately equal to 3.1399952105194355...
626 Iterations: Pi is approximately equal to 3.1431875489073047...
627 Iterations: Pi is approximately equal to 3.1400002979112887...
628 Iterations: Pi is approximately equal to 3.1431824777044470...... 2,453 to get three decimal places (3.141), 136,115 to get to four (3.1415) and about half a million to get five (3.14159).
EDIT: I've now found a way to make the code run more efficiently, without the
if statement.I like the way I thought about the original method so I'm leaving it up, but here's a short explanation of the new one:
Instead of checking to see if the divisor is in the set {D+} or {D-}, it's easier just to put -1 to the power of the iteration number, which will be wither 1 or -1 depending on whether it is even or not, then using this as the numerator of the aforementioned Leibniz formula for π, the denominator being double the iteration number plus one:
And then, as this generated the value for a quarter of π, instead of multiplying this value by four every time, I just added a multiplier of 4 to the numerator of the fraction:p = p + ((-1)^i)) / (2*i + 1)
The final program being:p = p + (4*((-1)^i)) / (2*i + 1)
p = 0;
for i = 0:1000000000000
fprintf('%i Iterations: Pi is approximately equal to %1.16f...\n',i,p)
end
The final program (.m file) can be downloaded here: pi.m
Labels: For Loop, Maths, MATLAB, Pi
Posted by Ben Nuttall at 9:21 PM 
