Monte Carlo Method and pi

[Dan Christensen]

I got to thinking about this problem.  The most direct implementation is
something like:

from random import random
from math import sqrt

def v1(n = 500000):
    rand = random
    sqr  = sqrt
    sm = 0
    for i in xrange(n):
        sm += sqr(1-rand()**2)
    return 4*sm/n

This takes about 0.68s on a 2.66GHz P4 (using Python 2.4a) and
about 0.29s with psyco (using Python 2.3.4).

One of the great things about python is that it's possible (and fun)
to move loops into the C level, e.g. with this implementation:

def v2(n = 500000):
    a = numarray.random_array.uniform(0, 1, n)
    return 4 * numarray.sum(numarray.sqrt(1 - a**2)) / n

which clocks in at 0.16s with Python 2.4a.  (psyco doesn't help.)

The problem with the numarray approach is that it consumes large 
amounts of memory for no good reason.  So, now that generator
expressions have arrived, the next thing to try is:

def v6(n = 500000):
    rand = random
    sqr  = sqrt
    sm = sum(sqr(1-rand()**2) for i in xrange(n))
    return 4*sm/n

Unfortunately, this takes 0.64s, not much better than the most direct
method.  (I don't have psyco installed for 2.4, so I haven't tried
that, but I wouldn't be surprised if it was worse than the direct
method with psyco.)

Maybe the reason is that the loop is still in python?  I tried some
awkward ways of eliminating the for loop, e.g.

from itertools import *
from operator import pow, sub, add

def v4(n = 500000):
    rand = random
    pw   = pow
    sqr  = sqrt
    sm  = sum(imap(lambda dummy: sqr(1-rand()**2), repeat(None, n)))
    return 4*sm/n

def v3(n = 500000):
    rand = random
    pw   = pow
    sqr  = sqrt
    sm  = sum(imap(sqr,
                   imap(sub, repeat(1),
                        imap(pw,
                             imap(lambda dummy: rand(), repeat(None, n)),
                             repeat(2)))))
    return 4*sm/n

but they are both slow (0.88s and 0.97s respectively).  

Is there a faster way to get Python to compute a sequence of repeated
operations like this and sum them up?

Is there a more direct way to make an iterator like this?

    itrand = imap(lambda dummy: rand(), repeat(None, n)) 

Wouldn't it be great if one could write iterator code as concisely as
the numarray code, e.g.

    sm = sum(sqrt(1-itrand**2)) 

i.e. if operations like +, - and ** were overloaded to handle
iterators producing numbers ("numiterators"?)  Could a method like
this be fast?

Dan