The scheme described in the following paragraphs is fast but not very flat. For a flat distribution choose n2 random normal deviates as matrix elements and do Gram-Schmidt on the matrix (like this C code). A Scheme version of this.

This C code computes random orthogonal matrices and a few more tricks. The defined preprocessor symbol M is the number of dimensions. The M by M matrix named m is global and most routines read or modify it. I write here as if the matrix were a row of M orthonormal vectors.

twist(int i, int j, r theta) rotates the vectors about axes i and j by angle theta. Vectors orthogonal to both i and j remain fixed.

void to() tests m for orthogonality.

void pm(char * y) prints m with label y.

void stir() rotates the space about random pairs of axes by random amounts, 1024 times.

The above goes if sig=0. If 0≤N≤M and sig = 2N−1 then this chooses an element from SO(M−N, N). If (sig>>n)&1 then the nth dimension is imaginary, or better said the quadratic form that this ‘orthogonal’ matrix preserves as a transform is Σ(−1)Ki(xi)2 where Ki = 0 or 1 and sig = ΣKi2i.

Such a matrix transformation preserves the Minkowski (or Lorentz) metric when sig = 1. When sig = 0 it generates random elements of SO(M) but then the constant 1024 should be increased to perhaps 100000. There are faster and somewhat more complex ways to choose a random member of SO(M) that do not merely depend on a large number of random twists. The Lie group of Lorentz transforms, SO(3, 1), is not compact and there is thus no finite Haar measure on it and thus no sense of a uniform probability distribution over such transforms.

See this about applying these matrices to defining isometries in curved spaces. See this for a connection to Clifford algebras.