function U = rand_unitary(n, special)

%Generates a random n-dimensional unitary matrix U.
%U has determinant 1 if special==true.

C = 2*pi*(randn(n,n) + i*randn(n,n));
U = expm(C-C');
if(special)
	U = U/(det(U))^(1/n);
end