PurposeThe purpose of this project is to explore a Monte Carlo
method of estimating the value of pi (p),
which was first discovered by Count Buffon in 1733.

Mathematics Background
Take a number of toothpicks and randomly throw them on a tablecloth, a hardwood floor, a
brick sidewalk, or on anything that has a number of equally spaced parallel lines. The
parallel lines must be at least as far apart as the length of a toothpick. Count the
number of tosses and the number of toothpicks that cross the parallel lines. The following
formula gives you an estimate for p:

If you can find some parallel lines that are exactly as far apart as the
length of a toothpick, the formula simplifies to:

This is the equation well use for the example in the Figure 1:

Figure 1. Tossing Toothpicks to Estimate the Value of Pi

Figure 1 shows 7 tosses and 4 crossings, the formula tells us:

p» 2 ´ 7 / 4 = 3.5

So our estimate is in the right neighborhood  its close to
3.141592653589793... Can we just keep throwing more toothpicks and get a more
precise answer? Not really. You can use more precision on your calculator, but
youll only ever be accurate to a few decimal digits with this method.
Nevertheless, this "Monte Carlo" method of estimating a value is still useful
for certain kinds of scientific calculations. Molecular chemists can use this method
on a computer instead of beakers to explore how molecules react. Nuclear engineers
studying how a nuclear reactor works can use Monte Carlo calculations to answer questions
that have no other answers.

How can throwing toothpicks on parallel lines possibly be used to estimate p? Can you believe that Count Buffon first experimented with
this method in the 1733 when he dropped needles onto parallel lines?

Consider the geometry of the toothpicks when tossed onto the parallel lines in the
Figure 2:

Figure 2. Geometry of Toothpicks Tossed Onto Parallel Lines

In the diagram above, the parallel lines are all W units apart. The toothpicks
are all L units long. Every time a toothpick is thrown onto the parallel
lines, lets take a look at two values, called D and angleq:

Lets consider how far the midpoint of the toothpick is from the closest
parallel line. Lets call this value D. If the midpoint is exactly on one of the
parallel lines, then D = 0. The toothpick at the left in Figure 2 has D = 0 since its
midpoint is exactly on one of the parallel lines. The value of D for the toothpick at the
right in Figure 2 is not zero. The value of D for this toothpick very close to W/2 since
it is nearly halfway between two of the parallel lines.

The value of D must range between 0 and W/2. The maximum value is W/2 since once the
midpoint is greater than W/2 from one parallel line, it is less than W/2 from one of the
other parallel lines.

The angleq is the angle the toothpick makes with a
horizontal line through its midpoint. This angle must range from 0 to 180 degrees.

The value H shown in Figure 2 can be computed from knowledge of L and q:

H = (L/2) sin q

When D £ H, that is, whenever the distance of
the midpoint of the toothpick to the closest parallel line is less than H, the toothpick
crosses a line. When D > H the toothpick does not cross a line.

As we throw toothpicks, we can form a dot plot of D versus q.
We plot points that cross a line in red and points that do not cross a line in black.
Figure 3 is such a dot plot from a computer simulation of 100,000 tosses:

The top part of Figure 3 explains how to compute the probability of a line
crossing a parallel line. This probability is simply the ratio of the Number
Crossing a Parallel Line to the Number Tossed. This value can be computed by simply
observing the number of toothpicks tossed and the number that cross one of the parallel
lines.

But this ratio of crossings versus tosses is also mathematically the
ratio of the red area above to the area of the rectangle. That is:

The Red area can computed using some calculus:

The Area of the Rectangle is much simpler to compute: Area of
Rectangle = width x height = p x W/2.

If desired, change the number of horizontal lines by changing the "Areas"
spinbox, and then pressing the Reset button.

Press the Toss button with "1 time" selected to see
"needles" tossed one-by-one.

Use the Up/Down arrows to change the Toss "times" from 1 to 10,000,000 by
factors of 10. Try 100,000 times on your processor. (On a 120 MHz Pentium this
may take about 47 seconds. On a 450 MHz Pentium II this may take only about 7
seconds.)

