Code golf a random orthogonal matrix An orthogonal matrix is a square matrix with real entries whose columns and rows are orthogonal unit vectors (i.e., orthonormal vectors). 

This means that M^T M = I, where I is the identity matrix and ^T signifies matrix transposition.

Note that this is orthogonal not "special orthogonal" so the determinant of M can be 1 or -1.

The aim of this challenge is not  machine precision so if M^T M = I to within 4 decimal places that will be fine.

The task is to write code that takes a positive integer n &gt; 1 and outputs a random orthogonal n by n matrix.  The matrix should be randomly and uniformly chosen from all  n by n orthogonal matrices.   In this context, "uniform" is defined in terms of the Haar measure, which essentially requires that the distribution does not change if multiplied by any freely chosen orthogonal matrix.  This means the values of the matrix will be floating point values in the range -1 to 1.

The input and output can be any form you find convenient.

Please show an explicit example of your code running.

You may not use any existing library function which creates orthogonal matrices. This rule is a little subtle so I  will explain more. This rule bans the use of any existing function which takes in some (or no) input and outputs a matrix of size at least n by n which is guaranteed to be orthogonal.  As an extreme example, if you want the n by n identity matrix, you will have to create it yourself.

You can use any standard random number generator library for choosing random numbers of your choosing.

Your code should complete within at most a few seconds for n &lt; 50.
Haskell, 169 150 148 141 132 131 bytes

import Numeric.LinearAlgebra
z=(unitary.flatten&lt;$&gt;).randn 1
r 1=asRow&lt;$&gt;z 1
r n=do;m&lt;-r$n-1;(&lt;&gt;diagBlock[m,1]).haussholder 2&lt;$&gt;z n


Recursively extend an orthogonal matrix of size n-1 by adding 1 to lower right corner and apply a random Householder reflection. randn gives a matrix with random values from a gaussian distribution and z d gives a uniformly distributed unit vector in d dimensions.

haussholder tau v returns the matrix I - tau*v*vᵀ which isn't orthogonal when v isn't a unit vector.

Usage:

*Main&gt; m &lt;- r 5
*Main&gt; disp 5 m
5x5
-0.24045  -0.17761   0.01603  -0.83299  -0.46531
-0.94274   0.12031   0.00566   0.29741  -0.09098
-0.02069   0.30417  -0.93612  -0.13759   0.10865
 0.02155  -0.83065  -0.35109   0.32365  -0.28556
-0.22919  -0.41411   0.01141  -0.30659   0.82575
*Main&gt; (&lt;1e-14) . maxElement . abs $ tr m &lt;&gt; m - ident 5
True