function value = determinant ( A )
%
%  function value = determinant ( A )
%
%% DETERMINANT computes the determinant of a matrix using the definition.
%
%  Discussion:
%
%    For a square matrix A, the determinant is defined as 
%
%      det ( A ) = sum ( all permutations P ) 
%        Sign(P) * Product ( 1 <= I <= N ) A(I,P(I))
%
%    The determinant is very expensive and time-consuming to compute,
%    even for relatively small values of N.
%
%  Modified:
%
%    01 February 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(N,N), a matrix whose determinant is desired.
%
%    Output, real VALUE, the determinant of A.
%
[m,n] = size ( A );

if ( m ~= n )
  fprintf ( '\n' );
  fprintf ( 'DETERMINANT - Fatal error!\n' );
  fprintf ( '  The matrix is not square!\n' );
  value = 0.0;
  return
end

p = zeros ( 1, n );
nfact = prod ( [1:n] );

value = 0.0;
for i = 1 : nfact
  p = next_perm ( p );
  s = perm_sign ( p );
  x = diag ( A([1:n],p) );
  value = value + s * prod ( x );
end

function q = next_perm ( p )
%
%***********************************************************************
%
%  function q = next_perm ( p )
%
%% NEXT_PERM computes the lexicographic successor of a permutation.
%
%  Discussion:
%
%    Before calling this routine the first time, set
%
%      p = zeros ( n, 1 );
%
%    You may then make N factorial calls to this routine:
%
%      p = next_perm ( p )
%
%    and each time you will receive a new permutation.  
%    
%  Reference:
%
%    Algorithm 2.14,
%    Donald Kreher and Douglas Simpson,
%    Combinatorial Algorithms,
%
%  Modified:
%
%    01 February 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P(N), describes the permutation.
%    P(I) is the item which is permuted into the I-th place
%    by the permutation.  
%
%    Output, integer Q(N), describes the next permutation.  If
%    P was the very last permutation, then Q starts over with
%    the first one.
%
  n = length ( p );

  q = p;
%
%  Return the only permutation.
%
  if ( n == 1 ) 

    q = 1;
%
%  Return the first permutation.
%
  else if ( q == 0 ) 

    q = [1:n];
%
%  Return the next permutation.
%
  else
%
%  Identify the last index I for which the permutation value increases.
%
    i = n - 1;

    while ( q(i) > q(i+1) )
      i = i - 1;
      if ( i == 0 ) 
        break;
      end
    end
%
%  If no such I, we've reached the last permutation.
%
    if ( i == 0 ) 

      q = [1:n];
%
%  Seek the last index J whose permutation value is greater than I's.
%
    else

      j = n;
      while ( q(j) < q(i) )
        j = j - 1;
      end
%
%  Swap elements I and J.
%
      t = q(j);
      q(j) = q(i);
      q(i) = t;
%
%  Reverse elements I+1 through N.
%
      q(i+1:n) = q(n:-1:i+1);

    end

  end

end

function p_sign = perm_sign ( p )
%
%  function p_sign = perm_sign ( p )
%
%% PERM_SIGN returns the sign of a permutation.
%
%
%  Discussion:
%
%    A permutation can always be replaced by a sequence of pairwise
%    transpositions.  A given permutation can be represented by
%    many different such transposition sequences, but the number of
%    such transpositions will always be odd or always be even.
%    If the number of transpositions is even or odd, the permutation is
%    said to be even or odd.
%
%  Example:
%
%    Input:
%
%      P = 2, 3, 9, 6, 7, 8, 5, 4, 1
%
%    Output:
%
%      P_SIGN = +1
%
%  Reference:
%
%    A Nijenhuis and H Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Modified:
%
%    01 February 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P(N), describes the permutation.
%    The object in position I was permuted to position P(I).
%
%    Output, integer P_SIGN, the "sign" of the permutation.
%    +1, the permutation is even,
%    -1, the permutation is odd.
%
  n = length ( p );

  p_sign = 1;
%
%  Put each item 1 through N-1 back in its original position.
%
  for i = 1 : n-1
%
%  J is the current position of item I.
%
    j = i;

    while ( p(j) ~= i )
      j = j + 1;
    end
%
%  Unless the item is already in the correct place, restore it.
%
    if ( j ~= i )

      temp = p(i);
      p(i) = p(j);
      p(j) = temp;

      p_sign = - p_sign;

    end

  end