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)

enter image description here

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

• Fatalize: 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?
• NonlinearFruit: @Fatalize if it is a limitation of the language, sure. Otherwise, you need to handle N=10^6.
• Luis Mendo: Related
• Stack Exchange Broke The Law: The goal is misleading, it states only one pair of integers is checked.
• Sok: Does the upper limit to the random numbers generated need to be exactly 1000000? Would a larger upper limit be acceptable?

Answer 1

APL, 23 bytes

×ばつ⍵÷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.
• ×ばつ⍵÷: (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÷36÷1⁄2

(Extremely slow due to the random function)

How?

Ḥȷ6xX€g2/ċ1÷36÷1⁄2 - 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
 3 - first input, n
 6 - literal 6
 ÷ - divide
 1⁄2 - square root

TryItOnline

• PurkkaKoodari: ḤRµȷ6Xµ€g2/ċ1÷³6÷½ saves 2 bytes. (ȷ6 is 10^6 in a single nilad, ċ1 counts ones)
• Jonathan Allan: 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)
• PurkkaKoodari: 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
• Dennis: Ḥȷ6xX€ should work for the random sampling.

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

• mbomb007: According to OP, 1 <= i , j <= 10^6, so you need to use randrange(1,1e6+1).
• mbomb007: 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.
• Kade: @mbomb007 Good point, I've fixed it :) Been a while!
• Jonathan Allan: k=lambda:r.randrange(1e6)+1 saves two bytes

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&

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

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

• asac: You can simply use sample(10^6,N). Not only is it shorter, it is also much more efficient.
• MickyT: 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.
• rturnbull: @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.
• asac: Both methods are equally accurate when N<<10^6. To handle arbitrarily big N, you have to sample with replacement, good catch.

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

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

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 % ×ばつN matrix of uniformly random integer values between 1 and 1e6
Z} % Split into its two rows. Gives two ×ばつN arrays
Zd % GCD, element-wise. Gives a ×ばつ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
 )
}

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

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

• Linnea Gräf: In the Ungolfed Version gcd calls the the function g
• aross: r=_=> is that code or a drawing?
• l4m2: 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: n=>(n*6/(q=_=>n--&&q(r=_=>Math.random()*1e6)+g(~r(),~r()))(g=(a,b)=>b?g(b,a%b):a>-2))**.5

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

• MickyT: If you're not worried about warnings, 1:x will work.

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

• TNT: 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.

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..106).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 1⁄2

• CalculatorFeline: That's 90 bytes and 88 chars...?
• maybeso: Thanks for catching that. I didn't count newlines, and while ², ³ are only 1 byte 6 is more. I fixed it to 89 bytes with no final newline.
• CalculatorFeline: You have not fixed the verbose code.
• maybeso: It's not a one-to-one match anyway with spacing, the quotes and numbers, etc.

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*6+Đ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*6+Đ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ÆÒL3öÃrjÃx)¬

Try it