Monte Carlo Method and pi
beliavsky at aol.com
beliavsky at aol.com
Mon Jul 12 19:54:56 CEST 2004
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
More information about the Python-list
mailing list