Buffon’s Needle
March 15, 2013
To keep things simple, assume that the length of the needle is 1 (we don’t care about the units) which is also the distance between two lines. Trigonometry tells us that the needle will intersect a line if y < sin(θ) / 2, where y is the distance from the center of the needle to the nearest line and θ is the acute angle (from 0 to 90 degrees) the needle makes with the line. The simulation repeatedly “drops” a needle by choosing a random y on the range 0 <= y < 1 and a θ on the range 0 <= θ < 90, with both y and θ real numbers (not integers). Since our sine function takes its arguments in radians, not degrees, we change θ to the range 0 <= θ < 1.571.
(define (buffon n)
(let loop ((drops 0) (hits 0))
(if (= drops n) (/ drops hits 1.0)
(if (< (rand) (/ (sin (* 1.571 (rand))) 2.0))
(loop (+ drops 1) (+ hits 1))
(loop (+ drops 1) hits)))))
Our simulation isn’t really all that close:
> (buffon 1000000)
3.130899789290444
We used the random number generator from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/SOYXlD1f.
Pages: 1 2
[…] 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)))))