/*******************************************************************************
+
+  LEDA 3.5
+
+  _integer_matrix.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/
//---------------------------------------------------------------------
// file generated by notangle from integer_matrix.lw
// please debug or modify LEDA web file
// mails and bugs: leda@mpi-sb.mpg.de
// based on LEDA architecture by S. Naeher, C. Uhrig
// coding: K. Mehlhorn, M. Seel
//---------------------------------------------------------------------


#include <LEDA/integer_matrix.h>


template<class T>
inline void swap(T& x, T& y) 
{ 
  T h = x; 
  x = y; 
  y = h; 
}


integer_matrix::
integer_matrix(int dim1, int dim2)  
{ 
  LEDA_OPT_PRECOND((dim1>=0 && dim2>=0), "integer_matrix::constructor: \
  negative dimension.") 

  d1=dim1; 
  d2=dim2; 

  if (d1 > 0) { 
    allocate_mat_space(v,d1);
    for (int i=0; i<d1; i++) 
      v[i] = new integer_vector(d2); 
  }
  else 
    v = nil; 
}


integer_matrix::
integer_matrix(int dim)  
{ 
  LEDA_OPT_PRECOND((dim >= 0), "integer_matrix::constructor: \
  negative dimension.") 

  d1 = d2 = dim;

  if (d1 > 0)
  { 
    allocate_mat_space(v,d1);
    for (int i=0; i<d1; i++) 
      v[i] = new integer_vector(d2); 
  }
  else 
    v = nil; 
}


integer_matrix::
integer_matrix(const array< integer_vector >& A)  
{  
  int al = A.low(); 
  d2 = A.high() - al + 1; 
  d1 = A[al].dim(); 
    
  if (d1 > 0) 
  { 
    allocate_mat_space(v,d1);
    for (int i=0; i<d1; i++) 
    { 
      v[i] = new integer_vector(d2); 
      for (int j = 0; j < d2; j++) 
        v[i]->v[j] = A[al + j][i]; 
    }
  }
  else 
    v = nil; 
}


integer_matrix::
integer_matrix(const integer_matrix& p)  
{ 
  d1 = p.d1; 
  d2 = p.d2; 
    
  if (d1 > 0) 
  {  
    allocate_mat_space(v,d1);
    for (int i=0; i<d1; i++) 
      v[i] = new integer_vector(*p.v[i]); 
  }
  else 
    v = nil; 
}


integer_matrix::
integer_matrix(const integer_vector& vec)
{   
  d1 = vec.d; 
  d2 = 1; 
  if (d1>0)
    allocate_mat_space(v,d1);
  else
    v = nil;
  for(int i = 0; i < d1; i++)
  { 
    v[i] = new integer_vector(1); 
    elem(i,0) = vec[i]; 
  }
}


integer_matrix::
integer_matrix(int dim1, int dim2, RTINT** p)  
{ 
  LEDA_OPT_PRECOND((dim1 >= 0 && dim2 >= 0), "integer_matrix::constructor: \
  negative dimension.")
  d1=dim1; 
  d2=dim2; 

  if (d1 > 0)
  {
    allocate_mat_space(v,d1);
    for(int i=0; i<d1; i++) 
    { 
      v[i] = new integer_vector(d2); 
      for(int j=0; j<d2; j++) 
        elem(i,j) = p[i][j]; 
    }
  }
  else 
    v = nil;
}


integer_matrix& integer_matrix::
operator=(const integer_matrix& mat)
{ 
  register int i,j; 

  if (d1 != mat.d1 || d2 != mat.d2)
  { 
    for(i=0; i<d1; i++) 
      delete v[i]; 
    if (v)
      deallocate_mat_space(v,d1);

    d1 = mat.d1; 
    d2 = mat.d2; 

    if (d1>0)
      allocate_mat_space(v,d1);
    for(i = 0; i < d1; i++) 
      v[i] = new integer_vector(d2); 
  }

  for(i = 0; i < d1; i++)
    for(j = 0; j < d2; j++) 
      elem(i,j) = mat.elem(i,j); 
  return *this; 
}


integer_matrix::
~integer_matrix()  
{ 
  if (v) {
    for (int i=0; i<d1; i++) 
      delete v[i];  
    deallocate_mat_space(v,d1);
  }
}


integer_matrix integer_matrix::
identity(int n)
{
  integer_matrix M(n);
  for (int i=0; i<n; i++)
    M(i,i) = (RTINT)1;
  return M;
}


integer_vector integer_matrix::
col(int i)  const
{ 
  LEDA_OPT_PRECOND((i>=0 && i<d2), "integer_matrix::col:\
  index out of range.") 

  integer_vector result(d1); 
  int j = d1; 
  while (j--) 
    result.v[j] = elem(j,i); 
  return result; 
}


integer_vector integer_matrix::
to_vector() const
{ 
  LEDA_OPT_PRECOND((d2==1), "integer_matrix::to_vector: \
  cannot make vector from matrix.") 
  return col(0); 
}


int integer_matrix::
operator==(const integer_matrix& x) const
{ 
  register int i,j; 
  if (d1 != x.d1 || d2 != x.d2) 
    return false; 

  for(i = 0; i < d1; i++)
    for(j = 0; j < d2; j++)
      if (elem(i,j) != x.elem(i,j)) 
        return false; 
  return true; 
}


integer_matrix integer_matrix::
operator+ (const integer_matrix& mat)
{ 
  register int i,j; 
  check_dimensions(mat); 
  integer_matrix result(d1,d2); 
  for(i=0; i<d1; i++)
    for(j=0; j<d2; j++)
      result.elem(i,j) = elem(i,j) + mat.elem(i,j); 
  return result; 
}


integer_matrix integer_matrix::
operator- (const integer_matrix& mat)
{ 
  register int i,j; 
  check_dimensions(mat); 
  integer_matrix result(d1,d2); 
  for(i=0; i<d1; i++)
    for(j=0; j<d2; j++)
      result.elem(i,j) = elem(i,j) - mat.elem(i,j); 
  return result; 
}


integer_matrix integer_matrix::
operator- ()  // unary
{ 
  register int i,j; 
  integer_matrix result(d1,d2); 
  for(i=0; i<d1; i++)
    for(j=0; j<d2; j++)
      result.elem(i,j) = -elem(i,j); 
  return result; 
}


integer_matrix& integer_matrix::
operator-= (const integer_matrix& mat) 
{ 
  register int i,j; 
  check_dimensions(mat); 
  for(i=0; i<d1; i++)
    for(j=0; j<d2; j++)
      elem(i,j) -= mat.elem(i,j); 
  return *this; 
}


integer_matrix& integer_matrix::
operator+= (const integer_matrix& mat) 
{ 
  register int i,j; 
  check_dimensions(mat); 
  for(i=0; i<d1; i++)
    for(j=0; j<d2; j++)
      elem(i,j) += mat.elem(i,j); 
  return *this; 
}


integer_matrix integer_matrix::
operator*(const integer_matrix& mat) const
{ 
  LEDA_PRECOND((d2==mat.d1), "integer_matrix::operator*: \
  incompatible matrix types.") 
  
  integer_matrix result(d1, mat.d2); 
  register int i,j; 

  for (i=0; i<mat.d2; i++)
  for (j=0; j<d1; j++) 
    result.elem(j,i) = *v[j] * mat.col(i); 
  return result; 
}


integer_matrix integer_matrix::
compmul(const RTINT& f) const
{ 
  register int i,j; 
  integer_matrix result(d1,d2); 
  for(i=0; i<d1; i++)
    for(j=0; j<d2; j++)
      result.elem(i,j) = elem(i,j) *f; 
  return result; 
}


ostream&  
operator<<(ostream& O, const integer_matrix& M)
{ 
  /* syntax: d1 d2 
             x_0,0  ... x_0,d1-1
                  d2-times
             x_d2,0 ... x_d2,d1-1 */
  O << M.d1 << ' ' << M.d2 << ' ' << '\n';
  for (register int i=0; i<M.d1; i++) {
    for (register int j=0; j<M.d2; j++) 
      O << M(i,j) << " ";
    if (i<(M.d1-1)) O << "\n";
  }
  return O; 
}

