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.

About these ads

Pages: 1 2

11 Responses to “Buffon’s Needle”

  1. [...] today’s Programming Praxis exercise, our goal is to approximate pi by using Georges-Louis Leclerc’s [...]

  2. 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)
    
  3. jpverkamp said

    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.

  4. spainmc said

    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/

  5. spainmc said

    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;
    }

  6. spainmc said

    Ugh. . . nothing’s posting right. One last try: JS Fiddle

  7. Mike said

    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?

  8. a1911 said
    
    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
    
    
  9. jpverkamp said

    @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)/j
    
  10. Mirko said

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 630 other followers

%d bloggers like this: