## Copyright (C) 2017 Michael Leitner
## Copyright (C) 2001-2012 Paul Kienzle
## Copyright (C) 2009 VZLU Prague
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {} perms (@var{v})
##
## Generate all permutations of @var{v}, one row per permutation.  The
## result has size @code{factorial (@var{n}) * @var{n}}, where @var{n}
## is the length of @var{v}.
##
## As an example, @code{perms([1, 2, 3])} returns the matrix
##
## @example
## @group
##   1   2   3
##   2   1   3
##   1   3   2
##   2   3   1
##   3   1   2
##   3   2   1
## @end group
## @end example
## @end deftypefn

#CHECKME should in principle be more efficient when done with uint8 index type, but is not. perhaps in later octave versions integer addressing will become more efficient, then this should be changed here to get large increases in efficiency.

function A = perms (v)
  if (nargin != 1)
    print_usage ();
  endif
  v = sort(v(:))(end:-1:1)';
  n = length (v);

  if (n == 0)
    A = [];
  elseif (n==1)
    A=v;
  else
    --n;
    
    idx=(zeros(factorial(n),n));
    
    idx(1,n)=n;
    
    f=1;
    for j=1:n-1
      b=1:n;
      f*=j;
      perm=idx(1:f,n-(j-1):n);
      idx(1:(j+1)*f,n-j)=(n-j:n)(ones(f,1),:)(:);
      for i=0:j
        --b(i+n-j);
        idx((1:f)+i*f,n-(j-1):n)=b(perm);
      endfor
    endfor
    
    ++n;
    f*=n-1;
    A= zeros (factorial(n), n, class (v));
    A(:,1)=v(ones(f,1),:)(:);
    
    for i=1:n
      b=v([1:i-1 i+1:n]);
      A((1:f)+(i-1)*f,2:end)=b(idx);
    end

  endif
endfunction

%!error perms ();
%!error perms (1, 2);

%!assert (perms ([1,2,3]), [1,2,3;2,1,3;1,3,2;2,3,1;3,1,2;3,2,1]);
%!assert (perms (1:3), perms ([1,2,3]));

%!assert (perms (int8([1,2,3])), int8([1,2,3;2,1,3;1,3,2;2,3,1;3,1,2;3,2,1]));
