In MatLab I generate Haar distributed matrices like this:
m = randn(m,m)% all elements are normally distributed
u= qr(m) % make qr decomposition and what you get is Haar measure on "u"
So the math. statement is that if "m" is normally distributed, then "u" is Haar.
The reason is quite trivial - normal distribution is preserved by unitary transformations.
However righting this I begin to doubt myself about tiny details - depending how they implement qr algorithm the matrix u is not unique, it can be multiplied diag(+-1 ). Neverthelss most probably everything should be correct.