This is an interesting problem, because the difficulty is not the concept, but rather how to compute it (efficiently). Given points $(X_i,Y_i)$ distributed Uniformly on the unit square, we are interested in $$P\big[\sqrt{(X_2-X_1)^2 + (Y_2-Y_1)^2} \; > \; 1\big] $$
Let $X = X_2 - X_1$ denote the difference of two standard Uniform random variables, which is well-known to be a symmetric Triangular distribution on (-1,1). Similarly, let $Y = Y_2-Y_1$. By independence, the joint pdf of $(X,Y)$, say $f(x,y)$ is then:
f = (1-Abs[x]) (1-Abs[y]); domain[f] = {{x,-1,1}, {y,-1,1}};
We seek:

(source: tri.org.au)
All done. This takes just 2 seconds to evaluate, starting from a fresh kernel, where I am using the Prob
function from the mathStatica package for Mathematica (how I roll, being one of the authors), but one can equally use in-built Mathematica functions to the same effect:
dist = ProbabilityDistribution[(1 - Abs[x]) (1 - Abs[y]), {x, -1, 1}, {y, -1, 1}];
Probability[Sqrt[x^2 + y^2] > 1, Distributed[{x, y}, dist]]
(NIntegrate[Boole[(x1 - x2)^2 + (y1 - y2)^2 > 1], {x1, 0, 1}, {y1, 0, 1}, {x2, 0, 1}, {y2, 0, 1}] + Pi // RootApproximant) - Pi
evaluates to19/6 - Pi
$\endgroup$