Generally speaking, the Monte Carlo method refers to a technique that attempts to find an approximate solution to a numerical problem by making use of randomly generated inputs. I would like to present as follows an implementation based on such a technique that approximates the value of the constant pi (the ratio of the circumference of circle over its diameter).
The meaning of the following declarations should be clear:
// #define N 1024 // typedef point = $tup(double, double) // typedef points = list0(point) typedef pointss = list0(points) typedef pointsss = list0(pointss) //
Assume that there are K points on the grid. A randomly chosen point generates one hit if it resides within the unit disk centered at one point among the K points. Accumulatively, a randomly chosen point can generate n hits if it resides within the unit disks centered at n points among the K points.
Let create a matrix of the dimension N by N for storing lists of points:
We are to play N2 rounds of experiment. At round K for each K between 1 and N2, there are already K-1 points on the grid chosen randomly in the previous rounds. We choose a new point on the grid randomly for round K and records the number of hits it generates (with respect to the K-1 points on the grid). We then add this new point to the list stored at the matrix cell (i, j), where i and j are the integer parts of the x-coordinate and y-coordinate, respectively.The following functions are for computing the number of hits generated by a point with respect to a collection of points:
// fun hit_point_point ( p1: point, p2: point ) : int = ( if dist(p1, p2) <= 1.0 then 1 else 0 ) // fun hit_point_points ( p0: point, ps: points ) : int = ( list0_foldleft<int><point> ( ps, 0 , lam(res, p1) => res + hit_point_point(p0, p1) ) (* list0_foldleft *) ) // fun hit_point_pointss ( p0: point, pss: pointss ) : int = ( list0_foldleft<int><points> ( pss, 0 , lam(res, ps) => res + hit_point_points(p0, ps) ) (* list0_foldleft *) ) // fun hit_point_pointsss ( p0: point, psss: pointsss ) : int = ( list0_foldleft<int><pointss> ( psss, 0 , lam(res, ps) => res + hit_point_pointss(p0, ps) ) (* list0_foldleft *) ) //
A unit square in the grid is located at (i, j) if the lower-left corner of the unit square is at the point (i, j). Two unit squares are neighbors if they locate at (i1, j1) and (i2, j2) such that abs(i1-i2) and abs(j1-j2) are at most 1, where abs refers to the absolute value function (on integers). Given a point inside the unit square located at (i, j), the point can only generate hits with respect to those residing inside the neighbors of the unit square. Given a unit square at (i, j), the following function generate a collection of the entire neighbors of the unit square:
// fun theGrid_get_neighbors ( i: int, j: int ) : pointsss = let // fun isvalid_row (i: int): bool = (0 <= i && i < N) fun isvalid_col (j: int): bool = (0 <= j && j < N) // fun fopr ( i: int , j: int ) : points = if ( isvalid_row(i) && isvalid_col(j) ) then theGrid[i, j] else list0_nil() // in // int_list0_map<pointss> ( 3 , lam(i') => int_list0_map<points> ( 3 , lam(j') => fopr(i + i' - 1, j + j' - 1) ) ) (* end of [int_list0_map] *) // end // end of [theGrid_get_neighbors] //
// fun do_one(): int = let // val x = N * randfloat() val y = N * randfloat() // val p = $tup(x, y) // val i = g0float2int_double_int(x) val j = g0float2int_double_int(y) // val nhit = hit_point_pointsss (p, theGrid_get_neighbors(i, j)) // in // theGrid[i, j] := list0_cons(p, theGrid[i, j]); nhit // end // end of [do_one] //
// fun do_all ( K: int ) : int = ( fix loop ( i: int, res: int ) : int =<cloref1> if i < K then loop(i+1, do_one()+res) else res )(0, 0) // end of [do_all] //
It can be proven that do_all(N*N)/(N*N) converges to pi/2 (with probability 1) as N approaches the positive infinity. Please find on-line the entirety of the code used in this chapter.