Buffon's Needle

Buffon's Needle is an embarrassingly parallel method for approximating π (that's the greek letter Pi, the ratio of the circumference of a circle and its diameter). It is a good test of whether the Fortress implementation du jour is scaling, or not. This code can also be found at trunk/ProjectFortress/demos/buffons.fss.

Buffon's needle is a parallel Monte Carlo simulation for estimating π. Given a sheet of paper with lines a unit length apart, if we drop a needle of the same unit length onto the paper, the probability of it intersecting a line is 2/π. We can take advantage of this to write a simulation to estimate π.

   probability = hits/n
   pi_est = 2.0/probability
   printTime(6.0)
   println("")
   print("estimated Pi = ")
   println(pi_est)

The "variables" needleLength, numRows, and tableHeight are all immutable; they may not be further assigned to. This is also all the declaration that they need. The variables hits and n may be assigned to, and their declarations require a type to indicate what may be assigned into them. Notice the expression initializing tableHeight; as in mathematical notation, juxtaposition (of non-functions) indicates multiplication.

needleLength = 20
numRows = 10
tableHeight = needleLength numRows

...
   var hits : RR64 := 0.0
   var n : RR64 := 0.0

Choosing (delta_X, delta_Y) randomly in the top half of the unit circle allows us to specify an angle between 0 and π without actually specifying π. Note that we discard any (delta_X, delta_Y) pairs which fall outside the unit circle. If we randomly choose a y1 coordinate, and an angle determined by (delta_X, delta_Y), we can calculate the y2 coordinate. If ceiling(y_L/needleLength) = floor(y_H/needleLength), then the needle crosses a line. Iteration of the for loop is controlled by its generator, and the generator here is parallel, so the loop runs in parallel. If the implementation were more mature, the atomic update at the end of the loop could be replaced with a naked increment, which would automatically be transformed into a reduction, but we're not there yet.

   for i <- 1#3000 do
      delta_X = random(2.0) - 1
      delta_Y = random(2.0) - 1
      rsq = delta_X^2 + delta_Y^2
      if 0 < rsq < 1 then
         y1 = tableHeight random(1.0)
         y2 = y1 + needleLength (delta_Y / SQRT rsq)
         (y_L, y_H) = (y1 MIN y2, y1 MAX y2)
         if |/ y_L/needleLength \| = |\ y_H/needleLength /| then
                        atomic do hits += 1.0 end
         end
         atomic do n += 1.0 end
      end
   end

The trunk/Fortify tool can render the code as follows: