Friday, March 14, 2014

Pi is for Everyone

Source: Wikipedia
The formulas, jargon, and code in this post may not be the best way to prove that STEM is for everyone, but it's Pi Day, and I think it's good to keep in mind that things we take for granted now were once considered extremely difficult.  According to CNN, it's cool to like pi.

So today, in honor of Pi Day, I decided to refresh my memory and sharpen my programming skills by implementing Viète's Forumla. Viète's formula, developed in 1593, expands on the Archimedes method of approximating pi developed around 250 B.C.  This method relies on computing the circumference of regular polygons with many sides. As you increase the sides of the polygon, it gets closer and closer to being a circle, and thus, the ratio of the perimeter to the "diameter" gets closer to the value of pi.

This can be represented mathematically as: 
\pi = \lim_{k\to\infty} 2^k \underbrace{\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\cdots}}}}}}_{k\ \mathrm{square}\ \mathrm{roots}}
Source: Wikipedia
While this can be solved iteratively, implementing this is a good way to get a handle on recursion.  Implementing all of those "2+sqrt(2+sqrt(2+..." terms recursively is pretty straightforward (I call this the "inner term."  After that, you just need to take the square root of 2 minus the inner term, and then multiply that times 2 raised to the power of however many square roots you've taken.

If you understand the code below, you can do many powerful things.  I've written it in perl, so no compiling is necessary and many systems have perl installed already (especially if you have a Mac or Linux).  But you should be able to easily write a version of this program for any programming language.


Some notes/questions about this code:
  1. This program will never stop running on its own.  It is an infinite series after all.  Hit Ctrl-C when you've had enough.  A simple change (adding 3 characters) would make it stop running after 10 of iterations - do you know what those 3 characters are, and where they would go?
  2. I've broken some things out to make it a little more clear to novice programmers what is happening here.  I could eliminate several lines of code at the expense of clarity.  If I were really trying to demonstrate the power of recursion, I would not be storing intermediate values inside the subroutine.  How would you re-write the subroutine to use less memory?
  3. All of this talk of pi is making me hungry.  What kind of pie should I make?
  4. As Larry points out in the comments, this version starts to run into problems due to lack of precision after a fairly small number of rounds.  Below the first version is a version that uses perl's BigFloat module, which does not have this problem (at least not for a much larger number of rounds). 


#!/usr/bin/perl

my $k=1;

while($k++){  

  my $sumterm = sqrt2sum(0$k-1);  #calculate the inner term
  my $pi = 2**$k * sqrt(- $sumterm); # final square root and multiplication

  print "$k: $pi\n";

  sleep 1;  #this slows it down so you can watch it run.  Delete this line for speed
}


sub sqrt2sum {
  my ($val, $k) = @_;
  if($k > 0){
    $val = sqrt(2+$val);
    $k--;
    sqrt2sum($val, $k)
  }else{
    return $val;
  }

}


Version with Precision Issues Fixed:

#!/usr/bin/perl -wT
use strict;
no warnings 'recursion';
use Math::BigFloat;
Math::BigFloat->accuracy(100);

my $k=1;

while($k++ < 150){
  my $sumterm = Math::BigFloat->new(sqrt2sum(0,($k-1)));
  $sumterm = $sumterm->bsub(2) * -1;
  my $root = $sumterm->bsqrt;
  my $pi = $root->bmul(2**$k);

  print "$k: $pi\n";

  #sleep 1;
}


sub sqrt2sum {
  my ($val, $k) = @_;
  my $bigval = Math::BigFloat->new($val);
  if($k > 0){
    $bigval = $bigval->badd(2);
    $bigval = $bigval->bsqrt;
    $k--;
    sqrt2sum($bigval, $k)
  }else{
    return $bigval;
  }


}

4 comments:

  1. I look forward to the followup post that explains why the estimate of pi gets worse and then starts printing 0.

    ReplyDelete
  2. I hadn't ever let it run that long, actually. Are you volunteering to write a guest post on floating point precision for me? This is why I don't do this stuff for a living (and why nobody does this kind of stuff in perl).

    ReplyDelete
    Replies
    1. Using perl's Math::BigFloat module takes care of this problem, but yields much less readable code.

      Delete
    2. Ha, awesome! I can now calculate many more digits of pi.

      Delete