double det( matrix A) /* determinant */
{
  int j, n = cols(A), sign = 1;  double r = 0;  vector k;

  if( n == 1) return A(1,1);

  for( j = 1; j <= n; ++j) /* go across row 1 */
  {
    if( j == 1) k = [2:n]; /* select all columns except column j */
    else if( j == n) k = [1:n-1];
    else k = [1:j-1,j+1:n];

    r += sign*A(1,j)*det(A(2:n,k)); /* add on sub-determinant */

    sign = -sign;
  }

  return r;
}
