you are asking for the two-point correlation function of a Poisson process with unit density, which is just unity: $g(u)=1$.
the support for this is about as strong as for the Riemann zeroes: there is extensive numerical evidence but no conclusive theorem; see Soundarajan's 2006 paper cited above, or more recent papers on arXiv:0708.2567 and arXiv:1102.3648
the esssential difference between the function $g_R(u)=1-sinc^2(u)$ for the Riemann zeroes and $g(u)=1$ for the prime numbers, is that the former vanishes $\propto u^2$ for small $u$ ("level repulsion"), while the latter remains constant. This is the difference between the (conjectured) Gaussian unitary ensemble of Riemann zeroes and the Poisson ensemble of prime numbers. For large $u$ all correlations decay and $g_R(u)\rightarrow g(u)$.
First of all, note that there is no canonical notion of equidistribution on a countable set
like the integers.
When asking for the probability that $k$ 'randomly chosen' integers are coprime,
it is more-or-less intuitively clear what is meant by 'randomly chosen'.
This is basically because the density of integers divisible by a given prime $p$ is
(up to 'differences from rounding') the same in any interval $\{1, \dots, n\}$,
namely $1/p$.
In your question this is not the case: for example in $S_4 \cap \{1, \dots, 10\}$,
5 of 10 numbers (= 50 percent) are even, while in $S_4 \cap \{1, \dots, 10^6\}$,
1070 of 1273 integers are even (which is about 84 percent) and in
$S_4 \cap \{1, \dots, 10^{30}\}$, already 445064 of 462692 are even, which is
about 96 percent.
So in order to make your question well-defined you need to be explicit about what
you mean by 'randomly chosen'.
Best Answer
Well, first of all, $\pi$ is not just a random real number. Almost every real number is transcendental so how can we make the notion "$\pi$ is special" (in a number-theoretical sense) more precise?
Start by noticing that $$\pi=\int_{-\infty}^{\infty}\frac{dx}{1+x^2}$$ This already tells us that $\pi$ has something to do with rational numbers. It can be expressed as "a complex number whose real and imaginary parts are values of absolutely convergent integrals of rational functions with rational coefficients, over domains in $\mathbb{R}^n$ given by polynomial inequalities with rational coefficients." Such numbers are called periods. Coming back to the identity $$\zeta(2)=\frac{\pi^2}{6}$$ There is a very nice proof of this (that at first seems very unnatural) due to Calabi. It shows that $$\frac{3\zeta(2)}{4}=\int_0^1\int_0^1\frac{dx\,dy}{1-x^2y^2}$$ by expanding the corresponding geometric series, and then evaluates the integral to $\pi^2/8$. (So yes, $\pi^2$ and all other powers of $\pi$ are periods.) But the story doesn't end here as it is believed that there are truly deep connections between values of zeta functions (or L-functions) and certain evaluations involving periods, such as $\pi$. Another famous problem about primes is Sylvester's problem of which primes can be written as a sum of two rational cubes. So one studies the elliptic curve $$E_p: p=x^3+y^3$$ and one wants to know if there is one rational solution, the central value of the corresponding L-function will again involve $\pi$ up to some integer factor and some Gamma factor. Next, periods are also values of multiple zeta functions: $$\zeta(s_1,s_2,\dots,s_k)=\sum_{n_1>n_2>\cdots>n_k\geq 1}\frac{1}{n_1^{s_1}\cdots n_k^{s_k}}$$ And they also appear in other very important conjectures such as the Birch and Swinnerton-Dyer conjecture. But of course all of this is really hard to explain without using appropriate terminology, the language of motives etc. So, though, this answer doesn't mean much, it's trying to show that there is an answer to your question out there, and if you study a lot of modern number theory, it might just be satisfactory :-).