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.

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: