Solve this probability problem symbolically

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

What's with all the heavy lifting and machinations?

d = ProductDistribution[{TriangularDistribution[{-1, 1}], 2}];
Probability[a^2 + b^2 > 1, {a, b} \[Distributed] d]

$\frac{19}{6}-\pi$

Finishes in a few seconds on a netbook...

This is not quick (includes J.M. comment):

pdf = UniformDistribution[2];
td = TransformedDistribution[(x - y)^2, {x, y} \[Distributed] pdf];
zd = TransformedDistribution[
   a + b, {a \[Distributed] td, b \[Distributed] td}];

then

ans = 1 - FullSimplify[CDF[zd,1]]

yields the desired result.