function Q = randrot(n, N)
% Generates uniformly random rotation matrices.
%
% function Q = randrot(n, N)
%
% Q is an n-by-n-by-N array such that each slice Q(:, :, i) is a random
% orthogonal matrix of size n of determinant +1 (i.e., a matrix in SO(n)),
% sampled from the Haar measure (uniform distribution).
%
% By default, N = 1.
%
% Complexity: N times O(n^3).
% Theory in Diaconis and Shahshahani 1987 for the uniformity on O(n);
% With details in Mezzadri 2007,
% "How to generate random matrices from the classical compact groups."
%
% To ensure matrices in SO(n), we permute the two first columns when
% the determinant is -1.
%
% See also: randskew qr_unique randunitary

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Sept. 25, 2012.
% Contributors: 
% Change log:
%   June 18, 2019 (NB)
%       Now generating all initial random matrices in one shot (which
%       should be more efficient) and calling qr_unique.


    if nargin < 2
        N = 1;
    end
    
    if n == 1
        Q = ones(1, 1, N);
        return;
    end
    
    % Generated as such, Q is uniformly distributed over O(n): the group
    % of orthogonal matrices; see Mezzadri 2007.
    Q = qr_unique(randn(n, n, N));
    
    for k = 1 : N
        
        % If a slice of Q is in O(n) but not in SO(n), we permute its two
        % first columns to negate its determinant. This ensures the new
        % slice is in SO(n), uniformly distributed.
        if det(Q(:, :, k)) < 0
            Q(:, [1 2], k) = Q(:, [2 1], k);
        end
        
    end

end
