Buffon's Needles | Lab Report |

**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:

- 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
*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 RequirementsWindows 95/98/NT

Delphi 3 or 4 (to recompile)

Buffon.EXE program

Hardware RequirementsVGA display

**Procedure**

- Double click on the
*Buffon.EXE*icon to start the program.

- 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.)

**Results**

**Estimates of p for Various Random
Direction Generation Methods**

(166 MHz Pentium)

Graphics Checked (Enabled) |
||||

Method |
1,000 tosses |
Time [sec] |
100,000 tosses |
Time [sec] |

1 | 3.044140 | 0.289 | 3.083374 | 19.048 |

2 | 3.044140 | 0.305 | 3.128715 | 19.310 |

3 | 3.100775 | 0.302 | 3.125244 | 18.293 |

4 | 3.149606 | 0.288 | 3.139570 | 18.828 |

Graphics Not Checked (Disabled) |
||||

Method |
1,000,000 tosses |
Time [sec] |
10,000,000 tosses |
Time [sec] |

1 | 3.085158 | 3.872 | 3.085784 | 39.002 |

2 | 3.137368 | 4.600 | 3.140937 | 46.390 |

3 | 3.132474 | 3.440 | 3.131164 | 34.347 |

4 | 3.141878 | 2.031 | 3.140465 | 19.842 |

**Discussion
**The

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:

// update statistics INC (TossCount); IF ABS(D) <= 0.5*SinTheta THEN INC (CrossCount); |

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:

// define midpoint of needle in throwing space xMid := Random(BitmapTosses.Width+1); // determine placement of needle midpoint yMid := PixelsBetweenParallelLines * Random(1+ParallelLineCount) + ROUND(D*PixelsBetweenParallelLines); |

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

xDelta := ROUND(HalfLength*CosTheta); yDelta := ROUND(HalfLength*SinTheta); BitmapTosses.Canvas.Pen.Color := RGB( Random(250), Random(250), Random(250) ); BitmapTosses.Canvas.MoveTo(xMid-xDelta, yMid-yDelta); BitmapTosses.Canvas.LineTo(xMid+xDelta, yMid+yDelta) |

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:

angle := ArcTan2(y,x); SimplePantograph.MapRealToPixel(angle, ABS(D), i,j); IF ABS(D) <= 0.5*SinTheta THEN BitmapGraph.Canvas.Pixels[i,j] := clRed ELSE BitmapGraph.Canvas.Pixels[i,j] := clBlack; |

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:

// UpdateDisplay ImageTosses.Picture.Graphic := BitmapTosses; ImageGraph.Picture.Graphic := BitmapGraph; |

The estimate for p is updated in the *UpdatePIvalue*
method:

PROCEDURE TFormBuffon.UpdatePiValue; VAR PiApprox: DOUBLE; BEGIN LabelTosses.Caption := FormatFloat(',0', TossCount) + Plural(TossCount, ' Toss', ' Tosses') + ', ' + FormatFloat(',0', CrossCount) + Plural(CrossCount, ' Crossing', ''); IF CrossCount > 0 // can't divide by zero THEN BEGIN PiApprox := 2.0*TossCount/CrossCount; LabelPIApproxValue.Caption := Format('%.6f', [PiApprox]); LabelPIApprox.Visible := TRUE END ELSE BEGIN LabelPIApproxValue.Caption := ''; LabelPIApprox.Visible := FALSE END END; |

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.

type TDirectionVectorRoutine = PROCEDURE (VAR x, y, SinTheta, CosTheta: Single); |

The *TFormBuffon* class includes a private *DirectionVectorRoutine*
variable

private DirectionVectorRoutine : TDirectionVectorRoutine; ... |

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

// Setup default random direction generator DirectionVectorRoutine := DirectionMethod2; RadioButtonMethod2.Checked := TRUE; |

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.

procedure TFormBuffon.RadioButtonMethodClick(Sender: TObject); begin CASE (Sender AS TRadioButton).Tag OF 1: DirectionVectorRoutine := DirectionMethod1; 2: DirectionVectorRoutine := DirectionMethod2; 3: DirectionVectorRoutine := DirectionMethod3; 4: DirectionVectorRoutine := DirectionMethod4; ELSE // do nothing (should never happen) END; ResetState end; |

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.

// My original, naive and INCORRECT APPROACH. // Assume PI can't be used anywhere since that's what we're calculating. // Avoid use of SIN and COS. PROCEDURE DirectionMethod1 (VAR x, y, SinTheta, CosTheta: Single); VAR hypotenuse : SINGLE; BEGIN // define direction vector (x,y) to assign needle orientation x := 2.0*(Random - 0.5); // -1.0 <= x <= 1.0 y := Random; // 0.0 <= y <= 1.0 IF x = 0.0 THEN BEGIN SinTheta := 1.0; CosTheta := 0.0 END ELSE BEGIN hypotenuse := SQRT( SQR(x) + SQR(y) ); SinTheta := y / hypotenuse; CosTheta := x / hypotenuse // needed for plotting needle END END {DirectionMethod1}; |

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:

// Gustavo Massaccesi's suggested improvement. Force (x,y)
value to be // inside unit semicircle within rectangle. PROCEDURE DirectionMethod2 (VAR x, y, SinTheta, CosTheta: Single); VAR hypotenuse : SINGLE; BEGIN // Generate random numbers until one is within semicircle. This // ensures all directions have an equal probability. REPEAT // define direction vector (x,y) to assign needle orientation x := 2.0*(Random - 0.5); // -1.0 <= x <= 1.0 y := Random; // 0.0 <= y <= 1.0 hypotenuse := SQRT( SQR(x) + SQR(y) ); UNTIL hypotenuse <= 1.0; IF x = 0.0 THEN BEGIN SinTheta := 1.0; CosTheta := 0.0 END ELSE BEGIN SinTheta := y / hypotenuse; CosTheta := x / hypotenuse // needed for plotting needle END END {DirectionMethod2}; |

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.

// Rainer Elsing's suggested improvement. ASSUME that we
don't know // the value of PI since we're trying to compute it! But allow // evaluation of SIN and COS functions. PROCEDURE DirectionMethod3 (VAR x, y, SinTheta, CosTheta: Single); CONST AnyNumber = 100; // Assumes value is >> Circumference VAR ArcLength: Single; BEGIN ArcLength := Random*AnyNumber; // Force from 3rd or 4th quadrants back to 1st or 2nd SinTheta := ABS( SIN(ArcLength) ); CosTheta := COS(ArcLength); x := CosTheta; y := SinTheta END {DirectionMethod3}; |

**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^{3}
= 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 2^{3}-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:

PROCEDURE DirectionMethod4 (VAR x, y, SinTheta, CosTheta:
Single); VAR index: INTEGER; BEGIN index := Random(LookupArrayCount+1); x := xArray[index]; y := yArray[index]; // 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.

**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,

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