Natural Pi #1 - Sand

Question

Goal

Generate (N) random line segments of uniform length (l), check if they cross the equidistant (t) parallel lines.

Simulation

What are we simulating? Buffon's needle. Smooth out the sand in your sandbox, draw a set of equally spaced parallel lines (call the distance in between t). Take a straight stick of length l and drop it N times into the sandbox. Let the number of times it crossed a line be c. Then Pi = (2 * l * n) / (t * c)!

How are we simulating this?

• Take input N,t,l
• With N, t, l all being positive integers
• Do the following N times:
  • Generate a uniformly random integer coordinate x,y
  • With 1 <= x, y <= 10^6
  • x,y is the center of a line segment of length l
  • Generate a uniformly random integer a
  • With 1 <= a <= 180
  • Let P be the point where the line segment would cross the x-axis
  • Then a is the angle (x,y), P, (inf,0)
• Count the number, c, of line segments that cross the line x = i*t for any integer i
• Return (2 * l * N) / (t * c)

Specification

• Input
  • Flexible, take input in any of the standard ways (eg function parameter,STDIN) and in any standard format (eg String, Binary)
• Output
  • Flexible, give output in any of the standard ways (eg return, print)
  • White space, trailing and leading white space is acceptable
  • Accuracy, please provide at least 4 decimal places of accuracy (ie 3.1416)
• Scoring
  • Shortest code wins!

Test Cases

Your output may not line up with these, because of random chance. But on average, you should get about this much accuracy for the given value of N, t, l.

Input (N,t,l)    ->  Output 
-----------        ------
10,10,5          -> ?.????
10,100,50        -> ?.????
1000,1000,600    -> 3.????
10000,1000,700   -> 3.1???
100000,1000,700  -> 3.14??

TL;DR

These challenges are simulations of algorithms that only require nature and your brain (and maybe some re-usable resources) to approximate Pi. If you really need Pi during the zombie apocalypse, these methods don't waste ammo! There are nine challenges total.

Answer 1

R, 113 100 75 70 68 67 65 59 63 57 bytes

As a statistical, functional programming language, it's not surprising that R is fairly well-suited to this kind of task. The fact that most functions can take vectorized input is really helpful for this problem, as rather than looping over N iterations, we just pass around vectors of size N. Thanks to @Billywob for some suggestions that lead to cutting off 4 bytes. Many thanks to @Primo for patiently explaining to me how my code wasn't working for cases where t > l, which is now fixed.

pryr::f(2*l*N/t/sum(floor(runif(N)+sinpi(runif(N))*l/t)))

Try it online!

Sample output:

N=1000, t=1000, l=500
3.037975

N=10000, t=1000, l=700
3.11943

N=100000, t=1000, l=700
3.140351

Explanation

The problem boils down to determining whether the two x values of the needle are on either side of a parallel line. This has some important consequences:

1. y-values are irrelevant
2. The absolute location on the x-axis is irrelevant, only the position relative to the closest parallel lines.

Essentially this is a task in a 1-dimensional space, where we generate a line with length in [0, l] (the angle a determines this length), and then we check to see how many times this length exceeds t. The rough algorithm is then:

1. Sample x1 values from [0, 1000000]. Since parallel lines occur at every tth point along the x-axis, the relative x-position is x modulo t.
2. Sample an angle a.
3. Calculate the x2 position based on a.
4. Check how many times x1+x2 fits into t, i.e. take the floor of (x1+x2)/t.

Sampling N numbers in [0, 1e6] modulo t is equivalent to simply sampling N numbers in [0, t]. Since (x1+x2)/t is equivalent to x1/t + x2/t, the first step becomes sampling from [0, t]/t, i.e. [0, 1]. Lucky for us, that is the default range for R's runif function, which returns N real numbers from 0 to 1 from a uniform distribution.

                          runif(N)

We repeat this step to generate a, the angle of the needle.

                                         runif(N)

These numbers are is interpreted as a half-turn (i.e. .5 is 90 degrees). (The OP asks for degrees from 1 to 180, but in comments it's clarified that any method is allowed if it is as or more precise.) For an angle θ, sin(θ) gives us the x-axis distance between the ends of the needle. (Normally you'd use the cosine for something like this; but in our case, we are considering the angle θ as being relative to the y-axis, not the x-axis (that is, a value of 0 degrees goes up, not right), and therefore we use the sine, which basically phase-shifts the numbers.) Multiplied by l this gives us the x location of the end of the needle.

                                   sinpi(runif(N))*l

