function [u, p] = polardecomp(a, side)
    % Compute the polar decomposition.
    %
    % Returns the factors of the polar decomposition [1] `u` and `p` such
    % that ``a = up`` (if `side` is "right") or ``a = pu`` (if `side` is
    % "left"), where `p` is positive semidefinite. Depending on the shape
    % of `a`, either the rows or columns of `u` are orthonormal. When `a`
    % is a square array, `u` is a square unitary array. When `a` is not
    % square, the "canonical polar decomposition" [2] is computed.
    %
    % Parameters:
    %   a: The array to be factored.
    %   side: (optional) Determines whether a right or left polar decomposition
    %         is computed. If `side` is "right", then ``a = up``. If `side` is
    %         "left", then ``a = pu``. The default is "right".
    %
    % Returns:
    %   u: If `a` is square, then `u` is unitary. If m > n, then the columns
    %      of `a` are orthonormal, and if m < n, then the rows of `u` are
    %      orthonormal.
    %   p: `p` is Hermitian positive semidefinite. If `a` is nonsingular, `p`
    %      is positive definite. The shape of `p` is (n, n) or (m, m), depending
    %      on whether `side` is "right" or "left", respectively.
    
    if nargin < 2
        side = 'right';
    end
    
    if ~strcmp(side, 'right') && ~strcmp(side, 'left')
        error('`side` must be either ''right'' or ''left''');
    end
    
    if ~ismatrix(a)
        error('`a` must be a 2-D array.');
    end
    
    [w, s, vh] = svd(a, 'econ');
    u = w * vh';
    
    if strcmp(side, 'right')
        % a = up
        p = (vh' * diag(diag(s))) * vh;
    else
        % a = pu
        p = (w * diag(diag(s))) * w';
    end
end

%!test
%! a = [1,-1; 2, 4];
%! [u, p] = polardecomp(a);
%! u0 = [0.85749293,  -0.51449576 ;  0.51449576,   0.85749293];
%! p0 = [1.88648444,  1.2004901; 1.2004901,  3.94446746];
%! assert(u,u0, 1e-8)
%! assert(p,p0, 1e-8)
%! assert(u*u', eye(2), 1e-15)