A probability problem by Gabriel Chacón

If $p_2, p_3, p_5, \ldots, p_n$ are the probabilities that $a$ and $b$ do not have the prime factors $2,3,5,\ldots,n$ in common respectively, the probability that they do not have any of these prime factors in common is:

$P=p_2 \cdot p_3 \cdot p_5 \cdot \ldots \cdot p_m \cdot \ldots \cdot p_s \cdot \ldots \quad (1)$

Let us first find the probability $q_m=1-p_m$ that $a$ and $b$ do have the common factor $m$ . Any integer number divided by $m$ yields as remainders the numbers $0,1,2,3,4,\ldots,(m-1)$ . When picking an integer at random, the probability of getting the remainder $0$ is $\frac{1}{m}$ . Therefore the probability that $a$ and $b$ both admit $m$ as a factor is:

$q_m=1-p_m = \dfrac{1}{m} \cdot \dfrac{1}{m}=\dfrac{1}{m^2} \quad \implies \quad p_m=1-\dfrac{1}{m^2}$

Expression $(1)$ then becomes:

$P=\left( 1-\dfrac{1}{2^2} \right) \left( 1-\dfrac{1}{3^2} \right) \left( 1-\dfrac{1}{5^2} \right) \ldots$

To find the infinite product, we work with the reciprocal of $P$ :

$\dfrac{1}{P}=\dfrac{1}{ 1-\frac{1}{2^2} }\cdot \dfrac{1}{ 1-\frac{1}{3^2} }\cdot \dfrac{1}{ 1-\frac{1}{5^2} }\cdot \dots$

Each of the factors in the expression above is the sum of a geometric progression:

$\dfrac{1}{P}=\left( 1+\dfrac{1}{2^2} +\dfrac{1}{(2^2)^2} + \ldots \right) \left( 1+\dfrac{1}{3^2} +\dfrac{1}{(3^2)^2} + \ldots \right) \left( 1+\dfrac{1}{5^2} +\dfrac{1}{(5^2)^2} + \ldots \right) \ldots$

Evaluate the product above to obtain:

$\dfrac{1}{P}= 1 + \dfrac{1}{2^2} + \dfrac{1}{3^2} + \dfrac{1}{4^2} + \dfrac{1}{5^2} + \ldots = \dfrac{\pi ^2}{6}$ (The first one to evaluate this sum was Leonhard Euler in the 18th century.)

Finally, the sought for probability is $\boxed{ P=\dfrac{6}{\pi ^2}=0.607927... }$

It is fair to say that neither the problem nor the solution are mine, of course. I found them in a 1948 book published in Madrid, Problemas de cálculo de probabilidades by J. Gallego Díaz, Ediciones Hispano-Argentinas. The problem is attributed to Chebyshev.

tested with b going up to 999 and every possible value for a given that value b, got a close enough answer to be considered correct

<?php
$t=0;
$c=0;
for($b=2;$b<=999;$b++){
for($a=1;$a<$b;$a++){
$c++;
if(gcd($a,$b)==1){
$t++;
}
}
}
echo $t/$c; //0.609

function gcd($a,$b) {
return ($a % $b) ? gcd($b,$a % $b) : $b;
}
?>