[Question]
  [
# Introduction
Number theory is full of wonders, in the form of unexpected connections. Here's one of them.
Two integers are [co-prime](https://en.wikipedia.org/wiki/Coprime_integers) 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} / {\pi^2} \approx 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 = {11} / {16} = 0.6875\$ for \$N = 4\$.
The *exact* value of \$p\$ needs to be computed with at least \$4\$ 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]
# [Jelly](http://github.com/DennisMitchell/jelly), ~~12~~ 8 bytes
```
RÆṪSḤ’÷²
```
[Try it online!](http://jelly.tryitonline.net/#code=UsOG4bmqU-G4pOKAmcO3wrI&input=&args=MTAwMA)
The following binary code works with [this version of the Jelly interpreter](https://github.com/DennisMitchell/jelly/tree/0ea6f804d01b9da15070c77e8e7ba3ffa0f9d114).
```
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]
# Mathematica ~~43~~ 42 bytes
I found [Lattice points visible from the origin](https://shreevatsa.wordpress.com/2008/11/07/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?
[![grid](https://i.stack.imgur.com/AUnvH.png)](https://i.stack.imgur.com/AUnvH.png)
---
```
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]
# Pyth - ~~13~~ 11 bytes
```
.OqR1iM^SQ2
```
[Test Suite](http://pyth.herokuapp.com/?code=.OqR1iM%5ESQ2&test_suite=1&test_suite_input=1%0A2%0A4%0A10&debug=0).
[Answer]
## 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]
# 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]
# 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]
## 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](http://sagecell.sagemath.org/?z=eJwdw0EKgzAQQNF9IHeY5Uw6RSd0JXiSUmTExAZMGtIKHl_p5_04bprnRUGHgvjdM2atGPYttKm-Ezcta0BhvQkROX8X6tQ5T9ZYEz8NDkgFnsKeHyz95X__GqyBq9pS-UHEg06BNByF&lang=sage)
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E), 8 bytes
```
LÕO·<sn/
```
[Try it online!](https://tio.run/##yy9OTMpM/f/f5/BU/0PbbYrz9P//NzQwMAAA "05AB1E – Try It Online")
Explanation:
```
L         Push inclusive range of (implicit) input (4 -> [1, 2, 3, 4])
 Õ        Perform Euler's Totient on each element ([1, 2, 3, 4] -> [1, 1, 2, 2])
  O       Sum list elements ([1, 1, 2, 2] -> 6)
   ·<     Double result and subtract 1 (2*6 - 1 = 11)
     sn   Square input (4 -> 16)
       /  Divide sum by input squared (11 / 16 = 0.6875)
          (Implicit output)
```
It might be golf-able, I just copied [Dennis' Jelly method](https://codegolf.stackexchange.com/a/67778/89342),
and I haven't really bothered trying to make it any smaller.
[Answer]
# 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]
## Seriously, 15 bytes
```
,;ª)u1x`▒`MΣτD/
```
Hex Dump:
```
2c3ba62975317860b1604de4e7442f
```
[Try It Online](http://seriouslylang.herokuapp.com/link/code=2c3ba62975317860b1604de4e7442f&input=1000)
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]
# 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](https://codegolf.stackexchange.com/a/67778/7311) gives a nice explanation how can we use the totient function.
[Try it online here.](http://tryj.tk/)
[Answer]
## 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]
# 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]
# [PARI/GP](http://pari.math.u-bordeaux.fr/), 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]
# [MATL](https://github.com/lmendo/MATL), ~~20~~ 17 bytes
This uses [current version (5.0.0)](https://github.com/lmendo/MATL/releases/tag/5.0.0) of the language.
Approach based on [@flawr's answer](https://codegolf.stackexchange.com/a/67818/36398).
```
i:tg!X*t!Zd1=Y)Ym
```
*Edit (April 28, 2015)*: [Try it online!](http://matl.tryitonline.net/#code=aTp0ZyFYKnQhWmQxPVg6WW0&input=MTAw) 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]
# [Samau](https://github.com/AlephAlpha/Samau/tree/master/OldSamau), 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]
# 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.
[Answer]
# Factor, ~~120~~ 113 bytes
I've spent class golfing this and I can't get it shorter.
Translation of: [Julia](https://codegolf.stackexchange.com/a/67876/46231).
```
[ [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]
# [Whispers v2](https://github.com/cairdcoinheringaahing/Whispers), 75 bytes
```
> 1
> Input
>> φ(L)
>> ∑ 1 2 3
>> 4+4
>> ≺5
>> 2²
>> 6÷7
>> Output 8
```
[Try it online!](https://tio.run/##K8/ILC5ILSo2@v/fTsGQy07BM6@gtITLzk7hfJuGjyaI8ahjooKhgpGCMYhjom0CFuvcZQqijQ5tAlFmh7ebg2j/0hKgbgWL//8NDQwMAA "Whispers v2 – Try It Online")
## How it works
Implements [Dennis](https://codegolf.stackexchange.com/users/12012/dennis)'s [algorithm](https://codegolf.stackexchange.com/a/67778/66833).
We calculate \$m = \sum\_{k=1}^n\phi(k)\$. We then double \$m\$ and decrement it to get \$2m-1\$, before squaring \$n\$ and dividing by it, to output \$\frac{2m-1}{n^2}\$
[Answer]
# [R](https://www.r-project.org/), 78 bytes
```
mean(mapply(function(x,y)sum(!x%%which(!y%%1:y))<2,1:(N=scan()),rep(1:N,e=N)))
```
[Try it online!](https://tio.run/##DctBCoAgEADAt3gQdsGD21HyC/5BxFBIE03K12/NfTpzib5C8a2dC45Zw52vCq9aOGYB8Ur5pBwSiCUlmYW4b4oMODvC/xBVjw3IOBWtQ0QmrTV/ "R – Try It Online")
Uses `which(!y%%1:y)` to calculate divisors of y (including 1), and `!x%%divisors_of_y` to count common divisors of x and y (including 1), so `sum(common_divisors)<2` is `TRUE` for co-prime pairs.
The R `mapply` function applies a custom function to two (or more) sets of arguments, here specified as `1:N` and `rep(1:N,each=N)` (which repeats each number in `1:N`, `N` times).
Conveniently, R 'recycles' the shorter argument list to the length of the longer one, so `1:N` gets repeated for us automaticaly.
]