Buffon’s Needle
March 15, 2013
Yesterday was pi day, so today we will estimate pi by a method that dates to 1777: if you randomly drop needles onto a flat surface with lines separated by the length of a needle, then the number of needles dropped, divided by the number of needles that intersect a line, will equal pi. You can determine if a needle intersects a line with a little bit of trigonometry. The method is called Buffon’s needle because it was discovered by the French naturalist Georges-Louis Leclerc, the Comte de Buffon. We did a similar exercise to this one a few years ago.
Your task is to write a program that calculates an approximate value of pi by simulating the dropping of a million needles. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
[…] today’s Programming Praxis exercise, our goal is to approximate pi by using Georges-Louis Leclerc’s […]
My Haskell solution (see http://bonsaicode.wordpress.com/2013/03/15/programming-praxis-buffons-needle/ for a version with comments):
import Control.Applicative import System.Random rnds :: IO [Double] rnds = fmap randoms newStdGen buffon :: Int -> IO Double buffon n = (fromIntegral n /) . sum . take n <$> (zipWith (\y t -> if y < sin (t*pi/2) / 2 then 1 else 0) <$> rnds <*> rnds)Had a bit of fun with this one. My solution is in JavaScript, with neat demonstration / visualization that can drop the needles at varying speeds.
Linky: Approximating Pi with Buffon’s Needle
// Drop a bunch of needles var x1, y1, theta, x2, y2; for (var i = 0; i < numberToDrop; i++) { // Drop a new needle x1 = (Math.random() * (width - 2 * needleLength)) + needleLength; y1 = (Math.random() * (height - 2 * needleLength)) + needleLength; theta = Math.random() * 2 * Math.PI; x2 = x1 + needleLength * Math.cos(theta); y2 = y1 + needleLength * Math.sin(theta); // Check if it crosses a line // Yes, this isn't the best way to do this var crossesLine = false; for (var x = needleLength / 2; x <= canvas.width; x += needleLength) { if ((x1 <= x && x <= x2) || (x2 <= x && x <= x1)) crossesLine = true; } // Record the toss tossed += 1; if (crossesLine) crossed += 1; } console.log("pi ~ " + (2 * tossed / crossed));If you run at the maximum speed for a few minutes, I’ve seen the approximation get out to four decimal places at least. Granted, that was after about a hundred billion needles, but still.
Also, side note: the way I worked it out pi equals two times the number dropped divided by the number crossing lines. That was confusing for a bit.
Also in JavaScript, also worked out to pi/2. Not sure why though? Any insight?
var crossed = 0;
var pi;
for (var i = 0; i = n || y2 <= 0) {
cross = 1;
}
return cross;
}
http://jsfiddle.net/zDGdB/3/
Didn’t post right. . .
if (confirm(“Run Program”)) {
var tosses = 1000000;
var n = 1; // pin length
Run();
}
function Run() {
var crossed = 0;
var pi;
for (var i = 0; i = n || y2 <= 0) {
cross = 1;
}
return cross;
}
Ugh. . . nothing’s posting right. One last try: JS Fiddle
If y is the distance from the center of the needle to the nearest line, and the lines are 1 unit apart. Then doesn’t y have the range 0 <= y <= 0.5?
import random import math j = 0 k = 100000 for i in range(k): if round(random.random(),4) <= 0.5*math.sin(round((random.random()*1000)%180)*math.pi/180): j+=1 print float(k)/j@spainmc / @Mike: The factor of two is actually in everyone’s code now that I look at it a bit more, it’s just a matter of where. The other three solutions all use a half circle centered on the midpoint and ours both use a full circle centered on an endpoint. You half to correct either of them by a factor of two somewhere, it’s just ours was at the end.
Also, check out the link for HOWTO: Posting Source Code up at the top of the page. :)
@a1911: Why are you converting from radians to degrees and then back to radians? For that matter, why are you rounding everything? This code does pretty much exactly the same thing:
import random import math j = 0 k = 100000 for i in range(k): if random.random() <= 0.5 * math.sin(random.random() * math.pi): j+=1 print float(k)/jNy entry in Common Lisp:
A comment regarding this method to estimate pi: This method seems to employ a circular argument, since we need to know pi in order to estimate it (pi is used in the calculation of the sine). This is unlike the other geometric method, of generating random points in a square, and counting the number of points that fall inside a circle. For that one, we do not need to know pi in advance.
It is interesting that the other method can be expanded to higher dimensions (like generating random points in a cube and counting points inside a sphere). But I could not figure out an equivalent method in one dimension.
Oops, my common lisp code did not make it. Here is a direct paste:
;; Common Lisp
(defconstant +pi/2+ (float (/ pi 2) 1.0))
(defun buffon-s-needle (n)
“Estimate PI via Buffon’s algorithm
http://en.wikipedia.org/wiki/Buffon%27s_needle
”
(loop for i below n
with sum = 0
when (< (random 1.0) (* 0.5
(sin (random +pi/2+))))
do (incf sum)
finally (return-from estimate-pi (/ (float n 1.0)
(float sum 1.0)))))