39
\$\begingroup\$

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??
\$\endgroup\$
  • 1
    \$\begingroup\$ 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? \$\endgroup\$ – Fatalize Oct 3 '16 at 13:41
  • \$\begingroup\$ @Fatalize if it is a limitation of the language, sure. Otherwise, you need to handle N=10^6. \$\endgroup\$ – NonlinearFruit Oct 3 '16 at 13:43
  • 1
    \$\begingroup\$ Related \$\endgroup\$ – Luis Mendo Oct 3 '16 at 14:01
  • 2
    \$\begingroup\$ The goal is misleading, it states only one pair of integers is checked. \$\endgroup\$ – user253751 Oct 4 '16 at 2:08
  • 1
    \$\begingroup\$ Does the upper limit to the random numbers generated need to be exactly 1000000? Would a larger upper limit be acceptable? \$\endgroup\$ – Sok Oct 4 '16 at 12:24

22 Answers 22

12
\$\begingroup\$

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
\$\endgroup\$
11
\$\begingroup\$

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

\$\endgroup\$
  • \$\begingroup\$ ḤRµȷ6Xµ€g2/ċ1÷³6÷½ saves 2 bytes. (ȷ6 is 10^6 in a single nilad, ċ1 counts ones) \$\endgroup\$ – PurkkaKoodari Oct 3 '16 at 15:44
  • \$\begingroup\$ 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) \$\endgroup\$ – Jonathan Allan Oct 3 '16 at 15:56
  • \$\begingroup\$ 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 \$\endgroup\$ – PurkkaKoodari Oct 3 '16 at 16:00
  • 1
    \$\begingroup\$ Ḥȷ6xX€ should work for the random sampling. \$\endgroup\$ – Dennis Oct 3 '16 at 16:08
9
\$\begingroup\$

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

\$\endgroup\$
  • \$\begingroup\$ According to OP, 1 <= i , j <= 10^6, so you need to use randrange(1,1e6+1). \$\endgroup\$ – mbomb007 Oct 3 '16 at 21:38
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – mbomb007 Oct 3 '16 at 21:42
  • \$\begingroup\$ @mbomb007 Good point, I've fixed it :) Been a while! \$\endgroup\$ – Kade Oct 4 '16 at 13:00
  • 1
    \$\begingroup\$ k=lambda:r.randrange(1e6)+1 saves two bytes \$\endgroup\$ – Jonathan Allan Oct 5 '16 at 19:05
  • 1
    \$\begingroup\$ @JonathanAllan good catch, thanks! \$\endgroup\$ – Kade Oct 6 '16 at 13:01
8
\$\begingroup\$

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&
\$\endgroup\$
  • 1
    \$\begingroup\$ You can save a byte: (6#/Count[GCD@@1*^6~RandomInteger~{2,#},1])^.5& \$\endgroup\$ – LegionMammal978 Oct 15 '16 at 15:42
  • 1
    \$\begingroup\$ Also, 1*^6 should be replaced with {1,1*^6} to ensure that i, j ≠ 0. \$\endgroup\$ – LegionMammal978 Oct 15 '16 at 15:47
8
\$\begingroup\$

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
\$\endgroup\$
  • 1
    \$\begingroup\$ You can simply use sample(10^6,N). Not only is it shorter, it is also much more efficient. \$\endgroup\$ – asac - Reinstate Monica Oct 3 '16 at 17:26
  • \$\begingroup\$ 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. \$\endgroup\$ – MickyT Oct 3 '16 at 19:02
  • \$\begingroup\$ @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. \$\endgroup\$ – rturnbull Oct 3 '16 at 19:21
  • \$\begingroup\$ Both methods are equally accurate when N<<10^6. To handle arbitrarily big N, you have to sample with replacement, good catch. \$\endgroup\$ – asac - Reinstate Monica Oct 4 '16 at 8:30
7
\$\begingroup\$

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
\$\endgroup\$
  • \$\begingroup\$ i, j aren't allowed to be 0 \$\endgroup\$ – isaacg Oct 3 '16 at 17:00
  • 1
    \$\begingroup\$ @isaacg They aren't. If you'll read the explanation, it says that the random values are selected from [0, 10**6), then incremented. \$\endgroup\$ – Mego Oct 3 '16 at 17:22
7
\$\begingroup\$

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
\$\endgroup\$
6
\$\begingroup\$

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
\$\endgroup\$
6
\$\begingroup\$

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
  )
}
\$\endgroup\$
  • \$\begingroup\$ I <3 Scala! Especially, because it sometimes really needs an explanation. \$\endgroup\$ – Roman Gräf Oct 3 '16 at 21:31
  • \$\begingroup\$ @RomanGräf To be honest, the only things I think might be unclear are 6f, Seq.fill and math.random. \$\endgroup\$ – corvus_192 Oct 3 '16 at 21:35
5
\$\begingroup\$

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
\$\endgroup\$
5
\$\begingroup\$

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));
}
\$\endgroup\$
  • \$\begingroup\$ In the Ungolfed Version gcd calls the the function g \$\endgroup\$ – Roman Gräf Oct 4 '16 at 5:54
  • \$\begingroup\$ r=_=> is that code or a drawing? \$\endgroup\$ – aross Oct 4 '16 at 14:36
  • \$\begingroup\$ 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 \$\endgroup\$ – l4m2 Dec 26 '17 at 10:00
  • \$\begingroup\$ n=>(n*6/(q=_=>n--&&q(r=_=>Math.random()*1e6)+g(~r(),~r()))(g=(a,b)=>b?g(b,a%b):a>-2))**.5 \$\endgroup\$ – l4m2 Dec 26 '17 at 10:05
5
\$\begingroup\$

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
\$\endgroup\$
4
\$\begingroup\$

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
\$\endgroup\$
4
\$\begingroup\$

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

\$\endgroup\$
  • 1
    \$\begingroup\$ If you're not worried about warnings, 1:x will work. \$\endgroup\$ – MickyT Oct 4 '16 at 20:30
4
\$\begingroup\$

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
\$\endgroup\$
3
\$\begingroup\$

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
\$\endgroup\$
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – TNT Oct 5 '16 at 17:25
1
\$\begingroup\$

Perl 6, 56 53 bytes

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

Try it online!

\$\endgroup\$
1
\$\begingroup\$

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 ½
\$\endgroup\$
  • \$\begingroup\$ That's 90 bytes and 88 chars...? \$\endgroup\$ – CalculatorFeline Jul 9 '17 at 21:50
  • \$\begingroup\$ 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. \$\endgroup\$ – maybeso Jul 11 '17 at 0:59
  • \$\begingroup\$ You have not fixed the verbose code. \$\endgroup\$ – CalculatorFeline Jul 11 '17 at 3:35
  • \$\begingroup\$ It's not a one-to-one match anyway with spacing, the quotes and numbers, etc. \$\endgroup\$ – maybeso Jul 11 '17 at 7:16
1
\$\begingroup\$

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.

\$\endgroup\$
1
\$\begingroup\$

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
\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$
1
\$\begingroup\$

Japt, 23 bytes

*6/UÆ1+1e6ö)j1+1e6öÃx)¬

Try it

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.