function U = randunitary(n, N)
% Generates uniformly random unitary matrices.
%
% function U = randunitary(n, N)
%
% U is a n-by-n-by-N array such that each slice U(:, :, i) is a random
% unitary matrix of size n (i.e., a matrix in the unitary group U(n)),
% sampled from the Haar measure (uniform distribution).
% 
% By default, N = 1.
%
% Complexity: N times O(n^3).
% For details on the algorithm, see Mezzadri 2007,
% "How to generate random matrices from the classical compact groups."
%
% See also: randrot qr_unique

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, June 18, 2019.
% Contributors: 
% Change log: 

    if nargin < 2
        N = 1;
    end
    
    if n == 1
        U = sign(randn(1, 1, N) + 1i*randn(1, 1, N));
        return;
    end
    
    % Generated as such, the slides of U are uniformly distributed over
    % U(n), the set of unitary matrices: see Mezzadri 2007, p597.
    U = qr_unique(randn(n, n, N) + 1i*randn(n, n, N));

end
