**Random Normal Deviate**: `rnd()` returns a random normal deviate.
This is a normal (Gaussian) distribution with zero mean and standard deviation of 1.
The pdf (probability density function) is e^{−x2/2}/√(2π).

**Matrix Multiply**: `mm( m, n, j, A, B, P)`
places the product of the m by n matrix A, and the n by j matrix B, in the space P for an m by j matrix.

**Gram-Schmidt**: For 0 < m ≤ n, `gs(m, n, A)` transforms the m by n matrix A, into an orthogonal matrix.
For any k ≤ m the first k rows of the resulting matrix span the same subspace as the first k rows of the original matrix.
For 0 < i, j ≤ m the inner product of row i and row j is δ_{ij}.

**Random Orthogonal Matrix**: For 0 < j ≤ k, `grom(j, k, A)` places a random orthogonal j by k matrix in the space A.
The probability distribution is invariant under the Haar measure of the Lie group O(k).

A demo.

clang demo.c lob.c ./a.out