Unexpected distribution of ab-cd where a,b,c,d are independent and N(0,1) distributed I have discovered an unexpected curiosity that if
a,b,c,d are independent random variables and are N(0,1) distributed, then
|ab-cd| is exponentially distributed Exp(1). It's easy to verify that numerically.
Is there a proof of that ? or does it stem from some deeper fact ?
I noticed that it when I was considering cross products of random Gaussian distributed vectors and noticed that their components have \exp(-|x|/a) distribution.
Demonstration with python
N = 1000000
x = np.random.normal(size=(N,4))
z = np.abs(x[:,0] * x[:,1] - x[:,2] * x[:,3])
plt.hist(scipy.stats.expon.cdf(z), bins=100)



Writing
$$ab - cd = \left(\left[\left(\frac{a+b}{\sqrt 2}\right)^2 + \left(\frac{c+d}{\sqrt 2}\right)^2 \right] - \left[\left(\frac{a-b}{\sqrt 2}\right)^2  + \left(\frac{c-d}{\sqrt 2}\right)^2 \right]\right)/2$$
and noting that $(a+b, c+d, a-b, c-d)/\sqrt{2}$ has a standard Normal distribution, it is immediate (from the definitions of chi-square distributions arising from sums of squares of independent standard Normal variables) that the distribution of $ab -cd$ is twice that of the difference of two independent chi-squared(2) variables.
Since chi-squared(2) variables have Exponential$(1/2)$ distributions, $ab-cd$ is distributed as the difference of two Exponential$(1)$ distributions.  Since the characteristic function of an Exponential$(1)$ distribution is $\phi(t) = 1/(1 - it),$ the cf. of  the difference is
$$\phi(t)\phi(-t) = \frac{1}{1+t^2}.$$
The Fourier Transform of that is proportional to $f(\omega)=\exp(-|\omega|),$ showing that the absolute value $|ab-cd|$ must have a density proportional to $f(x)$ for $x\ge 0:$ that's the Exponential$(1)$ distribution.
The thread at Probability function for difference between two i.i.d. Exponential r.v.s gives a direct demonstration of this latter result via integration of the convolution.
Reference
(Courtesy Glen B: see comments to the question.)
The symmetry of standard Normal distributions around zero and the independence of $d$ from $(a,b,c)$ imply $(a,b,c,-d)$ has the same distribution as $(a,b,c,d),$ whence the distribution of $|ab-cd|$ is the same as that of $|ab-c(-d)|=|ab+cd|.$
A brief analysis of the question in this form appears in Stuart & Ord, Kendall's Advanced Theory of Statistics Volume I (3rd Ed. 1987) in the first half of exercise 11.21:

$x_r,y_r,$ $r=1,2,\ldots,$ are independent standardized normal variates.  Show that $z= x_1y_1 + x_2y_2$ is distributed exactly as $w_1-w_2,$ where the $w_j$ are independent with f.f. [probability density] $e^{-w_j}$ [NB: only for $w_j\ge 0$] and hence or otherwise that $z$ is distributed in the Laplace form $g(z)=\frac{1}{2}\exp\left(-|z|\right),$ and that $|z|$ is again exponentially distributed.