Now we divide by t and add the x1 value. This yields (x1+x2)/t, which is how far the needle protrudes from x1, in terms of number of parallel lines. To get the integer of how many lines were crossed, we take the floor.

                    floor(runif(N)+sinpi(runif(N))*l/t)

We calculate the sum, giving us the count c of how many lines are crossed by needles.

                sum(floor(runif(N)+sinpi(runif(N))*l/t))

The rest of the code is just implementing the formula for approximating pi, that is, (2*l*N)/(t*c). We save some bytes on parentheses by taking advantage of the fact that (2*l*N)/(t*c) == 2*l*N/t/c:

        2*l*N/t/sum(floor(runif(N)+sinpi(runif(N))*l/t))

And the whole thing is wrapped into an anonymous function:

pryr::f(2*l*N/t/sum(floor(runif(N)+sinpi(runif(N))*l/t)))

Answer 2

Perl, 97 bytes

#!perl -p
/ \d+/;$_*=2*$'/$&/map{($x=(1+~~rand 1e6)/$&)-$a..$x+($a=$'/$&/2*sin~~rand(180)*71/4068)-1}1..$_

Counting the shebang as one, input is taken from stdin, space separated. If non-integer random values were allowed, this could be somewhat shorter.

I have taken one liberty, approximating π/180 as 71/4068, which is accurate within 1.48·10^-9.

Sample Usage

$ echo 1000000 1000 70000 | perl pi-sand.pl
3.14115345174061

---

More-or-less mathematically equivalent substitutions

Assuming the x-coordinate represents the left-most point of the needle, rather than its middle, as specified in the problem description:

89 bytes

#!perl -p
/ \d+/;$_*=2*$'/$&/map{($x=(1+~~rand 1e6)/$&)..$x+($'/$&*sin~~rand(180)*71/4068)-1}1..$_

The problem specifies that x is to be sampled as a random integer. If we project the line spacing to a gap of one, this will leave us with values of the form n/t with 0 <= n < t, not necessarily uniform, if t does not evenly divide 1e6. Assuming that a uniform distribution is nonetheless acceptable:

76 bytes

#!perl -p
/ \d+/;$_*=2*$'/$&/map{($x=rand)..$x+($'/$&*sin~~rand(180)*71/4068)-1}1..$_

Note that since rand will always be less than one (and therefore truncated to zero), it is not necessary at the start of the range:

70 bytes

#!perl -p
/ \d+/;$_*=2*$'/$&/map{1..(rand)+($'/$&*sin~~rand(180)*71/4068)}1..$_

Assuming that the angle of the needle need not be an integer degree, but only uniformly random:

59 bytes

#!perl -p
/ \d+/;$_*=2*$'/$&/map{1..(rand)+abs$'/$&*sin rand$`}1..$_

Assuming that the angle may be any uniform distribution:

52 bytes

#!perl -p
/ \d+/;$_*=2*$'/$&/map{1..(rand)+abs$'/$&*sin}1..$_

The above is a mathematically correct simulation of Buffon's Needle. However, at this point I think most people would agree that this isn't actually what the question asked for.

---

Really pushin' it

We could just throw away half of the test cases, whenever the second endpoint is to the left of the first (rather than swapping them):

47 bytes

#!perl -p
/ \d+/;$_*=$'/$&/map{1..(rand)+$'/$&*sin}1..$_

Note that the values of t and l are inconsequential to the results of the experiment. We could just ignore them (implicitly assuming them to be equal):

28 bytes

#!perl -p
$_/=map{1..(rand)+sin}1..$_

Obviously non-competing, but you have to admit, it does have a certain elegance to it.

Answer 3

Python 2, 141 bytes

shameless port of rtumbull, already skipping y because totally not needed.

from math import*
from random import*
lambda N,t,l:(2.*l*N)/(t*sum(randint(1,1e6)%t+abs(cos(randint(1,180)*pi/180))*l>t for _ in range(N)))

Problem is only, that pi is already known in the program.

Here it is (golfable) with unknown pi and no trigonometric functions

def g(N,t,l):
 c=0
 for _ in range(N):
    x,y=gauss(0,1),gauss(0,1);c+=randint(1,1e6)%t+abs(x/sqrt(x*x+y*y))*l>t
 return(2.*l*N)/(t*c)

x,y in g is only for the direction.