• Coding
  • [Exercise] Approximating Pi

Almost 400 years ago, mathematician John Wallis has discovered the follwing formula:
pi/4 = 2/3 * 4/3 * 4/5 * 6/5 * 6/7 * 8/7 * ...
The exercise consists of calculating the value of Pi (to an arbitrary order of precision) by using the above formula.

Extra Credits:

The exercise is taken from SICP, and is meant to introduce the notion of higher order functions.
A good way to think of the solution would be to consider a function product that takes a function f as a parameter and computes the product of f(n) for n in [0; n] (or [a; b]).
g[n_] := 2 Quotient[n, 2] + 2;
f[n_] := g[n + 1] / (g[n] + 1);
pi[n_] := 4 Product[f[i], {i, 0, n}];
f /@ Range[0, 10]
pi /@ Table[10^i, {i, 6}] // N
{2/3, 4/3, 4/5, 6/5, 6/7, 8/7, 8/9, 10/9, 10/11, 12/11, 12/13}
{3.02317, 3.12638, 3.14003, 3.14144, 3.14158, 3.14159}
#include <iostream>

using namespace std;

int main()
{
    float pi=1;
    int a=2,b=3, N;
    cin >> N;
    for(int i=0;i<N;i++)
    {
        pi*=a*1.0/b;
        if(i%2==0)
        {
            a+=2;
        }
        else
        {
            b+=2;
        }
    }
    cout << pi*4;
    return 0;
}
Monte Carlo exercise
@Ra8: cool code!
You should edit your answer and post the Monte Carlo exercise as its own thread.

Here's my code in Scheme:
(define product 
  (lambda (term next prod a b)
    (if (> a b)
	prod
	(product term
		   next
		   (* prod (term a))
		   (next a)
		   b))))


(define even?
  (lambda (x) (= (remainder x 2) 0)))

(define approx-pi
  (lambda (n)
    (* 4 
       (product (lambda (x) (if (even? x)
			   (/ (+ x 2.0) (+ x 1))
			   (/ (+ x 1.0) (+ x 2)))
		    (lambda (x) (+ x 1))
		    1
		    1
		    n))))

(write (approx-pi 1000))
(write (approx-pi 10000))
(write (approx-pi 100000))
Result:
3.14316070553226
3.14174970573801
3.14160836127794
4 months later
Here's another solution taking advantage of python's generators and itertools:
import itertools as it
from operator import mul

def foo():
    n = 0
    a = 2.0
    b = 3.0

    while True:
        yield a/b
        n += 1
        if n%2:
            a += 2
        else:
            b += 2      

            
a = it.takewhile(lambda x: abs(1-x) > 1e-6, foo())

# Python badly needs a 'product()' function!
print 4*reduce(mul, a)
a year later
hey guys it has been a while since i last checked on you
It seems you're fine.
anyways here's what i got to, tell me if any good.. also tips are appreciated.
from fractions import Fraction
def calcPiToX(X): 
  listx = []
  for i in xrange(4,X+1):
    if i%2==0:
      listx.append(Fraction(i,i-1))
      listx.append(Fraction(i,i+1))
  listx.append(Fraction(2,3))
  return 4*reduce(lambda x,y: float(x)*float(y),map(float, listx))
@NuclearVision

I'm not going to comment on the exercise, just the very last line of your code:
return 4*reduce(lambda x,y: float(x)*float(y),map(float, listx))
You're mapping all the entries in listx to float, then in your reduce's applied function, you also float x and y. That's redundant.
return 4*reduce(lambda x,y: x*y, map(float, listx))
Ofcourse you can act all cool and replace lambda x,y: x*y with operator.mul You'll be left with:
import operator
...
return 4*reduce(operator.mul, map(float, listx))
@xterm thanks for the float correction.
as for the lambda function i wanted to make use of it xD. And it is indeed cooler and shorter than importing a whole module.
As for the answers im getting ones different from geeks's.
can any one define the 'order of precision' is it the n-th fraction?
thanks.
Style remarks
Some minor remarks about your style:
  • Use 4 spaces for indentation
  • Do not use camelCase for function names. Prefer underscores. A better name would be approximate_pi.
  • Insert more empty lines, in particular 2 empty lines after the import statement.
  • Spaces around operators and between function arguments.
  • Avoid single letter variable names as much as you can. When dealing with mathematical problems it's not always easy to do, so it's okay to keep them. However avoid upper cases, unless you have reasons not to.
Look at how style can make your code more readable.
from fractions import Fraction


def approximate_pi(n):
    listx = []

    for i in xrange(4, n+1):
        if i % 2 == 0:
            listx.append(Fraction(i, i-1))
            listx.append(Fraction(i, i+1))

    listx.append(Fraction(2, 3))
    return 4 * reduce(lambda x, y: float(x) * float(y), map(float, listx))
To understand better why style is such a big deal in the Python community, and how you should format your code, you should read PEP8.

No need for a list
The problem with your code is that you're building a useless list. This can be costly if you want to run it on big numbers, affecting both memory consumption and performance.

You can solve the problem with a generator.
from fractions import Fraction
from operator import mul

def approximate_pi_generator(n):
    yield Fraction(2, 3)

    for i in xrange(4, n+1):
        if i % 2 == 0:
            yield Fraction(i, i-1)
            yield Fraction(i, i+1)

def approximate_pi(n):
    return 4 * reduce(mul, map(float, approximate_pi_generator(n)))
If you don't really understand generators, or how the yield keyword behaves, you should read this StackOverflow answer that does a great job at explaining them.

No need for a container
You don't even need to hold all these fraction in one place. You can multiply the Fractions on the fly.
from fractions import Fraction

def approximate_pi(n):
    pi = float(Fraction(2, 3))

    for i in xrange(4, n+1):
        if i % 2 == 0:
            pi *= Fraction(i, i-1)
            pi *= Fraction(i, i+1)

    return 4 * pi
Notice that I just explicitly transform the first element into a float. All the subsequent Fraction objects are automatically coerced into floats by Python itself.

Note on functional programming
While map, lambda and reduce are fun, you should avoid using them when you can. Notice how the last code I write is more readable and quickly understandable.

You can read more about Guido's dislike of functional programming features in these links:
Thanks for that reply, I will try to remember those in the future.
and by the way rahmu I see you beat your old code :-)
6 months later
Simple way in C#.
        public static double getpi(int precision)
        {
            double _pi = (2d / 3d), _numrt = 2d, _denum = 3d;

            for (int i = 0; i < precision; i++)
            {
                _numrt += i % 2 == 0 ? 2 : 0;
                _denum += (i - 1) % 2 == 0 ? 2 : 0;

                _pi *= (_numrt / _denum);
            }

            return _pi * 4d;
        }
Usually have precision > 1000 so it can give at least the first few values correctly.
4*(n-1)*reduce(lambda x,y: x*y, [float(i**2)/(i+1)**2 for i in xrange(2,n,2)] )/2 #n%2!=0
One line!
20 days later
#include <iostream>
using namespace std;
long double picalc(int n)
{
  long double j=1,two=2,k=2;
  for(j;j<=n;j++) k*=(two*j*two*j/(two*two*j*j-1.0));
  return k;
}
int main()
{ cout.precision(71);
  cout<<picalc(1000)<<endl;
  return 0;
}