# Natural Pi #0 - Rock

## 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??

• Does our answer need to work for N = 1000000 or is it ok if the program returns e.g. a stack overflow if N is too big? Oct 3, 2016 at 13:41
• @Fatalize if it is a limitation of the language, sure. Otherwise, you need to handle N=10^6. Oct 3, 2016 at 13:43
• Related Oct 3, 2016 at 14:01
• The goal is misleading, it states only one pair of integers is checked. Oct 4, 2016 at 2:08
• Does the upper limit to the random numbers generated need to be exactly 1000000? Would a larger upper limit be acceptable?
– Sok
Oct 4, 2016 at 12:24

# 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

# 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

• ḤRµȷ6Xµ€g2/ċ1÷³6÷½ saves 2 bytes. (ȷ6 is 10^6 in a single nilad, ċ1 counts ones) Oct 3, 2016 at 15:44
• Ah I couldn't work out how to chain it up like that (I tried a few things), and forgot the count 1 trick - thanks (I think ȷ² is a tiny tiny bit faster than ȷ6) Oct 3, 2016 at 15:56
• Might be. Now that I think of it, ȷ² being two links doesn't hurt here, but would require an extra link or ¤ for some use cases Oct 3, 2016 at 16:00
• Ḥȷ6xX€ should work for the random sampling. Oct 3, 2016 at 16:08

# Python 2, 143140132124122124122 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 :)

• According to OP, 1 <= i , j <= 10^6, so you need to use randrange(1,1e6+1). Oct 3, 2016 at 21:38
• Also, it's really strange to have the repl.it link within the language name. A link in the lang name should be to the language's home page, if anything. Put your repl.it link as a separate link below your code. Oct 3, 2016 at 21:42
• @mbomb007 Good point, I've fixed it :) Been a while!
Oct 4, 2016 at 13:00
• k=lambda:r.randrange(1e6)+1 saves two bytes Oct 5, 2016 at 19:05
• @JonathanAllan good catch, thanks!
Oct 6, 2016 at 13:01

# Mathematica, 4948 51 bytes

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

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

• You can save a byte: (6#/Count[GCD@@1*^6~RandomInteger~{2,#},1])^.5& Oct 15, 2016 at 15:42
• Also, 1*^6 should be replaced with {1,1*^6} to ensure that i, j ≠ 0. Oct 15, 2016 at 15:47

# R, 10399959998 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

• You can simply use sample(10^6,N). Not only is it shorter, it is also much more efficient.
– asac
Oct 3, 2016 at 17:26
• I may be wrong, but shouldn't the sample be used with replace=T for a properly uniform random integers. For example sample(10,10) is guaranteed to return all the numbers in 1:10, whereas sample(10,10,T) will produce a random selection where numbers can be repeated. Oct 3, 2016 at 19:02
• @MickyT You're absolutely correct, I just realized this a few minutes ago myself. I'm not entirely sure how this plays out mathematically in this instance - as far as I can tell, both methods are roughly equally accurate. I'll edit my post to add this information. Oct 3, 2016 at 19:21
• Both methods are equally accurate when N<<10^6. To handle arbitrarily big N, you have to sample with replacement, good catch.
– asac
Oct 4, 2016 at 8:30

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

• i, j aren't allowed to be 0 Oct 3, 2016 at 17:00
• @isaacg They aren't. If you'll read the explanation, it says that the random values are selected from [0, 10**6), then incremented.
– user45941
Oct 3, 2016 at 17:22

# 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


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


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

• I <3 Scala! Especially, because it sometimes really needs an explanation. Oct 3, 2016 at 21:31
• @RomanGräf To be honest, the only things I think might be unclear are 6f, Seq.fill and math.random. Oct 3, 2016 at 21:35

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


# JavaScript (ES7), 10795 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));
}

• In the Ungolfed Version gcd calls the the function g Oct 4, 2016 at 5:54
• r=_=> is that code or a drawing? Oct 4, 2016 at 14:36
• n=>(n*6/(r=_=>Math.random()*1e6,g=(a,b)=>b?g(b,a%b):a>-2,q=n=>n&&g(~r(),~r())+q(n-1))(n))**.5 1B shorter
– l4m2
Dec 26, 2017 at 10:00
• n=>(n*6/(q=_=>n--&&q(r=_=>Math.random()*1e6)+g(~r(),~r()))(g=(a,b)=>b?g(b,a%b):a>-2))**.5
– l4m2
Dec 26, 2017 at 10:05

# PHP, 8277 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.

• Saved 3 bytes by using $argn # 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  # R, 94 bytes N=scan();a=replicate(N,{x=sample(1e6,2);q=1:x;max(q[!x%%q&!x%%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). • If you're not worried about warnings, 1:x will work. Oct 4, 2016 at 20:30 ## 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


# 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
• You don't need parentheses around the first n, t+=1 can become t++, you can condense your int declarations into one line, i.e. int c=n,t=0,x,y;, and !=0 (I think) can become >0. That should save 12 bytes overall. That's a neat way of finding the GCD of x and y, though.
– TNT
Oct 5, 2016 at 17:25

# Perl 6, 56 53 bytes

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


Try it online!

# Frink, 84 89

r[]:=random[10^6]+1
g=n=eval[input]
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,  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 ½

• That's 90 bytes and 88 chars...? Jul 9, 2017 at 21:50
• Thanks for catching that. I didn't count newlines, and while ², ³ are only 1 byte ⁶ is more. I fixed it to 89 bytes with no final newline. Jul 11, 2017 at 0:59
• You have not fixed the verbose code. Jul 11, 2017 at 3:35
• It's not a one-to-one match anyway with spacing, the quotes and numbers, etc. Jul 11, 2017 at 7:16

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

# 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


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

# Japt, 23 18 bytes

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


Try it