[Question]
  [
## Goal
Create a program/function that takes an input `N`, check if `N` random pairs of integers are relatively prime, and returns `sqrt(6 * N / #coprime)`.
## 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](https://medium.com/the-physics-arxiv-blog/how-mathematicians-used-a-pump-action-shotgun-to-estimate-pi-c1eb776193ef#.eld2u7w76)! There are eight more challenges to come. Checkout the [sandbox post](https://codegolf.meta.stackexchange.com/a/10035/43214) to make recommendations.
## Simulation
**What are we simulating?** Well, the probability that two random integers are relatively prime (ie coprime or gcd==1) is [`6/Pi/Pi`](https://mathoverflow.net/a/97045), so a natural way to calculate Pi would be to scoop up two buckets (or handfuls) of rocks; count them; see if their [gcd](https://en.wikipedia.org/wiki/Greatest_common_divisor) is 1; repeat. After doing this a ~~couple~~ lot of times, `sqrt(6.0 * total / num_coprimes)` will tend towards `Pi`. If calculating the square-root in post-apocalyptic world makes you nervous, don't worry! There is [Newton's Method](https://en.wikipedia.org/wiki/Newton%27s_method#Square_root_of_a_number) for that. 
**How are we simulating this?** 
* Take input `N`
* Do the following `N` times:
	+ Uniformly generate random positive integers, `i` and `j`
	+ With `1 <= i , j <= 10^6`
	+ If `gcd(i , j) == 1`: `result = 1`
	+ Else: `result = 0`
* Take the sum of the `N` results, `S`
* Return `sqrt(6 * N / S)`
[![enter image description here](https://i.stack.imgur.com/xjU5L.jpg)](https://i.stack.imgur.com/xjU5L.jpg)
## 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`.
```
Input     ->  Output 
-----         ------
100       ->  3.????
10000     ->  3.1???
1000000   ->  3.14??
```
  
[Answer]
# APL, 23 bytes
```
{.5*⍨6×⍵÷1+.=∨/?⍵2⍴1e6}
```
Explanation:
* `?⍵2⍴1e6`: generate a 2-by-⍵ matrix of random numbers in the range [1..106]
* `1+.=∨/`: get the GCD of each pair and see how many are equal to 1. This calculates S.
* `.5*⍨6×⍵÷`: (6 × ⍵ ÷ S)0.5
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~20 18~~ 16 [bytes](https://github.com/DennisMitchell/jelly/wiki/Code-page)
-2 bytes thanks to @Pietu1998 (chain & use count 1s, `ċ1` in place of less than two summed `<2S`)
-2 bytes thanks to @Dennis (repeat 1e6 multiple times before sampling to avoid chaining) 
```
Ḥȷ6xX€g2/ċ1÷³6÷½
```
(Extremely slow due to the random function)
How?
```
Ḥȷ6xX€g2/ċ1÷³6÷½ - Main link: n
 ȷ6              - 1e6
   x             - repeat
Ḥ                -     double, 2n
    X€           - random integer in [1,1e6] for each
       2/        - pairwise reduce with
      g          -     gcd
         ċ1      - count 1s
           ÷     - divide
            ³    - first input, n
             6   - literal 6
              ÷  - divide
               ½ - square root
```
**[TryItOnline](http://jelly.tryitonline.net/)**
[Answer]
# Python 2, ~~143~~ ~~140~~ ~~132~~ ~~124~~ ~~122~~ ~~124~~ [122 bytes](https://goo.gl/7RS871)
It has been quite some time since I have tried golfing, so I may have missed something here! Will be updating as I shorten this.
```
import random as r,fractions as f
n,s=input(),0
k=lambda:r.randrange(1e6)+1
exec's+=f.gcd(k(),k())<2;'*n
print(6.*n/s)**.5
```
Test me [here!](https://repl.it/DnAQ/2)
*thanks to Jonathan Allan for the two-byte save :)*
[Answer]
# Mathematica, ~~49~~ ~~48~~ 51 bytes
Saved one byte and fixed one bug thanks to [@LegionMammal978](https://codegolf.stackexchange.com/users/33208/legionmammal978).
```
(6#/Count[GCD@@{1,1*^6}~RandomInteger~{2,#},1])^.5&
```
[Answer]
# R, ~~103~~ ~~99~~ ~~95~~ ~~99~~ ~~98~~ 94 bytes
Can likely be golfed down a bit. Cut down 4 bytes due to @antoine-sac, and another 4 bytes by defining an alias for `sample`, using `^.5` instead of `sqrt`, and `1e6` instead of `10^6`. Added 4 bytes to ensure that the sampling of `i` and `j` is truly uniform. Removed one byte after I realized that `6*N/sum(x)` is the same as `6/mean(x)`. Used `pryr::f` instead of `function(x,y)` to save 4 bytes.
```
N=scan()
s=sample
g=pryr::f(ifelse(o<-x%%y,g(y,o),y))
(6/mean(g(s(1e6,N,1),s(1e6,N,1))==1))^.5
```
Sample output:
```
N=100     -> 3.333333
N=10000   -> 3.137794
N=1000000 -> 3.141709
```
[Answer]
## Actually, 19 bytes
```
`6╤;Ju@Ju┤`nkΣß6*/√
```
[Try it online!](http://actually.tryitonline.net/#code=YDbilaQ7SnVASnXilKRgbmvOo8OfNiov4oia&input=MTAw)
Explanation:
```
`6╤;Ju@Ju┤`nkΣß6*/√
`6╤;Ju@Ju┤`n         do this N times:
 6╤;                   two copies of 10**6
    Ju                 random integer in [0, 10**6), increment
      @Ju              another random integer in [0, 10**6), increment
         ┤             1 if coprime else 0
            kΣ       sum the results
              ß      first input again
               6*    multiply by 6
                 /   divide by sum
                  √  square root
```
[Answer]
# [MATL](http://github.com/lmendo/MATL), 22 bytes
```
1e6Hi3$YrZ}Zd1=Ym6w/X^
```
[Try it online!](http://matl.tryitonline.net/#code=MWU2SGkzJFlyWn1aZDE9WW02dy9YXg&input=MTAwMDA)
```
1e6      % Push 1e6
H        % Push 2
i        % Push input, N
3$Yr     % 2×N matrix of uniformly random integer values between 1 and 1e6
Z}       % Split into its two rows. Gives two 1×N arrays
Zd       % GCD, element-wise. Gives a 1×N array
1=       % Compare each entry with 1. Sets 1 to 0, and other values to 0
Ym       % Mean of the array
6w/      % 6 divided by that
X^       % Square root. Implicitly display
```
[Answer]
# Pyth, 21 bytes
```
@*6cQ/iMcmhO^T6yQ2lN2
```
[Try it online.](http://pyth.herokuapp.com/?code=%40%2a6cQ%2FiMcmhO%5ET6yQ2lN2&input=10000&debug=0)
### Explanation
```
                Q          input number
               y           twice that
         m                 map numbers 0 to n-1:
             T                 10
            ^ 6                to the 6th power
           O                   random number from 0 to n-1
          h                    add one
        c        2         split into pairs
      iM                   gcd of each pair
     /            lN       count ones
   cQ                      divide input number by the result
 *6                        multiply by 6
@                   2      square root
```
[Answer]
# Scala, ~~149~~ 126 bytes
```
val& =BigInt
def f(n: Int)={math.sqrt(6f*n/Seq.fill(n){val i,j=(math.random*99999+1).toInt
if(&(i).gcd(&(j))>1)0 else 1}.sum)}
```
Explanation:
```
val& =BigInt                //define & as an alias to the object BigInt, because it has a gcd method
def f(n:Int)={              //define a method
  math.sqrt(                //take the sqrt of...
    6f * n /                //6 * n (6f is a floating-point literal to prevent integer division)
    Seq.fill(n){            //Build a sequence with n elements, where each element is..
      val i,j=(math.random*99999+1).toInt //take 2 random integers
      if(&(i).gcd(&(j))>1)0 else 1        //put 0 or 1 in the list by calling
                                          //the apply method of & to convert the numbers to
                                          //BigInt and calling its bcd method
    }.sum                   //calculate the sum
  )
}
```
[Answer]
## Racket 92 bytes
```
(λ(N)(sqrt(/(* 6 N)(for/sum((c N))(if(= 1(gcd(random 1 1000000)(random 1 1000000)))1 0)))))
```
Ungolfed: 
```
(define f
  (λ (N)
    (sqrt(/ (* 6 N) 
            (for/sum ((c N))
              (if (= 1
                     (gcd (random 1 1000000)
                          (random 1 1000000)))
                  1 0)
              )))))
```
Testing: 
```
(f 100)
(f 1000)
(f 100000)
```
Output: 
```
2.970442628930023
3.188964020716403
3.144483068444827
```
[Answer]
# JavaScript (ES7), ~~107~~ ~~95~~ 94 bytes
```
n=>(n*6/(r=_=>Math.random()*1e6+1|0,g=(a,b)=>b?g(b,a%b):a<2,q=n=>n&&g(r(),r())+q(n-1))(n))**.5
```
The ES6 version is exactly 99 bytes, but the ES7 exponentiation operator `**` saves 5 bytes over `Math.sqrt`.
### Ungolfed
```
function pi(n) {
  function random() {
    return Math.floor(Math.random() * 1e6) + 1;
  }
  function gcd(a, b) {
    if (b == 0)
      return a;
    return gcd(b, a % b);
  }
  function q(n) {
    if (n == 0)
      return 0;
    return (gcd(random(), random()) == 1 ? 1 : 0) + q(n - 1));
  }
  return Math.sqrt(n * 6 / q(n));
}
```
[Answer]
# PHP, ~~82~~ ~~77~~ 74 bytes
```
for(;$i++<$argn;)$s+=2>gmp_gcd(rand(1,1e6),rand(1,1e6));echo(6*$i/$s)**.5;
```
Run like this:
```
echo 10000 | php -R 'for(;$i++<$argn;)$s+=2>gmp_gcd(rand(1,1e6),rand(1,1e6));echo(6*$i/$s)**.5;' 2>/dev/null;echo
```
# Explanation
Does what it says on the tin. Requires PHP\_GMP for `gcd`.
# Tweaks
* Saved 3 bytes by using `$argn`
[Answer]
# Perl, 64 bytes
```
sub r{1+~~rand 9x6}$_=sqrt$_*6/grep{2>gcd r,r}1..$_
```
Requires the command line option `-pMntheory=gcd`, counted as 13. Input is taken from stdin.
**Sample Usage**
```
$ echo 1000 | perl -pMntheory=gcd pi-rock.pl
3.14140431218772
```
[Answer]
# R, 94 bytes
```
N=scan();a=replicate(N,{x=sample(1e6,2);q=1:x[1];max(q[!x[1]%%q&!x[2]%%q])<2});(6*N/sum(a))^.5
```
Relatively slow but still works. Replicate N times a function that takes 2 random numbers (from 1 to 1e6) and checks if their gcd is less than 2 (using an [old gcd function of mine](https://codegolf.stackexchange.com/a/48845/6741)).
[Answer]
## PowerShell v2+, ~~118~~ 114 bytes
```
param($n)for(;$k-le$n;$k++){$i,$j=0,1|%{Random -mi 1};while($j){$i,$j=$j,($i%$j)}$o+=!($i-1)}[math]::Sqrt(6*$n/$o)
```
Takes input `$n`, starts a `for` loop until `$k` equals `$n` (implicit `$k=0` upon first entering the loop). Each iteration, get new `Random` numbers `$i` and `$j` (the `-mi`nimum `1` flag ensure we're `>=1` and no maximum flag allows up to `[int]::MaxValue`, which is allowed by the OP since it's larger than `10e6`).
We then go into a [GCD `while` loop](https://codegolf.stackexchange.com/a/83562/42963). Then, so long as the GCD is `1`, `$o` gets incremented. At the end of the `for` loop, we do a simple `[math]::Sqrt()` call, which gets left on the pipeline and output is implicit.
Takes about 15 minutes to run with input `10000` on my ~1 year old Core i5 laptop.
### Examples
```
PS C:\Tools\Scripts\golfing> .\natural-pi-0-rock.ps1 100
3.11085508419128
PS C:\Tools\Scripts\golfing> .\natural-pi-0-rock.ps1 1000
3.17820863081864
PS C:\Tools\Scripts\golfing> .\natural-pi-0-rock.ps1 10000
3.16756133579975
```
[Answer]
# Java 8, ~~164~~ 151 bytes
```
n->{int c=n,t=0,x,y;while(c-->0){x=1+(int)(Math.random()*10e6);y=1+(int)(Math.random()*10e6);while(y>0)y=x%(x=y);if(x<2)t++;}return Math.sqrt(6f*n/t);}
```
## Explanation
```
n->{
    int c=n,t=0,x,y;
    while(c-->0){                          // Repeat n times
        x=1+(int)(Math.random()*10e6);     // Random x
        y=1+(int)(Math.random()*10e6);     // Random y
        while(y>0)y=x%(x=y);               // GCD
        if(x<2)t++;                        // Coprime?
    }
    return Math.sqrt(6f*n/t);              // Pi
}
```
### Test Harness
```
class Main {
    public static interface F{ double f(int n); }
    public static void g(F s){
        System.out.println(s.f(100));
        System.out.println(s.f(1000));
        System.out.println(s.f(10000));
    }
    public static void main(String[] args) {
        g(
            n->{int c=n,t=0,y,x;while(c-->0){x=1+(int)(Math.random()*10e6);y=1+(int)(Math.random()*10e6);while(y>0)y=x%(x=y);if(x<2)t++;}return Math.sqrt(6f*n/t);}
        );
    }
}
```
**Update**
* **-13** [16-10-05] Thanks to [*@TNT*](https://codegolf.stackexchange.com/users/32515/tnt) and added test harness
[Answer]
# [Pyt](https://github.com/mudkip201/pyt), 24 bytes
```
Đ0⇹`⇹ɽɽǤ⁻?ŕ:ŕ⁺;⇹⁻łŕᵮ/6*√
```
[Try it online!](https://tio.run/##K6gs@f//yASDR@07E4D45N6Te48vedS42/7oVKujUx817rIGigL5R5uOTn24dZ2@mdajjln//xsaAAEA "Pyt – Try It Online")
```
Đ                           implicit input; Đuplicate on stack (N, i) where i is a counter
 0                          push 0 (S)
  ⇹                         swap top two items on stack
   `             ł          do... while top of stack is not 0
    ⇹                       swap top two items on stack
     ɽɽ                     get two ɽandom 32-bit integers
       Ǥ⁻?                  if the ǤCD of the two random integers is >1:
          ŕ                 ŕemove the GCD
           :ŕ⁺              otherwise ŕemove the GCD and increment S
              ;⇹⁻           either way, swap and decrement N
                 ŕ          ŕemove i (which is now 0)
                  ᵮ         cast S to a ᵮloat
                   /        N/S
                    6*√     sqrt(6*N/S); implicit print
```
[Answer]
# [Perl 6](https://perl6.org), ~~56~~ 53 bytes
```
{sqrt 6*$_/(1..10⁶).roll(*).map(*gcd*==1)[^$_].sum}
```
[Try it online!](https://tio.run/nexus/perl6#Ky1OVSgz00u25sqtVFBLU7D9X11cWFSiYKalEq@vYainZ2jwqHGbpl5Rfk6OhpamXm5igYZWenKKlq2toWZ0nEp8rF5xaW7t/@LESoU0BUMDay4Yy8DA@j8A "Perl 6 – TIO Nexus")
[Answer]
# Frink,  ~~84~~ 89
```
r[]:=random[10^6]+1
g=n=eval[input[1]]
for a=1to n
g=g-1%gcd[r[],r[]]
println[(6*n/g)^.5]
```
I got lucky: **g=n=...** saves a byte over **g=0 n=...**; and **1%gcd()** gives (0,1) vs (1,0) so I can subtract. 
And unlucky: **n** is preassigned and **a** used because loop variables and their bounds are local and undefined outside the loop.
## Verbose
```
r[] := random[10^6] + 1     // function. Frink parses Unicode superscript!
g = n = eval[input[""]]     // input number, [1] works too
for a = 1 to n              // repeat n times
   g = g - 1%gcd[r[], r[]]  // subtract 1 if gcd(i, j) > 1
println[(6*n/g)^.5]         // ^.5 is shorter than sqrt[x], but no super ".", no ½
```
[Answer]
# [AWK](https://www.gnu.org/software/gawk/manual/gawk.html), 109 bytes
```
func G(p,q){return(q?G(q,p%q):p)}{for(;i++<$0;)x+=G(int(1e6*rand()+1),int(1e6*rand()+1))==1;$0=sqrt(6*$0/x)}1
```
[Try it online!](https://tio.run/##ZcixCoMwEADQX@mQwp0JmFscjKFjvkOsghRi7owoiN@egqtvfP3@K2Xa4vAKkAzjKWPeJAJ/ArBJb8Y24XVOi4Cbte6UdXhoH2COGWhsKunjF1ATmseg9@SU9StLhqZStj7wolLI3v4 "AWK – Try It Online")
I'm surprised that it runs in a reasonable amount of time for 1000000.
[Answer]
# [Pyt](https://github.com/mudkip201/pyt), ~~37~~ 35 bytes
```
←Đ0⇹`25*⁶⁺Đ1⇹ɾ⇹1⇹ɾǤ1=⇹3Ș+⇹⁻łŕ⇹6*⇹/√
```
Explanation:
```
←Đ                                              Push input onto stack twice
  0                                             Push 0
   ⇹                                            Swap top two elements of stack
    `                      ł                    Repeat until top of stack is 0
     25*⁶⁺Đ1⇹ɾ⇹1⇹ɾ                              Randomly generate two integers in the range [1,10^6]
                  Ǥ1=                           Is their GCD 1?
                     ⇹3Ș                        Reposition top three elements of stack
                        +                       Add the top 2 on the stack
                         ⇹⁻                     Swap the top two and subtract one from the new top of the stack
                            ŕ                   Remove the counter from the stack
                             ⇹                  Swap the top two on the stack
                              6*                Multiply top by 6
                                ⇹               Swap top two
                                 /              Divide the second on the stack by the first
                                  √             Get the square root
```
[Answer]
# J, 27 Bytes
```
3 :'%:6*y%+/(1:=?+.?)y#1e6'
```
### Explanation:
```
3 :'                      '  | Explicit verb definition
                     y#1e6   | List of y copies of 1e6 = 1000000
            (1:=?+.?)        | for each item, generate i and j, and test whether their gcd is 1
          +/                 | Sum the resulting list
      6*y%                   | Divide y by it and multiply by six
    %:                       | Square root
```
Got pretty lucky with a `3.14157` for `N = 10000000`, which took `2.44` seconds.
[Answer]
# [Japt](https://github.com/ETHproductions/japt), ~~23~~ 18 [bytes](https://en.wikipedia.org/wiki/ISO/IEC_8859-1)
```
*6/UÆ2ÆÒL³öÃrjÃx)¬
```
[Try it](https://petershaggynoble.github.io/Japt-Interpreter/?v=1.4.6&code=KjYvVcYyxtJMs/bDcmrDeCms&input=MTAwMDA)
]