Observe the Graph Tabsheet and the distribution of tosses that resulted in line
crossings.

Experiment with the various Random Direction Generation Methods (see discussion below)
by selecting a Radio Button on the Setup tab and repeating Steps 3-5.

With the Graphics checkbox unchecked, try 10,000,000 tosses. (On a 120
MHz Pentium this may take about 58 seconds. On a 450 MHz Pentium II this may take
about 10.5 seconds.)

DiscussionThe FormCreate method creates the in-memory bitmaps used for the
"toss" and "graph" areas. These bitmaps are only freed when the
program terminates. The FormCreate sets up the TSimplePantograph,
which is a tool for mapping from "real world" coordinates to a particular pixel
in a bitmap. FormCreate also establishes DirectionVectorRoutine
:= DirectionMethod2 (on the Setup tab) as the default method of creating
random directions (discussed more below).

The ResetState method is called on FormCreate and
whenever the Reset button is pressed, or when the Setup algorithm is
changed.. ResetState clears the "toss" and "graph"
areas and draws the axes on the graph. The parallel lines are drawn using the value
specified in the "Areas" TSpinBox. When there are N
areas, there are N+1 lines, including the top and bottom lines of the
"toss" area.

A ProgressBar and Cancel button appears when more than
100 needles are tossed. Pressing the Cancel button sets a flag that is
recognized inside the computation loop. Application.ProcessMessages is called when
each 1% of the calculations have been completed. (Actually, it's called after each
0.2% of calculations are completed when more than 1,000,000 tosses are made.) The
computations can be completed much more quickly by unchecking the Graphics checkbox.

The ButtonTossClick method is called when the Toss
button is pressed. This method loops for the specified number of times and calls the
TossBuffonNeedle method.

The TossBuffonNeedle procedure assumes L = W = 1 to avoid
needless computations. The value of D, which is the distance from a parallel line to
the midpoint of the needle, is computed by:

D := Random - 0.5;

D ranges from -0.5 (below a line) to 0.5 (above a line). The random
orientation angle, q, ranges from 0 to 180 degrees (p radians). q is defined by calling
the DirectionVectorRoutine, which is a procedure variable that can have a value
of Method1, Method2, Method3, or Method 4. These direction methods
will be discussed later, but each routine defines the values x, y, Sin(q),
and Cos(q). A random direction vector is from (0,0) to
(x,y) and forms an angle q with the positive X axis.
(See the diagrams below explaining the various random direction generation
methods.)

The Cos(q) value is only computed to
display the "toss." Only Sin(q)
is needed to compute whether a needle crosses a line:

Using the Pixels property (instead of Scanline) for setting the red
and black dots in the graph using an in-memory bitmap wasn't that slow in this program.
Once a few hundred thousand needles are "tossed," every pixel is covered
and the graphical image doesn't change that much anyway.

The in-memory BitmapTosses and BitmapGraph are assigned in ButtonTossClick
to update the display:

Now that the assignment of the DirectionVectorRoutine has been discussed,
let's look at each of the Random Direction Generation Methods. Remember that since
we are trying to compute p that none of the algorithms
can depend directly on p.

Method 1. Naive and Incorrect Approach
This was my original, naive, and INCORRECT approach to create a random angle --
it's still educational to consider.

A random (x,y) value is chosen in the rectangle ranging from -1 to 1
in the X direction and 0 to 1 in the Y direction.

The Sin and Cos of the angle, q, between the
X-axis and the line from (0,0) to (x,y) were easily computed.

This algorithm is simple but is not correct. Any of the (x,y) points
along the line shown below will result in the same angle:

The lines from the origin (0,0) to the bounding rectangle as a function of
angle are not all the same length. The resulting graph shows a higher relative
probability near p/4 (45 degrees) and 3p/4
(135 degrees).

A plot of D versus q using Method 1
shows this bias:

Because of this minor bias, the estimates of p
using Method 1 are usually lower than the true value.

