Mathematics

Buffon's Needles

Lab Report

1,000 Needles Tossed Onto Parallel Lines Yield Approximation of p (Using Method 4 -- see text)

Purpose
The 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 angle q:

1. 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.
   
2. The angle q 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:

Figure 3.  Plot of D versus q for 100,000 Tosses
(Using Method 3 -- see below)

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.

Rearrangement of this equation gives us the equation we want:

Or, when L = W:

See "The Monte Carlo Method to approximate Pi " for a different Monte Carlo method for calculating p:  www.daimi.aau.dk/~u951581/pi/MonteCarlo/pimc.html .  (For additional links about p see efg's Math page).

Materials and Equipment

Software Requirements
Windows 95/98/NT
Delphi 3 or 4 (to recompile)
Buffon.EXE program

Hardware Requirements
VGA display

Procedure

1. Double click on the Buffon.EXE icon to start the program.
   
2. If desired, change the number of horizontal lines by changing the "Areas" spinbox, and then pressing the Reset button.
   
3. Press the Toss button with "1 time" selected to see "needles" tossed one-by-one.
   
4. 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.)
   
5. Observe the Graph Tabsheet and the distribution of tosses that resulted in line crossings.
   
6. 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.
   
7. 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.)

Results

Estimates of p for Various Random Direction Generation Methods
(166 MHz Pentium)

Discussion
The 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:

Once the D and q values are computed, and the statistics are updated, the "toss" is moved to a random location on the display.  First the midpoint in the display space is determined randomly:

With this midpoint location, the "needle" is drawn with a random color (except white):

To update the "dot" graph, the ArcTan function is used to find the angle and the SimplePantograph is used to map the (x, y) real-world coordinates to (i, j) pixel coordinates:

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:

The estimate for p is updated in the UpdatePIvalue method:

Note the Plural string function that displays the correct words:  "toss" or "tosses" and "crossing" or "crossings."

Methods of Creating Random Direction Vectors

The TDirectionVectorRoutine type is defined in the ScreenBuffon unit so that a variety of methods can be explored to compute the random direction angle.  

The TFormBuffon class includes a private DirectionVectorRoutine variable

As mentioned earlier, the TFormBuffon.FormCreate routine assigns the default DirectionMethodRoutine:

This assignment of a DirectionVectorRoutine is a much more elegant way in Delphi than the syntax used with function pointers in C/C++.

To force consistency the RadioButtonMethod2 RadioButton is Checked when the DirectionVectorRoutine assignment is made. 

The only other way to change the DirectionMethodRoutine is by directly clicking on a RadioButton on the Setup tab.

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:

A plot of D versus q using Method 2 shows no angular bias.

Method 3.  Ingenieurbüro Dr. Elsing's Random Arc Length Method
Ingenieurbü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 2^N 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 2³ = 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:

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 2³-interval case:

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

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

Conclusions
Buffon'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."

Keywords
pi, Monte Carlo method, Count Buffon, Random, Plural, TBitmap, TImage, Pixels, MoveTo, LineTo, FillRect, TextOut,  TSimplePantograph, ShellExecute, Application.ProcessMessages, Cancel button, FormatFloat, ArcTan2, RGB,   MessageBeep(MB_ICONEXCLAMATION), Procedure variable, Half-Angle Formula, Rotation Matrix

References
Beckmann, Peter Beckmann, A History of Pi, St. Martin's Press, 1971.  pp. 159-161. 

Hamming, Richard W.  Introduction to Applied Numerical Analysis, New York: McGraw Hill, 1971, pp. 311-312.

Ralston, Anthony Ralston (editor).  Encyclopedia of Computer Science and Engineering (second edition), New York: Van Nostrand Reinhold, 1983, pp. 997-998.

Zocchi, F., "Accuracy of Buffon's 200-year-old experimental data," Nature, 1988, 336:318.

Links

Buffon's Needle:  An Analysis and Simulation
www.mste.uiuc.edu/reese/buffon/buffon.html

Lazzarini Eats Humble Pi (Posthumously)
www.knowledge.co.uk/frontiers/sf096/sf096m18.htm

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.

Files
Delphi 3/4 Source (155 KB):  Buffon.Zip

Updated 18 Feb 2002

since 28 Mar 1999