# 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} / {\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

• Are there any bounds on the range of inputs? – Eric Towers Dec 27 '15 at 9:24
• @EricTowers The program should work for any reasonable size up to limitations of memory and data type. At least 1000 – Luis Mendo Dec 27 '15 at 17:21
• Are rational numbers as return values (not strings) allowed? My language has a native rational type, in which 63/100 is a valid literal, usable in computation. (Langs I'm talking about: Factor, Racket) – cat Apr 25 '16 at 13:54
• @cat Sure, go ahead! Take into account the required precision, though – Luis Mendo Apr 25 '16 at 14:11

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

• Oh, Jelly includes the totient function!? Nice idea! – Luis Mendo Dec 26 '15 at 20:01
• Countdown until MATL includes a totient command at T-1 day... – quintopia Dec 27 '15 at 7:43
• @quintopia (I finally included the totient function) :-D – Luis Mendo Jan 19 '16 at 20:10

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

# Pyth - 13 11 bytes

.OqR1iM^SQ2


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

• Can you post an example run and explanation, for those of us who don't have Mathematica? – Luis Mendo Dec 26 '15 at 16:56
• This unites our submissions: Count[Array[GCD,{#, #}],1,2]/#^2.& Be my guest. – DavidC Dec 26 '15 at 21:40

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

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

• kron! Good idea! – Luis Mendo Dec 27 '15 at 4:22
• I first used meshgrid, but then I noticed I could do the same with kron inline! (->anonymous function) – flawr Dec 27 '15 at 9:31

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

• Very nice! You seem to be using Sage quite often lately :-) – Luis Mendo Apr 25 '16 at 9:53
• @LuisMendo Sage is great and does all things. It's very nice to use on math-based challenges because it has a huge builtin library like Mathematica, but the syntax is better (by virtue of a) not being Mathematica, and b) being built on Python). – Mego Apr 25 '16 at 19:52

# 05AB1E, 8 bytes

LÕO·<sn/


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, and I haven't really bothered trying to make it any smaller.

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

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

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

## Mathematica, 35 bytes

Implements @Dennis's algorithm.

(24Plus@@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 "14".) 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=(24Plus@@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[#] - 14)/#^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.

# 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

• Far as I can tell, you don't need to collect a lazy object to take its length. – cat Apr 25 '16 at 14:50
• @cat You do for some where length doesn't have a method defined, which is the case here for the filtered combinations iterator. Similarly endof wouldn't work because there's no method for that type to getindex. – Alex A. Apr 25 '16 at 20:56
• But the language disagrees. – cat Apr 25 '16 at 20:59
• @cat range doesn't return the same kind of object as combinations. The latter returns an iterator over combinations which is a separate type with no defined length method. See here. (Also the : syntax is preferred over range for constructing ranges ;)) – Alex A. Apr 25 '16 at 21:13

# MATL, 20 17 bytes

This uses current version (5.0.0) of the language.

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


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

• Couldnt you be even shorter with an approach like the one I used in octave? – flawr Dec 27 '15 at 0:11
• Indeed! Thank you! 3 bytes less. You should have done it yourself in MATL :-) – Luis Mendo Dec 27 '15 at 4:33
• I would have tried if it wasn't way past my bedtime=) – flawr Dec 27 '15 at 9:30

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

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


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

• Add an example run maybe? – Luis Mendo Apr 25 '16 at 16:43
• @LuisMendo done! – cat Apr 25 '16 at 17:06

# Whispers v2, 75 bytes

> 1
> Input
>> φ(L)
>> ∑ 1 2 3
>> 4+4
>> ≺5
>> 2²
>> 6÷7
>> Output 8


Try it online!

## How it works

Implements Dennis's algorithm.

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}\$$

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

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.

• That recycling feature is nice for golfing – Luis Mendo Sep 6 at 12:11
• You can save 3 bytes by replacing the end with N<-scan():1,rep(N,e=N). – Robin Ryder Sep 18 at 14:00