Co-primality and the number pi

Question

Introduction

Number theory is full of wonders, in the form of unexpected connections. Here's one of them.

Two integers are co-prime if they have no factors in common other than 1. Given a number N, consider all integers from 1 to N. Draw two such integers at random (all integers have the same probability of being selected at each draw; draws are independent and with replacement). Let p denote the probability that the two selected integers are co-prime. Then p tends to 6/π² ≈ 0.6079... as N tends to infinity.

The challenge

The purpose of this challenge is to compute p as a function of N.

As an example, consider N = 4. There are 16 possible pairs obtained from the integers 1,2,3,4. 11 of those pairs are co-prime, namely (1,1), (1,2), (1,3), (1,4), (2,1), (3,1), (4,1), (2,3), (3,2), (3,4), (4,3). Thus p is 11/16 = 0.6875 for N = 4.

The exact value of p needs to be computed with at least four decimals. This implies that the computation has to be deterministic (as opposed to Monte Carlo). But it need not be a direct enumeration of all pairs as above; any method can be used.

Function arguments or stdin/stdout may be used. If displaying the output, trailing zeros may be omitted. So for example 0.6300 can be displayed as 0.63. It should be displayed as a decimal number, not as a fraction (displaying the string 63/100 is not allowed).

Winning criterion is fewest bytes. There are no restrictions on the use of built-in functions.

Test cases

Input / output (only four decimals are mandatory, as indicated above):

1    / 1.000000000000000
2    / 0.750000000000000
4    / 0.687500000000000
10   / 0.630000000000000
100  / 0.608700000000000
1000 / 0.608383000000000

Answer 1

Jelly, 12 8 bytes

RÆṪSḤ’÷²

Try it online!

The following binary code works with this version of the Jelly interpreter.

0000000: 52 91 b0 53 aa b7 9a 8a  R..S....

Idea

Clearly, the number of pairs (j, k) such that j ≤ k and j and k are co-prime equals the number of pairs (k, j) that satisfy the same conditions. Also, if j = k, j = 1 = k.

Thus, to count the number of co-prime pairs with coordinates between 1 and n, it suffices to calculate the amount m of pairs (j, k) such that j ≤ k, then compute 2m - 1.

Finally, since Euler's totient function φ(k) yields the number integers between between 1 and k that are co-prime to k, we can calculate m as φ(1) + … + φ(n).

Code

RÆṪSḤ’÷²    Input: n

R           Yield [1, ..., n].
 ÆṪ         Apply Euler's totient function to each k in [1, ..., n].
   S        Compute the sum of all results.
    Ḥ       Multiply the result by 2.
     ’      Subtract 1.
      ÷²    Divide the result by n².

Answer 2

Mathematica 43 42 bytes

I found Lattice points visible from the origin, from which the picture below is taken, to be helpful in reframing the problem through the following questions regarding a given square region of the unit lattice:

• What fraction of the unit-lattice points have co-prime coordinates?
• What fraction of unit-lattice points can be seen from the origin?

---

