This question is not urgent; just a matter of curiosity...
It is relatively easy to generate an arbitrary 3D or even 4D rotation matrix using conjugation (i.e. YXY−1) of orthogonal rotations. I expect that there are ways to choose the contributing orthogonal angles of rotation in order to get a uniform random distribution of the resulting axis (and angle). (3D rotations are also related to quaternions.)
There's a wonderfully symmetrical 3D rotation matrix presentation, given in the first edition of D. F. Rogers and J. A. Adams book "Mathematical Elements for Computer Graphics". The elements of the matrix are:
u2+(1-u2)c; uv(1-c)+ws; uw(1-c)-vs
uv(1-c)-ws; v2+(1-v2)c; vw(1-c)+us
uw(1-c)+vs; vw(1-c)-us; w2+(1-w2)c
Here, (u,v,w) is a unit vector along the chosen axis of rotation, and s and c are the sine and cosine of the chosen angle of rotation.
[I mention this on my blog page.]
Now, a 4D rotation must be about a "2D-axis", or plane (where a 3D rotation is about a "1D-axis", or line).
I wonder if there's an equally elegant 4x4 matrix, in terms of a pair of mutually orthogonal unit vectors (defining the plane of rotation) and the sine and cosine of the angle of rotation in 4D.
Any on-line references regarding anything mentioned here would we gratefully received.