Natural Pi #0 - Rock

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! There are eight more challenges to come. Checkout the sandbox post 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, so a natural way to calculate Pi would be to scoop up two buckets (or handfuls) of rocks; count them; see if their gcd 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 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)

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 1

APL, 23 bytes

{.5*⍨6×⍵÷1+.=∨/?⍵2⍴1e6}

Explanation:

• ?⍵2⍴1e6: generate a 2-by-⍵ matrix of random numbers in the range [1..10⁶]
• 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 2

Jelly, 20 18 16 bytes

-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

Answer 3

Python 2, 143 140 132 124 122 124 122 bytes

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

^{thanks to Jonathan Allan for the two-byte save :)}

Answer 4

Mathematica, 49 48 51 bytes

Saved one byte and fixed one bug thanks to @LegionMammal978.

(6#/Count[GCD@@{1,1*^6}~RandomInteger~{2,#},1])^.5&

Answer 5

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 6

Actually, 19 bytes

`6╤;Ju@Ju┤`nkΣß6*/√

Try it online!

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 7

MATL, 22 bytes

1e6Hi3$YrZ}Zd1=Ym6w/X^

Try it online!

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 8

Pyth, 21 bytes

@*6cQ/iMcmhO^T6yQ2lN2

Try it online.

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 9

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 10

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 11

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 12

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 13

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 14

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

Answer 15

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 -minimum 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. 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 16

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 and added test harness

Answer 17

Pyt, 24 bytes

Đ0⇹`⇹ɽɽǤ⁻?ŕ:ŕ⁺;⇹⁻łŕᵮ/6*√

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 18

Perl 6, 56 53 bytes

{sqrt 6*$_/(1..10⁶).roll(*).map(*gcd*==1)[^$_].sum}

Try it online!

Answer 19

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 20

AWK, 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!

I'm surprised that it runs in a reasonable amount of time for 1000000.

Answer 21

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 22

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 23

Japt, 23 18 bytes

*6/UÆ2ÆÒL³öÃrjÃx)¬

Try it