N@Mean[Mean/@Boole@Array[CoprimeQ,{#,#}]]&

---

Examples

Trailing zeros are omitted.

N@Mean[Mean/@Boole@Array[CoprimeQ,{#,#}]]&/@Range@10

{1., 0.75, 0.777778, 0.6875, 0.76, 0.638889, 0.714286, 0.671875, 0.679012, 0.63}

---

Timing

The timing, in seconds, precedes the answer.

N@Mean[Mean/@Boole@Array[CoprimeQ,{#,#}]]&[1000]// AbsoluteTiming

{0.605571, 0.608383}

Answer 3

Pyth - 13 11 bytes

.OqR1iM^SQ2

Test Suite.

Answer 4

Mathematica, 42 32 bytes

Count[GCD~Array~{#,#},1,2]/#^2.&

Anonymous function, uses simple brute force. The last case runs in about .37 seconds on my machine. Explanation:

                               &   A function taking N and returning
Count[               , , ]           the number of
                      1               ones
                                     in the
                        2             second
                                     level of
         ~Array~                      a matrix of
      GCD                              the greatest common denominators of
                {#,#}                 all 2-tuples in [1..N]
                          /         divided by
                           #          N
                            ^2.      squared.

Answer 5

Dyalog APL, 15 bytes

(+/÷⍴)∘∊1=⍳∘.∨⍳

Pretty straightforward. It's a monadic function train. Iota is the numbers from 1 to input, so we take the outer product by gcd, then count the proportion of ones.

Answer 6

Octave, 49 47 bytes

Just calculating the gcd of all pairs and counting.

@(n)mean(mean(gcd(c=kron(ones(n,1),1:n),c')<2))

The kronecker product is awesome.

Answer 7

MATL, 20 17 bytes

This uses current version (5.0.0) of the language.

Approach based on @flawr's answer.

i:tg!X*t!Zd1=Y)Ym

Edit (April 28, 2015): Try it online! After this answer was posted, function Y) was renamed to X:; the link includes that change.

Example

>> matl i:tg!X*t!Zd1=Y)Ym
> 100
0.6087

Explanation

i:         % vector 1, 2, ... up to input number
tg!        % copy, convert into ones, transpose
X*         % Kronecker product. Produces a matrix
t!         % copy, transpose
Zd         % gcd of all pairs
1=         % is it equal to 1?
Y)         % linearize into column array
Ym         % compute mean

---

Old answer: 20 bytes

Oi:t"t@Zd1=sb+w]n2^/

Explanation

O             % produce literal 0. Initiallizes count of co-prime pairs.
i             % input number, say N
:             % create vector 1, 2, ... N
t             % duplicate of vector
"             % for loop
    t         % duplicate of vector
    @         % loop variable: each element from vector
    Zd        % gcd of vector and loop variable. Produces a vector
    1=s       % number of "1" in result. Each "1" indicates co-primality
    b+w       % accumulate into count of co-prime pairs
]             % end
n2^/          % divide by N^2

Answer 8

JavaScript (ES6), 88 bytes

n=>eval(`p=0;for(i=n;i;i--)for(j=n;j;j--,p+=a)for(a=1,k=j;k>1;k--)a=i%k||j%k?a:0;p/n/n`)

Explanation

n=>
  eval(`                     // use eval to enable for loop without {} or return
    p=0;                     // p = number of pairs
    for(i=n;i;i--)           // i = 1 to n
      for(j=n;j;j--,p+=a)    // j = i to n, a will equal 1 if i and j are coprime, else 0
        for(a=1,k=j;k>1;k--) // for each number from 0 to j
          a=i%k||j%k?a:0;    // if i%k and j%k are both 0, this pair is not coprime
    p/n/n                    // return result (equivalent to p/(n*n))
  `)

Test

Takes a while for large (>100) values of n.

var solution = n=>eval(`p=0;for(i=n;i;i--)for(j=n;j;j--,p+=a)for(a=1,k=j;k>1;k--)a=i%k||j%k?a:0;p/n/n`)

N = <input type="text" id="input" value="100" />
<button onclick="result.textContent=solution(+input.value)">Go</button>
<pre id="result"></pre>

Answer 9

Seriously, 15 bytes

,;ª)u1x`▒`MΣτD/

Hex Dump:

2c3ba62975317860b1604de4e7442f

Try It Online

I'm not going to bother explaining it since it literally uses exactly the same algorithm as Dennis's Jelly solution (although I derived it independently).

Answer 10

Mathematica, 35 bytes

Implements @Dennis's algorithm.

(2`4Plus@@EulerPhi@Range[#]-1)/#^2&

Compute the sum of the totient (Euler's phi function) over the range from 1 to the input value. Multiply by floating point two (with four digits of precision) and subtract one. (More precision can be retained by using instead "2" and "1`4".) This is the total number of coprime pairs, so divide by the square of the input to get the desired fraction. Because the two is an approximate number, so is the result.

Testing (with timing data in the left column since at least one of us thinks that's interesting), with the function assigned to f so that the test line is more easily read.:

f=(2`4Plus@@EulerPhi@Range[#]-1)/#^2&
RepeatedTiming[f[#]] & /@ {1, 2, 4, 10, 100, 1000}
(* {{5.71*10^-6, 1.000}, 
    {5.98*10^-6, 0.750}, 
    {0.000010  , 0.6875}, 
    {0.0000235 , 0.6300}, 
    {0.00028   , 0.6087}, 
    {0.0033    , 0.6084} }  *)

Edit: Showing extent of the range of inputs (swapping the precision to the one instead of the two because otherwise the results get rather monotonous) and challenging others to do the same...

f = (2 Plus @@ EulerPhi@Range[#] - 1`4)/#^2 &
{#}~Join~RepeatedTiming[f[#]] & /@ {1, 2, 4, 10, 100, 1000, 10^4, 10^5, 10^6, 10^7}
(*  Results are {input, wall time, output}
    {{       1,  5.3*10^-6, 1.000}, 
     {       2,  6.0*10^-6, 0.7500}, 
     {       4,  0.0000102, 0.68750}, 
     {      10,  0.000023 , 0.63000}, 
     {     100,  0.00028  , 0.6087000}, 
     {    1000,  0.0035   , 0.608383000}, 
     {   10000,  0.040    , 0.60794971000}, 
     {  100000,  0.438    , 0.6079301507000}, 
     { 1000000,  4.811    , 0.607927104783000}, 
     {10000000, 64.0      , 0.60792712854483000}}  *)

RepeatedTiming[] performs the calculation multiple times and takes an average of the times, attempting to ignore cold caches and other effects causing timing outliers. The asymptotic limit is

N[6/Pi^2,30]
(*  0.607927101854026628663276779258  *)

so for arguments >10^4, we can just return "0.6079" and meet the specified precision and accuracy requirements.

Answer 11

J, 19 17 bytes

*:%~1-~5+/@p:|@i:

Usage:

   (*:%~1-~5+/@p:|@i:) 4
0.6875

Explanation:

*:%~1-~5+/@p:|@i:
               i: [-n..n]
             |@   absolute value of each element ([n..1,0,1,..n])
       5+/@p:     sum of the totient function for each element
    1-~           decreased by one, giving the number of co-prime pairs
*:%~              divided by N^2

Dennis' solution gives a nice explanation how can we use the totient function.

Try it online here.

Answer 12

Julia, 95 bytes

n->(c=combinations([1:n;1:n],2);((L=length)(collect(filter(x->gcd(x...)<2,c)))÷2+1)/L(∪(c)))

Just the straightforward approach for now; I'll look into shorter/more elegant solutions soon. This is an anonymous function that accepts an integer and returns a float. To call it, assign it to a variable.

Ungolfed:

function f(n::Integer)
    # Get all pairs of the integers from 1 to n
    c = combinations([1:n; 1:n], 2)

    # Get the coprime pairs
    p = filter(x -> gcd(x...) == 1, c)

    # Compute the probability
    return (length(collect(p)) ÷ 2 + 1) / length(unique(c))
end

Answer 13

Sage, 55 bytes

lambda a:n((sum(map(euler_phi,range(1,a+1)))*2-1)/a**2)

Thanks to Sage computing everything symbolically, machine epsilon and floating-point issues don't crop up. The tradeoff is, in order to follow the output format rule, an additional call to n() (the decimal approximation function) is needed.

Try it online

Answer 14

Samau, 12 bytes

Disclaimer: Not competing because I updated the language after the question was posted.

▌;\φΣ2*($2^/

Hex dump (Samau uses CP737 encoding):

dd 3b 5c ad 91 32 2a 28 24 32 5e 2f

Using the same algorithm as Dennis's answer in Jelly.

Answer 15

PARI/GP, 25 bytes

Making the function anonymous would save a byte, but then I'd have to use self making it more expensive overall.

f(n)=n^2-sum(j=2,n,f(n\j))

Answer 16

Factor, 120 113 bytes

I've spent class golfing this and I can't get it shorter.

Translation of: Julia.

[ [a,b] dup append 2 <combinations> [ members ] keep [ first2 coprime? ] filter [ length ] bi@ 2 /i 1 + swap /f ]

Example runs on the first 5 test cases (a value of 1000 causes it to freeze the editor, and I can't be bothered to compile an executable right now):

! with floating point division
IN: scratchpad auto-use {
      1    
      2    
      4    
      10   
      100  
    }
    [ 
      [1,b] dup append 2 <combinations> [ members ] keep 
      [ first2 coprime? ] filter [ length ] bi@ 2 /i 1 + swap /f 
    ]
    map

--- Data stack:
{ 1.0 0.75 0.6875 0.63 0.6087 }
! with rational division
IN: scratchpad auto-use {
      1    
      2    
      4    
      10   
      100  
    }
    [ 
      [1,b] dup append 2 <combinations> [ members ] keep 
      [ first2 coprime? ] filter [ length ] bi@ 2 /i 1 + swap / 
    ]
    map

--- Data stack:
{ 1.0 0.75 0.6875 0.63 0.6087 }
{ 1 3/4 11/16 63/100 6087/10000 }

Answer 17

Python2/Pypy, 178 bytes

The x file:

N={1:set([1])}
n=0
c=1.0
e=input()
while n<e:
 n+=1
 for d in N[n]:
  m=n+d
  try:N[m].add(d)
  except:N[m]=set([d,m])
 for m in range(1,n):
  if N[m]&N[n]==N[1]:c+=2
print c/n/n

Running:

$ pypy x <<<1
1.0
$ pypy x <<<10
0.63
$ pypy x <<<100
0.6087
$ pypy x <<<1000
0.608383

The code counts the co-prime pairs (n,m) for m<n twice (c+=2). This ignores (i,i) for i=1..n which is ok except for (1,1), thus being corrected by initialising the counter with 1 (1.0 to prepare for float division later) to compensate.