istream&  
operator>>(istream& in, integer_matrix& M) 
{ 
  /* syntax: d1 d2 
             x_0,0  ... x_0,d1-1
                  d2-times
             x_d2,0 ... x_d2,d1-1 */

  int dim1 = 0, dim2 = 0;
  if (!(in >> dim1 >> dim2)) 
    return in;
  if (M.d1 != dim1 || M.d2 != dim2)
    M = integer_matrix(dim1,dim2);

  for (register int i=0; i<dim1; i++)
    for (register int j=0; j<dim2; j++) 
      if (!(in >> M(i,j)))
        return in;
  return in; 
}


int integer_matrix::
cmp(const integer_matrix& M1, const integer_matrix& M2) 
{ 
  register int i;
  int res;
  M1.check_dimensions(M2);
  for(i=0; i < M1.dim1() && 
           (res = compare(M1.row(i),M2.row(i))) != 0; i++);
  return res;
}
 

bool integer_matrix::
gauss_solver(const integer_vector& b, 
             integer_vector& x, 
             RTINT& D, 
             integer_matrix& spanning_vectors, 
             integer_vector& c) const
{ 
  const integer_matrix& A = *this;
  bool solvable = true; 
  
    int i,j,k; // indices to step through the matrix
    int rows = A.dim1(); 
    int cols = A.dim2(); 

    /* at this point one might want to check whether the computation can
       be carried out with doubles, see section \ref{ optimization }.
    */

    if (b.dim() != rows) 
      error_handler(1,"linear_solver: b has wrong dimension"); 

    integer_matrix C(rows,cols + 1); 
    // the matrix in which we will calculate ($C = (A\vert b)$)
    
    /* copy |A| and |b| into |C| and L becomes the identity matrix */
    integer_matrix L(rows); // zero initialized
    for(i=0; i<rows; i++) { 
      for(j=0; j<cols; j++) 
        C(i,j)=A(i,j); 
      C(i,cols)=b[i]; 
      L(i,i) = 1; // diagonal elements are 1
    }


  
    array<int> var(0,cols - 1); 
    // column $j$ of |C| represents the |var[j]| - th variable
    // the array is indexed between |0| and |cols - 1|

    for(j=0; j<cols; j++)
      var[j]= j; // at the beginning, variable $x_j$ stands in column $j$

    RTINT denom = 1; // the determinant of an empty matrix is 1
    int sign = 1; // no interchanges yet
    int rank = 0; // we have not seen any non-zero row yet

    /* here comes the main loop */
    for(k=0; k<rows; k++)
    {   
      
         bool non_zero_found = false; 
         for(i =k; i<rows; i++)  // step through rows $k$ to |rows - 1|
         {    
           for (j = k ; j < cols && C(i,j) == 0; j++) ; 
           // step through columns |k| to |cols - 1|
           if (j < cols) 
           {  
             non_zero_found = true; 
             break; 
           }
         }


      if (non_zero_found)
      { 
        rank++; //increase the rank
        
          if (i != k) 
          { 
            sign = - sign; 
            /* we interchange rows |k| and |i| of |L| and |C| */
            swap(L[i],L[k]); 
            swap(C[i],C[k]); 
          }

          if (j != k)
          { 
            sign = - sign; 
            /* We interchange columns |k| and |j| */
            for(int ih = 0; ih < rows; ih++) 
              swap(C(ih,k),C(ih,j)); 

            /* We store the exchange of variables in |var| */
            swap(var[k],var[j]); 
          }
         
          
          for(i = k + 1; i < rows; i++)  
            for (j = 0; j <  rows; j++)  //and all columns of |L|
              L(i,j) = (L(i,j)*C(k,k) - C(i,k)*L(k,j))/denom; 

          for(i = k + 1; i < rows; i++) 
          {      
            /* the following iteration uses and changes |C(i,k)| */
            RTINT temp = C(i,k); 
            for (j = k; j <= cols; j++)
              C(i,j) = (C(i,j)*C(k,k) - temp*C(k,j))/denom; 
          }
          denom = C(k,k); 
          
          #ifdef LA_SELFTEST

            for(i = 0; i < rows; i++)
            { 
              for (j = 0; j < cols; j++)
              { 
                RTINT Sum = 0; 
                for (int l = 0; l < rows; l++) 
                  Sum += L(i,l)*A(l, var[j]); 
                if (Sum  != C(i,j)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 
              }
              RTINT Sum = 0; 
              for (int l = 0; l < rows; l++) 
                Sum += L(i,l)*b[l];
                if (Sum  != C(i,cols)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 

            }
          #endif


      }
      else 
        break; 
    }


  
    for(i =rank ; i <rows && C(i,cols) == 0; i++); // no body

    if (i < rows) 
    { 
      solvable = false; 
      c = integer_vector(rows); 
          
      for (j = 0; j <rows; j++) 
        c[j] = L(i,j); 
    }

  if (solvable) 
  { 
    
      x = integer_vector(cols); 
      D = denom; 
      for(i=rank - 1 ; i>= 0; i--) 
      { 
        RTINT h =C(i,cols) * D; 
        for (j = i + 1; j<rank; j++)  
        { 
          h -= C(i,j)*x[var[j]]; 
        }
        x[var[i]]= h / C(i,i); 
      }

    #ifdef LA_SELFTEST
      /* we check whether |x| is a solution */
      { 
        for (i = 0; i < rows; i++)
        { 
          RTINT sum = 0; 
          for (j =0; j < cols; j++) 
            sum += A(i,j)*x[j]; 
          if (sum != D * b[i])
            error_handler(1,"linear_solver: base is not a solution"); 
        }
      }
    #endif

      int dimension = cols - rank; //dimension of solution
      spanning_vectors = integer_matrix(cols,dimension); 
     
      if (dimension > 0)
      { 
     
    /* In the $l$ - th spanning vector, $0 \le l < |dimension|$ we set
    variable |var[rank + l]| to $1 = |denom|/|denom|$ and then the
    dependent variables as dictated by the $|rank| + l$ - th column of
    |C|.*/

        for(int l=0; l < dimension; l++) 
        { 
          spanning_vectors(var[rank + l],l)=D; 
          for(i = rank - 1; i >= 0 ; i--)
          { 
            RTINT h = - C(i,rank + l)* D; 
            for ( j= i + 1; j<rank; j++)  
              h -= C(i,j)*spanning_vectors(var[j],l); 
            spanning_vectors(var[i],l)=h / C(i,i); 
          }

    #ifdef LA_SELFTEST
          /* we check whether the $l$ - th spanning vector is a solution 
             of the homogeneous system */
          { 
            integer_vector zero(rows); 
            if (A *spanning_vectors.col(l) != zero)
              error_handler(1,"linear_solver: spanning_vector is not a solution."); 
          }
    #endif
        }
      }

  }
  return solvable; 
}


RTINT  
determinant(const integer_matrix& A)
{ 
  if (A.dim1() != A.dim2())
    error_handler(1,"determinant: only square matrices are legal inputs."); 
  integer_vector b(A.dim1()); // zero - vector
  
    int i,j,k; // indices to step through the matrix
    int rows = A.dim1(); 
    int cols = A.dim2(); 

    /* at this point one might want to check whether the computation can
       be carried out with doubles, see section \ref{ optimization }.
    */

    if (b.dim() != rows) 
      error_handler(1,"linear_solver: b has wrong dimension"); 

    integer_matrix C(rows,cols + 1); 
    // the matrix in which we will calculate ($C = (A\vert b)$)
    
    /* copy |A| and |b| into |C| and L becomes the identity matrix */
    integer_matrix L(rows); // zero initialized
    for(i=0; i<rows; i++) { 
      for(j=0; j<cols; j++) 
        C(i,j)=A(i,j); 
      C(i,cols)=b[i]; 
      L(i,i) = 1; // diagonal elements are 1
    }


  
    array<int> var(0,cols - 1); 
    // column $j$ of |C| represents the |var[j]| - th variable
    // the array is indexed between |0| and |cols - 1|

    for(j=0; j<cols; j++)
      var[j]= j; // at the beginning, variable $x_j$ stands in column $j$

    RTINT denom = 1; // the determinant of an empty matrix is 1
    int sign = 1; // no interchanges yet
    int rank = 0; // we have not seen any non-zero row yet

    /* here comes the main loop */
    for(k=0; k<rows; k++)
    {   
      
         bool non_zero_found = false; 
         for(i =k; i<rows; i++)  // step through rows $k$ to |rows - 1|
         {    
           for (j = k ; j < cols && C(i,j) == 0; j++) ; 
           // step through columns |k| to |cols - 1|
           if (j < cols) 
           {  
             non_zero_found = true; 
             break; 
           }
         }


      if (non_zero_found)
      { 
        rank++; //increase the rank
        
          if (i != k) 
          { 
            sign = - sign; 
            /* we interchange rows |k| and |i| of |L| and |C| */
            swap(L[i],L[k]); 
            swap(C[i],C[k]); 
          }

          if (j != k)
          { 
            sign = - sign; 
            /* We interchange columns |k| and |j| */
            for(int ih = 0; ih < rows; ih++) 
              swap(C(ih,k),C(ih,j)); 

            /* We store the exchange of variables in |var| */
            swap(var[k],var[j]); 
          }
         
          
          for(i = k + 1; i < rows; i++)  
            for (j = 0; j <  rows; j++)  //and all columns of |L|
              L(i,j) = (L(i,j)*C(k,k) - C(i,k)*L(k,j))/denom; 

          for(i = k + 1; i < rows; i++) 
          {      
            /* the following iteration uses and changes |C(i,k)| */
            RTINT temp = C(i,k); 
            for (j = k; j <= cols; j++)
              C(i,j) = (C(i,j)*C(k,k) - temp*C(k,j))/denom; 
          }
          denom = C(k,k); 
          
          #ifdef LA_SELFTEST

            for(i = 0; i < rows; i++)
            { 
              for (j = 0; j < cols; j++)
              { 
                RTINT Sum = 0; 
                for (int l = 0; l < rows; l++) 
                  Sum += L(i,l)*A(l, var[j]); 
                if (Sum  != C(i,j)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 
              }
              RTINT Sum = 0; 
              for (int l = 0; l < rows; l++) 
                Sum += L(i,l)*b[l];
                if (Sum  != C(i,cols)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 

            }
          #endif


      }
      else 
        break; 
    }



  if (rank < rows) 
    return 0; 
  else
    return RTINT(sign) * denom; 
}


RTINT  
determinant(const integer_matrix& A, 
            integer_matrix& Ld, 
            integer_matrix& Ud, 
            array<int>& q, 
            integer_vector& c) 
{ 
  if (A.dim1() != A.dim2())
    error_handler(1,"determinant: only square matrices are legal inputs."); 
  integer_vector b(A.dim1()); // zero - vector
  
    int i,j,k; // indices to step through the matrix
    int rows = A.dim1(); 
    int cols = A.dim2(); 

    /* at this point one might want to check whether the computation can
       be carried out with doubles, see section \ref{ optimization }.
    */

    if (b.dim() != rows) 
      error_handler(1,"linear_solver: b has wrong dimension"); 

    integer_matrix C(rows,cols + 1); 
    // the matrix in which we will calculate ($C = (A\vert b)$)
    
    /* copy |A| and |b| into |C| and L becomes the identity matrix */
    integer_matrix L(rows); // zero initialized
    for(i=0; i<rows; i++) { 
      for(j=0; j<cols; j++) 
        C(i,j)=A(i,j); 
      C(i,cols)=b[i]; 
      L(i,i) = 1; // diagonal elements are 1
    }


  
    array<int> var(0,cols - 1); 
    // column $j$ of |C| represents the |var[j]| - th variable
    // the array is indexed between |0| and |cols - 1|

    for(j=0; j<cols; j++)
      var[j]= j; // at the beginning, variable $x_j$ stands in column $j$

    RTINT denom = 1; // the determinant of an empty matrix is 1
    int sign = 1; // no interchanges yet
    int rank = 0; // we have not seen any non-zero row yet

    /* here comes the main loop */
    for(k=0; k<rows; k++)
    {   
      
         bool non_zero_found = false; 
         for(i =k; i<rows; i++)  // step through rows $k$ to |rows - 1|
         {    
           for (j = k ; j < cols && C(i,j) == 0; j++) ; 
           // step through columns |k| to |cols - 1|
           if (j < cols) 
           {  
             non_zero_found = true; 
             break; 
           }
         }


      if (non_zero_found)
      { 
        rank++; //increase the rank
        
          if (i != k) 
          { 
            sign = - sign; 
            /* we interchange rows |k| and |i| of |L| and |C| */
            swap(L[i],L[k]); 
            swap(C[i],C[k]); 
          }

          if (j != k)
          { 
            sign = - sign; 
            /* We interchange columns |k| and |j| */
            for(int ih = 0; ih < rows; ih++) 
              swap(C(ih,k),C(ih,j)); 

            /* We store the exchange of variables in |var| */
            swap(var[k],var[j]); 
          }
         
          
          for(i = k + 1; i < rows; i++)  
            for (j = 0; j <  rows; j++)  //and all columns of |L|
              L(i,j) = (L(i,j)*C(k,k) - C(i,k)*L(k,j))/denom; 

          for(i = k + 1; i < rows; i++) 
          {      
            /* the following iteration uses and changes |C(i,k)| */
            RTINT temp = C(i,k); 
            for (j = k; j <= cols; j++)
              C(i,j) = (C(i,j)*C(k,k) - temp*C(k,j))/denom; 
          }
          denom = C(k,k); 
          
          #ifdef LA_SELFTEST

            for(i = 0; i < rows; i++)
            { 
              for (j = 0; j < cols; j++)
              { 
                RTINT Sum = 0; 
                for (int l = 0; l < rows; l++) 
                  Sum += L(i,l)*A(l, var[j]); 
                if (Sum  != C(i,j)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 
              }
              RTINT Sum = 0; 
              for (int l = 0; l < rows; l++) 
                Sum += L(i,l)*b[l];
                if (Sum  != C(i,cols)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 

            }
          #endif


      }
      else 
        break; 
    }



  if (rank < rows)
  { 
    c = L.row(rows - 1); 
    return 0; 
  }
  else
  { 
    Ld = L; 
    Ud = integer_matrix(rows); // quadratic
    for (i = 0; i < rows; i++) 
      for(j = 0; j <rows; j++) 
        Ud(i,j) = C(i,j); 
    q = var; 
    return RTINT(sign) * denom; 
  }
}


int   
sign_of_determinant(const integer_matrix& M)
{ 
  return sign(determinant(M)); 
}


bool  
verify_determinant(const integer_matrix& A, RTINT D, 
                   integer_matrix& L, 
                   integer_matrix& U, array<int> q, 
                   integer_vector& c) 
{ 
  if (q.low() != 0 || q.high() != A.dim2() - 1)
    error_handler(1,"verify_determinant: q should be a permutation array \
with index range [0,A.dim2() - 1]."); 
  int n = A.dim1(); 
  int i,j; 
  if (D == 0)
  { 
    /* we have $c^T \cdot A = 0$  */
    integer_vector zero(n); 
    return  (transpose(A) * c == zero); 
  }
  else
  { 
    /* we check the conditions on |L| and |U| */
    if (L(0,0) != 1) return false; 
    for (i = 0; i<n; i++)
    { 
      for (j = 0; j < i; j++) 
        if (U(i,j) != 0) 
          return false; 

      if (i > 0 && L(i,i) != U(i - 1,i - 1)) 
        return false; 

      for (j = i + 1; j < n; j++) 
        if (L(i,j) != 0) 
          return false; 
    }

    /* check whether $L \cdot A \cdot Q = U$ */
   integer_matrix LA = L * A; 
   for (j = 0; j < n; j++)
     if (LA.col(q[j]) != U.col(j)) 
       return false; 

    /* compute sign |s| of |Q| */
    int sign = 1; 

    /* we chase the cycles of |q|. An even length cycle contributes - 1
       and vice versa */

    array<bool> already_considered(0,n - 1); 

    for (i = 0; i < n; i++)  
      already_considered[i] = false; 

    for (i = 0; i < n; i++)  
      already_considered[q[i]] = true; 

    for (i = 0; i < n; i++) 
      if (! already_considered[i])
        error_handler(1,"verify_determinant:q is not a permutation."); 
      else 
        already_considered[i] = false; 

    for (i = 0; i < n; i++)
    { 
      if (already_considered[i]) 
        continue; 

      /* we have found a new cycle with minimal element $i$. */

      int k = q[i]; 
      already_considered[i] =true; 

      while (k != i)
      { 
        sign = - sign; 
        already_considered[k]= true; 
        k = q[k]; 
      }
    }

    return (D == RTINT(sign) * U(n - 1,n - 1)); 

  }
}


void  
independent_columns(const integer_matrix& A, array<int>& columns) 
{ 
  integer_vector b(A.dim1()); // zero - vector
  
    int i,j,k; // indices to step through the matrix
    int rows = A.dim1(); 
    int cols = A.dim2(); 

    /* at this point one might want to check whether the computation can
       be carried out with doubles, see section \ref{ optimization }.
    */

    if (b.dim() != rows) 
      error_handler(1,"linear_solver: b has wrong dimension"); 

    integer_matrix C(rows,cols + 1); 
    // the matrix in which we will calculate ($C = (A\vert b)$)
    
    /* copy |A| and |b| into |C| and L becomes the identity matrix */
    integer_matrix L(rows); // zero initialized
    for(i=0; i<rows; i++) { 
      for(j=0; j<cols; j++) 
        C(i,j)=A(i,j); 
      C(i,cols)=b[i]; 
      L(i,i) = 1; // diagonal elements are 1
    }


  
    array<int> var(0,cols - 1); 
    // column $j$ of |C| represents the |var[j]| - th variable
    // the array is indexed between |0| and |cols - 1|

    for(j=0; j<cols; j++)
      var[j]= j; // at the beginning, variable $x_j$ stands in column $j$

    RTINT denom = 1; // the determinant of an empty matrix is 1
    int sign = 1; // no interchanges yet
    int rank = 0; // we have not seen any non-zero row yet

    /* here comes the main loop */
    for(k=0; k<rows; k++)
    {   
      
         bool non_zero_found = false; 
         for(i =k; i<rows; i++)  // step through rows $k$ to |rows - 1|
         {    
           for (j = k ; j < cols && C(i,j) == 0; j++) ; 
           // step through columns |k| to |cols - 1|
           if (j < cols) 
           {  
             non_zero_found = true; 
             break; 
           }
         }


      if (non_zero_found)
      { 
        rank++; //increase the rank
        
          if (i != k) 
          { 
            sign = - sign; 
            /* we interchange rows |k| and |i| of |L| and |C| */
            swap(L[i],L[k]); 
            swap(C[i],C[k]); 
          }

          if (j != k)
          { 
            sign = - sign; 
            /* We interchange columns |k| and |j| */
            for(int ih = 0; ih < rows; ih++) 
              swap(C(ih,k),C(ih,j)); 

            /* We store the exchange of variables in |var| */
            swap(var[k],var[j]); 
          }
         
          
          for(i = k + 1; i < rows; i++)  
            for (j = 0; j <  rows; j++)  //and all columns of |L|
              L(i,j) = (L(i,j)*C(k,k) - C(i,k)*L(k,j))/denom; 

          for(i = k + 1; i < rows; i++) 
          {      
            /* the following iteration uses and changes |C(i,k)| */
            RTINT temp = C(i,k); 
            for (j = k; j <= cols; j++)
              C(i,j) = (C(i,j)*C(k,k) - temp*C(k,j))/denom; 
          }
          denom = C(k,k); 
          
          #ifdef LA_SELFTEST

            for(i = 0; i < rows; i++)
            { 
              for (j = 0; j < cols; j++)
              { 
                RTINT Sum = 0; 
                for (int l = 0; l < rows; l++) 
                  Sum += L(i,l)*A(l, var[j]); 
                if (Sum  != C(i,j)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 
              }
              RTINT Sum = 0; 
              for (int l = 0; l < rows; l++) 
                Sum += L(i,l)*b[l];
                if (Sum  != C(i,cols)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 

            }
          #endif


      }
      else 
        break; 
    }


  /* at this point we have:
     |C| has an $rank \times rank$ upper triangular matrix 
     in its left upper corner; 
     |var| tells us the columns of |A| corresponding to the
     dependent variables; */
     
  columns = array<int>(rank); 
  for(i = 0; i < rank; i++) 
    columns[i] = var[i]; 
}


int  
rank(const integer_matrix& A)
{ 
  integer_vector b(A.dim1()); // zero - vector
  
    int i,j,k; // indices to step through the matrix
    int rows = A.dim1(); 
    int cols = A.dim2(); 

    /* at this point one might want to check whether the computation can
       be carried out with doubles, see section \ref{ optimization }.
    */

    if (b.dim() != rows) 
      error_handler(1,"linear_solver: b has wrong dimension"); 

    integer_matrix C(rows,cols + 1); 
    // the matrix in which we will calculate ($C = (A\vert b)$)
    
    /* copy |A| and |b| into |C| and L becomes the identity matrix */
    integer_matrix L(rows); // zero initialized
    for(i=0; i<rows; i++) { 
      for(j=0; j<cols; j++) 
        C(i,j)=A(i,j); 
      C(i,cols)=b[i]; 
      L(i,i) = 1; // diagonal elements are 1
    }


  
    array<int> var(0,cols - 1); 
    // column $j$ of |C| represents the |var[j]| - th variable
    // the array is indexed between |0| and |cols - 1|

    for(j=0; j<cols; j++)
      var[j]= j; // at the beginning, variable $x_j$ stands in column $j$

    RTINT denom = 1; // the determinant of an empty matrix is 1
    int sign = 1; // no interchanges yet
    int rank = 0; // we have not seen any non-zero row yet

    /* here comes the main loop */
    for(k=0; k<rows; k++)
    {   
      
         bool non_zero_found = false; 
         for(i =k; i<rows; i++)  // step through rows $k$ to |rows - 1|
         {    
           for (j = k ; j < cols && C(i,j) == 0; j++) ; 
           // step through columns |k| to |cols - 1|
           if (j < cols) 
           {  
             non_zero_found = true; 
             break; 
           }
         }


      if (non_zero_found)
      { 
        rank++; //increase the rank
        
          if (i != k) 
          { 
            sign = - sign; 
            /* we interchange rows |k| and |i| of |L| and |C| */
            swap(L[i],L[k]); 
            swap(C[i],C[k]); 
          }

          if (j != k)
          { 
            sign = - sign; 
            /* We interchange columns |k| and |j| */
            for(int ih = 0; ih < rows; ih++) 
              swap(C(ih,k),C(ih,j)); 

            /* We store the exchange of variables in |var| */
            swap(var[k],var[j]); 
          }
         
          
          for(i = k + 1; i < rows; i++)  
            for (j = 0; j <  rows; j++)  //and all columns of |L|
              L(i,j) = (L(i,j)*C(k,k) - C(i,k)*L(k,j))/denom; 

          for(i = k + 1; i < rows; i++) 
          {      
            /* the following iteration uses and changes |C(i,k)| */
            RTINT temp = C(i,k); 
            for (j = k; j <= cols; j++)
              C(i,j) = (C(i,j)*C(k,k) - temp*C(k,j))/denom; 
          }
          denom = C(k,k); 
          
          #ifdef LA_SELFTEST

            for(i = 0; i < rows; i++)
            { 
              for (j = 0; j < cols; j++)
              { 
                RTINT Sum = 0; 
                for (int l = 0; l < rows; l++) 
                  Sum += L(i,l)*A(l, var[j]); 
                if (Sum  != C(i,j)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 
              }
              RTINT Sum = 0; 
              for (int l = 0; l < rows; l++) 
                Sum += L(i,l)*b[l];
                if (Sum  != C(i,cols)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 

            }
          #endif


      }
      else 
        break; 
    }


  return rank; 
}

bool  
inverse(const integer_matrix& A, integer_matrix& inverse, RTINT& D, integer_vector& c)
{ 
  if (A.dim1() != A.dim2())
    error_handler(1,"inverse: only square matrices are legal inputs."); 
  integer_vector b(A.dim1()); // zero - vector
  
    int i,j,k; // indices to step through the matrix
    int rows = A.dim1(); 
    int cols = A.dim2(); 

    /* at this point one might want to check whether the computation can
       be carried out with doubles, see section \ref{ optimization }.
    */

    if (b.dim() != rows) 
      error_handler(1,"linear_solver: b has wrong dimension"); 

    integer_matrix C(rows,cols + 1); 
    // the matrix in which we will calculate ($C = (A\vert b)$)
    
    /* copy |A| and |b| into |C| and L becomes the identity matrix */
    integer_matrix L(rows); // zero initialized
    for(i=0; i<rows; i++) { 
      for(j=0; j<cols; j++) 
        C(i,j)=A(i,j); 
      C(i,cols)=b[i]; 
      L(i,i) = 1; // diagonal elements are 1
    }


  
    array<int> var(0,cols - 1); 
    // column $j$ of |C| represents the |var[j]| - th variable
    // the array is indexed between |0| and |cols - 1|

    for(j=0; j<cols; j++)
      var[j]= j; // at the beginning, variable $x_j$ stands in column $j$

    RTINT denom = 1; // the determinant of an empty matrix is 1
    int sign = 1; // no interchanges yet
    int rank = 0; // we have not seen any non-zero row yet

    /* here comes the main loop */
    for(k=0; k<rows; k++)
    {   
      
         bool non_zero_found = false; 
         for(i =k; i<rows; i++)  // step through rows $k$ to |rows - 1|
         {    
           for (j = k ; j < cols && C(i,j) == 0; j++) ; 
           // step through columns |k| to |cols - 1|
           if (j < cols) 
           {  
             non_zero_found = true; 
             break; 
           }
         }


      if (non_zero_found)
      { 
        rank++; //increase the rank
        
          if (i != k) 
          { 
            sign = - sign; 
            /* we interchange rows |k| and |i| of |L| and |C| */
            swap(L[i],L[k]); 
            swap(C[i],C[k]); 
          }

          if (j != k)
          { 
            sign = - sign; 
            /* We interchange columns |k| and |j| */
            for(int ih = 0; ih < rows; ih++) 
              swap(C(ih,k),C(ih,j)); 

            /* We store the exchange of variables in |var| */
            swap(var[k],var[j]); 
          }
         
          
          for(i = k + 1; i < rows; i++)  
            for (j = 0; j <  rows; j++)  //and all columns of |L|
              L(i,j) = (L(i,j)*C(k,k) - C(i,k)*L(k,j))/denom; 

          for(i = k + 1; i < rows; i++) 
          {      
            /* the following iteration uses and changes |C(i,k)| */
            RTINT temp = C(i,k); 
            for (j = k; j <= cols; j++)
              C(i,j) = (C(i,j)*C(k,k) - temp*C(k,j))/denom; 
          }
          denom = C(k,k); 
          
          #ifdef LA_SELFTEST

            for(i = 0; i < rows; i++)
            { 
              for (j = 0; j < cols; j++)
              { 
                RTINT Sum = 0; 
                for (int l = 0; l < rows; l++) 
                  Sum += L(i,l)*A(l, var[j]); 
                if (Sum  != C(i,j)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 
              }
              RTINT Sum = 0; 
              for (int l = 0; l < rows; l++) 
                Sum += L(i,l)*b[l];
                if (Sum  != C(i,cols)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 

            }
          #endif


      }
      else 
        break; 
    }


  
  if (rank < rows)
  { 
    /* matrix is singular; 
       we return a vector $c$ with $c^T \cdot A = 0$.*/

    c = integer_vector(rows); 
    for (j = 0; j <rows; j++) 
      c[j] = L(rows - 1,j); 

    return false; 
  }
  
    D = denom; 
    inverse = integer_matrix(rows); //quadratic
    RTINT h; 
    for(i = 0; i <rows; i++) 
    {  // $i$-th column of inverse
      for (j = rows - 1; j >= 0; j--) 
      { 
        h = L (j,i) * D; 
        for (int l = j + 1; l<rows; l++)  
          h -= C(j,l)*inverse(var[l],i); 
        inverse(var[j],i) = h / C(j,j); 
      }
    }  

  #ifdef LA_SELFTEST
    if (!(A*inverse == integer_matrix::identity(rows)*D))
      error_handler(1,"inverse: matrix inverse computed incorrectly."); 
  #endif      

  return true; 
}


int  
homogeneous_linear_solver(const integer_matrix &A, 
                          integer_matrix& spanning_vectors)
/* returns the dimension of the solution space of the homogeneous system 
   $Ax = 0$. The columns of spanning\_vectors span the solution space. */
{ 
  integer_vector b(A.dim1()); // zero - vector
  RTINT D; 
  
    int i,j,k; // indices to step through the matrix
    int rows = A.dim1(); 
    int cols = A.dim2(); 

    /* at this point one might want to check whether the computation can
       be carried out with doubles, see section \ref{ optimization }.
    */

    if (b.dim() != rows) 
      error_handler(1,"linear_solver: b has wrong dimension"); 

    integer_matrix C(rows,cols + 1); 
    // the matrix in which we will calculate ($C = (A\vert b)$)
    
    /* copy |A| and |b| into |C| and L becomes the identity matrix */
    integer_matrix L(rows); // zero initialized
    for(i=0; i<rows; i++) { 
      for(j=0; j<cols; j++) 
        C(i,j)=A(i,j); 
      C(i,cols)=b[i]; 
      L(i,i) = 1; // diagonal elements are 1
    }


  
    array<int> var(0,cols - 1); 
    // column $j$ of |C| represents the |var[j]| - th variable
    // the array is indexed between |0| and |cols - 1|

    for(j=0; j<cols; j++)
      var[j]= j; // at the beginning, variable $x_j$ stands in column $j$

    RTINT denom = 1; // the determinant of an empty matrix is 1
    int sign = 1; // no interchanges yet
    int rank = 0; // we have not seen any non-zero row yet

    /* here comes the main loop */
    for(k=0; k<rows; k++)
    {   
      
         bool non_zero_found = false; 
         for(i =k; i<rows; i++)  // step through rows $k$ to |rows - 1|
         {    
           for (j = k ; j < cols && C(i,j) == 0; j++) ; 
           // step through columns |k| to |cols - 1|
           if (j < cols) 
           {  
             non_zero_found = true; 
             break; 
           }
         }


      if (non_zero_found)
      { 
        rank++; //increase the rank
        
          if (i != k) 
          { 
            sign = - sign; 
            /* we interchange rows |k| and |i| of |L| and |C| */
            swap(L[i],L[k]); 
            swap(C[i],C[k]); 
          }

          if (j != k)
          { 
            sign = - sign; 
            /* We interchange columns |k| and |j| */
            for(int ih = 0; ih < rows; ih++) 
              swap(C(ih,k),C(ih,j)); 

            /* We store the exchange of variables in |var| */
            swap(var[k],var[j]); 
          }
         
          
          for(i = k + 1; i < rows; i++)  
            for (j = 0; j <  rows; j++)  //and all columns of |L|
              L(i,j) = (L(i,j)*C(k,k) - C(i,k)*L(k,j))/denom; 

          for(i = k + 1; i < rows; i++) 
          {      
            /* the following iteration uses and changes |C(i,k)| */
            RTINT temp = C(i,k); 
            for (j = k; j <= cols; j++)
              C(i,j) = (C(i,j)*C(k,k) - temp*C(k,j))/denom; 
          }
          denom = C(k,k); 
          
          #ifdef LA_SELFTEST

            for(i = 0; i < rows; i++)
            { 
              for (j = 0; j < cols; j++)
              { 
                RTINT Sum = 0; 
                for (int l = 0; l < rows; l++) 
                  Sum += L(i,l)*A(l, var[j]); 
                if (Sum  != C(i,j)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 
              }
              RTINT Sum = 0; 
              for (int l = 0; l < rows; l++) 
                Sum += L(i,l)*b[l];
                if (Sum  != C(i,cols)) 
                  error_handler(1,"linear_solver: L*A*P different from C"); 

            }
          #endif


      }
      else 
        break; 
    }


  integer_vector x; 
  
    x = integer_vector(cols); 
    D = denom; 
    for(i=rank - 1 ; i>= 0; i--) 
    { 
      RTINT h =C(i,cols) * D; 
      for (j = i + 1; j<rank; j++)  
      { 
        h -= C(i,j)*x[var[j]]; 
      }
      x[var[i]]= h / C(i,i); 
    }

  #ifdef LA_SELFTEST
    /* we check whether |x| is a solution */
    { 
      for (i = 0; i < rows; i++)
      { 
        RTINT sum = 0; 
        for (j =0; j < cols; j++) 
          sum += A(i,j)*x[j]; 
        if (sum != D * b[i])
          error_handler(1,"linear_solver: base is not a solution"); 
      }
    }
  #endif

    int dimension = cols - rank; //dimension of solution
    spanning_vectors = integer_matrix(cols,dimension); 
   
    if (dimension > 0)
    { 
   
  /* In the $l$ - th spanning vector, $0 \le l < |dimension|$ we set
  variable |var[rank + l]| to $1 = |denom|/|denom|$ and then the
  dependent variables as dictated by the $|rank| + l$ - th column of
  |C|.*/

      for(int l=0; l < dimension; l++) 
      { 
        spanning_vectors(var[rank + l],l)=D; 
        for(i = rank - 1; i >= 0 ; i--)
        { 
          RTINT h = - C(i,rank + l)* D; 
          for ( j= i + 1; j<rank; j++)  
            h -= C(i,j)*spanning_vectors(var[j],l); 
          spanning_vectors(var[i],l)=h / C(i,i); 
        }

  #ifdef LA_SELFTEST
        /* we check whether the $l$ - th spanning vector is a solution 
           of the homogeneous system */
        { 
          integer_vector zero(rows); 
          if (A *spanning_vectors.col(l) != zero)
            error_handler(1,"linear_solver: spanning_vector is not a solution."); 
        }
  #endif
      }
    }
; 
  return dimension; 
}

bool  
homogeneous_linear_solver(const integer_matrix& A, integer_vector& x)
/* returns true if the homogeneous system $Ax = 0$ has a non - trivial
   solution and false otherwise. */
{  
  integer_matrix spanning_vectors; 
  int dimension = homogeneous_linear_solver(A,spanning_vectors); 

  if (dimension == 0) 
    return false; 

  /* return first column of |spanning_vectors| */
  for (int i = 0; i < spanning_vectors.dim1() ; i++)  
    x[i] = spanning_vectors(i,0); 
  return true; 
}



integer_matrix  
transpose(const integer_matrix& M)
{ 
  int d1 = M.dim1(); 
  int d2 = M.dim2(); 
  integer_matrix result(d2,d1); 
  for(int i = 0; i < d2; i++)
    for(int j = 0; j < d1; j++)
      result.elem(i,j) = M(j,i); 
  return result; 
}

 





