Monte Carlo Method and pi

[beliavsky at aol.com]

Dan Christensen <jdc at uwo.ca> wrote in message news:<87llhqbbh5.fsf at uwo.ca>...
> 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.)

For comparison, on a 2.8GHz P4 the Fortran code to compute pi using
simulation runs in about 0.09s (measured using timethis.exe on Windows
XP), REGARDLESS of whether the calcuation is done with a loop or an
array. If 5e6 instead of 5e5 random numbers are used, the time taken
is about 0.4s, and the scaling with #iterations is about linear beyond
that. Here are the codes, which were compiled with the -optimize:4
option using Compaq Visual Fortran 6.6c.

program xpi_sim 
! compute pi with simulation, using an array
implicit none
integer, parameter      :: n = 500000
real (kind=kind(1.0d0)) :: xx(n),pi
call random_seed()
call random_number(xx)
pi = 4*sum(sqrt(1-xx**2))/n
print*,pi
end program xpi_sim

program xpi_sim 
! compute pi with simulation, using a loop
implicit none
integer                 :: i
integer, parameter      :: n = 500000
real (kind=kind(1.0d0)) :: xx,pi
call random_seed()
do i=1,n
   call random_number(xx)
   pi = pi + sqrt(1-xx**2)
end do
pi = 4*pi/n
print*,pi
end program xpi_sim