Method 1 can be modified, as suggested by Gustavo Massaccesi from
Argentina, to correct for this uneven probability distribution by angle.

Method 2. Gustavo Massaccesi's "Fix" to Method 1
My naive approach can be fixed by requiring each (x, y) point to be within the unit
semicircle.

The extra check involves using the Pythagorean formula to calculate the
distance and a REPEAT .. UNTIL loop:

Method 3. Ingenieurbüro Dr. Elsing's
Random Arc Length MethodIngenieurbüro Dr. Rainer Elsing from Germany suggested a different
approach to find a uniformly distributed random angle. In finding a random angle, we
pick a random arc length that is between 0 and an assumed large value that is much greater
than the circumference of the circle. Because the Sin and Cos functions are
periodic, the Sin and Cos of this random arc length can be used to calculate the
(x,y) point like that used in Methods 1 and 2.

This is the only method that directly computes Sin and Cos values for
arbitrary angles.

Method 4. Division of Semicircle into
2N Intervals Using Half-Angle Formula and Rotation Formula
An easy way to get a random direction vector, and simultaneously the Sin and Cos values,
would be to pick a point at random from a set of equally-spaced points along the
semicircle described in Method 2.

Such a set of equally-spaced points along a semicircle can be found in two
steps (and without direct evaluation of Sin and Cos for arbitrary angles):

1. Recall the Half-Angle Trig Formulas:

Using these formulas, we start at 90 degrees with Cos(90 degrees) = 0 and
compute the Sin and Cos of a very small angle.

As an example, let's assume we want to divide the semicircle into 23
= 8 intervals (the same logic works for any power of 2). This means we first
need the Sin and Cos of 180/8 = 22.5 degrees. To compute this, we start with Cos(90
degrees) = 0, and apply the half-angle formulas twice:

Angle [degrees]

Cos(angle)

Sin(angle)

90.0

0.0

not needed

45.0

0.707107

0.707107

22.5

0.923880

0.382683

The diagram below shows the geometry of this final angle in this example:

2. Now starting at the point (1, 0) on the unit circle, rotate this
point (1,0) about the origin by the angle in Step 1 using the rotation matrix.

Repeat the rotation process until the desired number of points have been
computed.

Here a summary table and diagram for this 23-interval case:

Step

x

y

0

1.000000

0.000000

1

0.923880

0.382683

2

0.707107

0.707107

3

0.382683

0.923880

4

0.000000

1.000000

5

-0.382683

0.923880

6

-0.707107

0.707107

7

-0.923880

0.382683

8

-1.000000

0.000000

With the table of (x,y) values also being the Sin and Cos values, the DirectionMethod4
routine was quite simple:

// Since x and y were defined on unit semicircle, they are also the
// SinTheta and CosTheta values
CosTheta := x;
SinTheta := y
END {DirectionMethod4};

The xArray and yArray are assigned values in the ScreenBuffon
unit initialization. The current implementation has . LookupArrayCount =
4096. See the source code for details.

ConclusionsBuffon's method is an interesting way to compute p
using random numbers and "tossing" needles onto parallel lines.

The execution times when graphics is enabled is dominated by the graphics
drawing time instead of the algorithm itself.

When graphics are disabled, the algorithms dominate the execution times.
Surprisingly, the Sin and Cos functions evaluated in Method 3 were fairly quick on
a Pentium. Method 4 is the fastest method since it uses a lookup table and did no
other computations.

Method 4 appears to be the most accurate method, but I'm not sure how one
proves accuracy in a Monte Carlo simulation.

Method 1 with its non-uniform angular probability distribution is really
solving a slightly different problem than the one desired. Method 1 predictably
always had a lower value for p, but was still in the right
"neighborhood."

Acknowledgements
Thanks to Ingenieurbüro Dr. Rainer Elsing (home page in German
or English) and Gustavo Massaccesi (home
page in Spanish) for their assistance in refining random direction algorithms.
Thanks to Peter Hamer for some very useful links about Buffon's Needle.