// $Id: matrix.lib,v 1.9 1998-09-30 17:24:32 Singular Exp $
// GMG/BM, last modified: 30.9.98
///////////////////////////////////////////////////////////////////////////////

version="$Id: matrix.lib,v 1.9 1998-09-30 17:24:32 Singular Exp $";
info="
LIBRARY:  matrix.lib    PROCEDURES FOR MATRIX OPERATIONS
                        Authors: GMG/BM

 compress(A);           matrix, zero columns from A deleted
 concat(A1,A2,..);      matrix, concatenation of matrices A1,A2,...
 diag(p,n);             matrix, nxn diagonal matrix with entries poly p
 dsum(A1,A2,..);        matrix, direct sum of matrices A1,A2,...
 flatten(A);            ideal, generated by entries of matrix A
 genericmat(n,m[,id]);  generic nxm matrix [entries from id]
 is_complex(c);         1 if list c is a complex, 0 if not
 outer(A,B);            matrix, outer product of matrices A and B
 power(A,n);            matrix/intmat, n-th power of matrix/intmat A
 skewmat(n[,id]);       generic skew-symmetric nxn matrix [entries from id]
 submat(A,r,c);         submatrix of A with rows/cols specified by intvec r/c
 symmat(n[,id]);        generic symmetric nxn matrix [entries from id]
 tensor(A,B);           matrix, tensor product of matrices A nd B
 unitmat(n);            unit square matrix of size n
 gauss_col(A);          transform constant matrix A into col-reduced nf
 gauss_row(A);          transform constant matrix A into row-reduced nf
 addcol(A,c1,p,c2);     add p*(c1-th col) to c2-th column of matrix A, p poly
 addrow(A,r1,p,r2);     add p*(r1-th row) to r2-th row of matrix A, p poly
 multcol(A,c,p);        multiply c-th column of A with poly p
 multrow(A,r,p);        multiply r-th row of A with poly p
 permcol(A,i,j);        permute i-th and j-th columns
 permrow(A,i,j);        permute i-th and j-th rows
 rowred(A[,any]);       reduction of matrix A with elementary row-operations
 colred(A[,any]);       reduction of matrix A with elementary col-operations
 rm_unitrow(A);         remove unit rows and associated columns of A
 rm_unitcol(A);         remove unit columns and associated rows of A
          (parameters in square brackets [] are optional)
";

LIB "inout.lib";
LIB "ring.lib";
LIB "random.lib";
///////////////////////////////////////////////////////////////////////////////

proc compress (A)
"USAGE:   compress(A); A matrix/ideal/module/intmat/intvec
RETURN:  same type, zero columns/generators from A deleted
         (if A=intvec, zero elements are deleted)
EXAMPLE: example compress; shows an example
"
{
   if( typeof(A)=="matrix" ) { return(matrix(simplify(A,2))); }
   if( typeof(A)=="intmat" or typeof(A)=="intvec" )
   {
      ring r=0,x,lp;
      if( typeof(A)=="intvec" ) { intmat C=transpose(A); kill A; intmat A=C; }
      module m = matrix(A);
      intmat B[nrows(A)][size(m)];
      int i,j;
      for( i=1; i<=ncols(A); i=i+1 )
      {
         if( m[i]!=[0] )
         {
            j=j+1;
            B[1..nrows(A),j]=A[1..nrows(A),i];
         }
      }
      if( defined(C) ) { return(intvec(B)); }
      return(B);
    }
   return(simplify(A,2));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),ds;
   matrix A[3][4]=1,0,3,0,x,0,z,0,x2,0,z2,0;
   print(A);
   print(compress(A));
   module m=module(A); show(m);
   show(compress(m));
   intmat B[3][4]=1,0,3,0,4,0,5,0,6,0,7,0;
   compress(B);
   intvec C=0,0,1,2,0,3;
   compress(C);
}
///////////////////////////////////////////////////////////////////////////////
proc concat (list #)
"USAGE:   concat(A1,A2,..); A1,A2,... matrices
RETURN:  matrix, concatenation of A1,A2,... . Number of rows of result matrix
         is max(nrows(A1),nrows(A2),...)
EXAMPLE: example concat; shows an example
"
{
   int i;
   module B=module(#[1]);
   for( i=2; i<=size(#); i=i+1 ) { B=B,module(#[i]); }
   return(matrix(B));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),ds;
   matrix A[3][3]=1,2,3,x,y,z,x2,y2,z2;
   matrix B[2][2]=1,0,2,0; matrix C[1][4]=4,5,x,y;
   print(A);
   print(B);
   print(C);
   print(concat(A,B,C));
}
///////////////////////////////////////////////////////////////////////////////

proc diag (list #)
"USAGE:   diag(p,n); p poly, n integer
         diag(A);   A matrix
RETURN:  diag(p,n): diagonal matrix, p times unitmatrix of size n
         diag(A)  : n*m x n*m diagonal matrix with entries all the entries of
                   the nxm matrix A, taken from the 1st row, 2nd row etc of A
EXAMPLE: example diag; shows an example
"
{
   if( size(#)==2 ) { return(matrix(#[1]*freemodule(#[2]))); }
   if( size(#)==1 )
   {
      int i; ideal id=#[1];
      int n=ncols(id); matrix A[n][n];
      for( i=1; i<=n; i=i+1 ) { A[i,i]=id[i]; }
   }
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),ds;
   print(diag(xy,4));
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   print(A);
   print(diag(A));
}
///////////////////////////////////////////////////////////////////////////////

proc dsum (list #)
"USAGE:   dsum(A1,A2,..); A1,A2,... matrices
RETURN:  matrix, direct sum of A1,A2,...
EXAMPLE: example dsum; shows an example
"
{
   int i,N,a;
   list L;
   for( i=1; i<=size(#); i=i+1 ) { N=N+nrows(#[i]); }
   for( i=1; i<=size(#); i=i+1 )
   {
      matrix B[N][ncols(#[i])];
      B[a+1..a+nrows(#[i]),1..ncols(#[i])]=#[i];
      a=a+nrows(#[i]);
      L[i]=B; kill B;
   }
   return(concat(L));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),ds;
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   matrix B[2][2]=1,x,y,z;
   matrix C[1][4]=4,5,x,y;
   print(A);
   print(B);
   print(C);
   print(dsum(A,B,C));
}
///////////////////////////////////////////////////////////////////////////////

proc flatten (matrix A)
"USAGE:   flatten(A); A matrix
RETURN:  ideal, generated by all entries from A
EXAMPLE: example flatten; shows an example
"
{
   return(ideal(A));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),ds;
   matrix A[3][3]=1,2,3,x,y,z,7,8,9;
   print(A);
   flatten(A);
}
///////////////////////////////////////////////////////////////////////////////

proc genericmat (int n,int m,list #)
"USAGE:   genericmat(n,m[,id]);  n,m=integers, id=ideal
RETURN:  nxm matrix, with entries from id (default: id=maxideal(1))
NOTE:    if id has less than nxm elements, the matrix is filled with 0's,
         genericmat(n,m); creates the generic nxm matrix
EXAMPLE: example genericmat; shows an example
"
{
   if( size(#)==0 ) { ideal id=maxideal(1); }
   if( size(#)==1 ) { ideal id=#[1]; }
   if( size(#)>=2 ) { "// give 3 arguments, 3-rd argument must be an ideal"; }
   matrix B[n][m]=id;
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R=0,x(1..16),lp;
   print(genericmat(4,4));    // the generic 4x4 matrix
   changevar("R1",A_Z("a",4),R);
   matrix A=genericmat(4,5,maxideal(1)^3);
   print(A);
   int n,m=4,3;
   ideal i = ideal(randommat(1,n*m,maxideal(1),9));
   print(genericmat(n,m,i));  // matrix of generic linear forms
   kill R1;
}
///////////////////////////////////////////////////////////////////////////////

proc is_complex (list c)
"USAGE:   is_complex(c); c = list of size-compatible modules or matrices
RETURN:  1 if c[i]*c[i+1]=0 for all i, 0 if not, hence checking whether the
         list of matrices forms a complex
NOTE:    Ideals are treated internally as 1-line matrices
EXAMPLE: example is_complex; shows an example
"
{
   int i;
   module @test;
   for( i=1; i<=size(c)-1; i=i+1 )
   {
      c[i]=matrix(c[i]); c[i+1]=matrix(c[i+1]);
      @test=c[i]*c[i+1];
      if (size(@test)!=0)
      {
         if( voice==2 ) { "// argument is not a complex at position",i; }
         return(0);
      }
   }
   if( voice==2 ) { "// argument is a complex"; }
   return(1);
}
example
{ "EXAMPLE:";   echo = 2;
   ring r=32003,(x,y,z),ds;
   ideal i=x4+y5+z6,xyz,yx2+xz2+zy7;
   list L=nres(i,0);
   is_complex(L);
   L[4]=matrix(i);
   is_complex(L);
}
///////////////////////////////////////////////////////////////////////////////

proc outer (matrix A, matrix B)
"USAGE:   outer(A,B); A,B matrices
RETURN:  matrix, outer product of A and B
EXAMPLE: example outer; shows an example
"
{
   int i,j; list L;
   int triv = nrows(B)*ncols(B);
   if( triv==1 )
   {
     return(B[1,1]*A);
   }
   else
   {
     int N = nrows(A)*nrows(B);
     matrix C[N][ncols(B)];
     for( i=1; i<=ncols(A); i=i+1 )
     {
       for( j=1; j<=nrows(A); j=j+1 )
       {
          C[(j-1)*nrows(B)+1..j*nrows(B),1..ncols(B)]=A[j,i]*B;
       }
       L[i]=C;
     }
     return(concat(L));
   }
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),ds;
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   matrix B[2][2]=x,y,0,z;
   print(A);
   print(B);
   print(outer(A,B));
}
///////////////////////////////////////////////////////////////////////////////

proc power ( A, int n)
"USAGE:   power(A,n);  A a square-matrix of type intmat or matrix, n=integer
RETURN:  inmat resp. matrix, the n-th power of A
NOTE:    for A=intamt and big n the result may be wrong because of int overflow
EXAMPLE: example power; shows an example
"
{
//---------------------------- type checking ----------------------------------
   if( typeof(A)!="matrix" and typeof(A)!="intmat" )
   {
      "// no matrix or intmat!";
      return (A);
   }
   if( ncols(A) != nrows(A) )
   {
      "// not a suare matrix!";
      return();
   }
//---------------------------- trivial cases ----------------------------------
   int ii;
   if( n <= 0 )
   {
      if( typeof(A)=="matrix" )
      {
         return (unitmat(nrows(A)));
      }
      if( typeof(A)=="intmat" )
      {
         intmat B[nrows(A)][nrows(A)];
         for( ii=1; ii<=nrows(A); ii++ )
         {
            B[ii,ii] = 1;
         }
         return (B);
      }
   }
   if( n == 1 ) { return (A); }
//---------------------------- sub procedure ----------------------------------
   proc matpow (A, int n)
   {
      def B = A*A;
      int ii= 2;
      int jj= 4;
      while( jj <= n )
      {
         B=B*B;
         ii=jj;
         jj=2*jj;
      }
      return(B,n-ii);
   }
//----------------------------- main program ----------------------------------
   list L = matpow(A,n);
   def B  = L[1];
   ii     = L[2];
   while( ii>=2 )
   {
      L = matpow(A,ii);
      B = B*L[1];
      ii= L[2];
   }
   if( ii == 0) { return(B); }
   if( ii == 1) { return(A*B); }
}
example
{ "EXAMPLE:"; echo = 2;
   intmat A[3][3]=1,2,3,4,5,6,7,8,9;
   print(power(A,3));"";
   ring r=0,(x,y,z),dp;
   matrix B[4][4]=0,x,y,z,0,0,y,z,0,0,0,z,x,y,z,0;
   print(power(B,3));"";
   matrix C[3][3]=1,2,3,4,5,6,7,8,9;
   power(C,50);
}
///////////////////////////////////////////////////////////////////////////////

proc skewmat (int n, list #)
"USAGE:   skewmat(n[,id]);  n integer, id ideal
RETURN:  skew-symmetric nxn matrix, with entries from id
         (default: id=maxideal(1))
NOTE:    if id has less than n*(n-1)/2 elements, the matrix is filled with 0's,
         skewmat(n); creates the generic skew-symmetric matrix
EXAMPLE: example skewmat; shows an example
"
{
   matrix B[n][n];
   if( size(#)==0 ) { ideal id=maxideal(1); }
   else { ideal id=#[1]; }
   id = id,B[1..n,1..n];
   int i,j;
   for( i=0; i<=n-2; i=i+1 )
   {
      B[i+1,i+2..n]=id[j+1..j+n-i-1];
      j=j+n-i-1;
   }
   matrix A=transpose(B);
   B=B-A;
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R=0,x(1..5),lp;
   print(skewmat(4));    // the generic skew-symmetric matrix
   changevar("R1",A_Z("a",5),R);
   matrix A=skewmat(6,maxideal(1)^2);
   print(A);
   int n=4;
   ideal i = ideal(randommat(1,n*(n-1) div 2,maxideal(1),9));
   print(skewmat(n,i));  // skew matrix of generic linear forms
   kill R1;
}
///////////////////////////////////////////////////////////////////////////////

proc submat (matrix A, intvec r, intvec c)
"USAGE:   submat(A,r,c);  A=matrix, r,c=intvec
RETURN:  matrix, submatrix of A with rows specified by intvec r and columns
         specified by intvec c
EXAMPLE: example submat; shows an example
"
{
   matrix B[size(r)][size(c)]=A[r,c];
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R=32003,(x,y,z),lp;
   matrix A[4][4]=x,y,z,0,1,2,3,4,5,6,7,8,9,x2,y2,z2;
   print(A);
   intvec v=1,3,4;
   matrix B=submat(A,v,1..3);
   print(B);
}
///////////////////////////////////////////////////////////////////////////////

proc symmat (int n, list #)
"USAGE:   symmat(n[,id]);  n integer, id ideal
RETURN:  symmetric nxn matrix, with entries from id (default: id=maxideal(1))
NOTE:    if id has less than n*(n+1)/2 elements, the matrix is filled with 0's,
         symmat(n); creates the generic symmetric matrix
EXAMPLE: example symmat; shows an example
"
{
   matrix B[n][n];
   if( size(#)==0 ) { ideal id=maxideal(1); }
   else { ideal id=#[1]; }
   id = id,B[1..n,1..n];
   int i,j;
   for( i=0; i<=n-1; i=i+1 )
   {
      B[i+1,i+1..n]=id[j+1..j+n-i];
      j=j+n-i;
   }
   matrix A=transpose(B);
   for( i=1; i<=n; i=i+1 ) {  A[i,i]=0; }
   B=A+B;
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R=0,x(1..10),lp;
   print(symmat(4));    // the generic symmetric matrix
   changevar("R1",A_Z("a",5),R);
   matrix A=symmat(5,maxideal(1)^2);
   print(A);
   int n=3;
   ideal i = ideal(randommat(1,n*(n+1) div 2,maxideal(1),9));
   print(symmat(n,i));  // symmetric matrix of generic linear forms
   kill R1;
}
///////////////////////////////////////////////////////////////////////////////

proc tensor (matrix A, matrix B)
"USAGE:   tensor(A,B); A,B matrices
RETURN:  matrix, tensor product of A and B
EXAMPLE: example tensor; shows an example
"
{
   int i,j;
   matrix C=B;
   for( i=2; i<=nrows(A); i=i+1 ) { C=dsum(C,B); }
   matrix D[nrows(C)][ncols(A)*nrows(B)];
   for( j=1; j<=nrows(B); j=j+1 )
   {
      for( i=1; i<=nrows(A); i=i+1 )
      {
         D[(i-1)*nrows(B)+j,(j-1)*ncols(A)+1..j*ncols(A)]=A[i,1..ncols(A)];
      }
   }
   return(concat(C,D));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),(c,ds);
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   matrix B[2][2]=x,y,0,z;
   print(A);
   print(B);
   print(tensor(A,B));
}
///////////////////////////////////////////////////////////////////////////////

proc unitmat (int n)
"USAGE:   unitmat(n);  n integer >= 0
RETURN:  nxn unit matrix
NOTE:    needs a basering, diagonal entries are numbers (=1) in the basering
EXAMPLE: example unitmat; shows an example
"
{
   return(matrix(freemodule(n)));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   print(xyz*unitmat(4));
   print(unitmat(5));
}
///////////////////////////////////////////////////////////////////////////////

proc gauss_col (matrix A)
"USAGE:   gauss_col(A); A=matrix with constant coefficients
RETURN:  matrix = col-reduced lower-triagonal normal form of A
NOTE:    works fine for constant matrices, use proc colred for poly matrices
EXAMPLE: example gauss_col; shows an example
"
{
   def R=basering;
   intvec v = option(get);
   changeord("@R","ds,c",R);
   option(redSB); option(nointStrategy);
   matrix A = imap(R,A);
   A = matrix(std(A),nrows(A),ncols(A));
   setring R;
   A=imap(@R,A);
   option(set,v);
   kill @R;
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring S=0,x,dp;
   matrix A[5][4] =  3, 1,1,-1,
                    13, 8,6,-7,
                    14,10,6,-7,
                     7, 4,3,-3,
                     2, 1,0, 3;
   print(gauss_col(A));
}
///////////////////////////////////////////////////////////////////////////////

proc gauss_row (matrix A)
"USAGE:   gauss_row(A); A=matrix with constant coefficients
RETURN:  matrix = row-reduced upper-triangular normal form of A
NOTE:    may be used to solve a system of linear equations
         works fine for constant matrices, use proc rowred for poly matrices
EXAMPLE: example gauss_row; shows an example
"
{
   A = gauss_col(transpose(A));
   return(transpose(A));
}
example
{ "EXAMPLE:"; echo = 2;
   ring S=0,x,dp;
   matrix A[4][5] =  3, 1,1,-1,2,
                    13, 8,6,-7,1,
                    14,10,6,-7,1,
                     7, 4,3,-3,3;
   print(gauss_row(A));
}
///////////////////////////////////////////////////////////////////////////////

proc addcol (matrix A, int c1, poly p, int c2)
"USAGE:   addcol(A,c1,p,c2);  A matrix, p poly, c1, c2 positive integers
RETURN:  matrix,  A being modified by adding p times column c1 to column c2
EXAMPLE: example addcol; shows an example
"
{
   A[1..nrows(A),c2]=A[1..nrows(A),c2]+p*A[1..nrows(A),c1];
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   print(A);
   print(addcol(A,1,xy,2));
}
///////////////////////////////////////////////////////////////////////////////

proc addrow (matrix A, int r1, poly p, int r2)
"USAGE:   addcol(A,r1,p,r2);  A matrix, p poly, r1, r2 positive integers
RETURN:  matrix,  A being modified by adding p times row r1 to row r2
EXAMPLE: example addrow; shows an example
"
{
   A[r2,1..ncols(A)]=A[r2,1..ncols(A)]+p*A[r1,1..ncols(A)];
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   print(A);
   print(addrow(A,1,xy,3));
}
///////////////////////////////////////////////////////////////////////////////

proc multcol (matrix A, int c, poly p)
"USAGE:   addcol(A,c,p);  A matrix, p poly, c positive integer
RETURN:  matrix,  A being modified by multiplying column c with p
EXAMPLE: example multcol; shows an example
"
{
   A[1..nrows(A),c]=p*A[1..nrows(A),c];
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   print(A);
   print(multcol(A,2,xy));
}
///////////////////////////////////////////////////////////////////////////////

proc multrow (matrix A, int r, poly p)
"USAGE:   addcol(A,r,p);  A matrix, p poly, r positive integer
RETURN:  matrix,  A being modified by multiplying row r with p
EXAMPLE: example multrow; shows an example
"
{
   A[r,1..ncols(A)]=p*A[r,1..ncols(A)];
   return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,2,3,4,5,6,7,8,9;
   print(A);
   print(multrow(A,2,xy));
}
///////////////////////////////////////////////////////////////////////////////

proc permcol (matrix A, int c1, int c2)
"USAGE:   permcol(A,c1,c2);  A matrix, c1,c2 positive integers
RETURN:  matrix,  A being modified by permuting column c1 and c2
EXAMPLE: example permcol; shows an example
"
{
   matrix B=A;
   B[1..nrows(B),c1]=A[1..nrows(A),c2];
   B[1..nrows(B),c2]=A[1..nrows(A),c1];
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,x,3,4,y,6,7,z,9;
   print(A);
   print(permcol(A,2,3));
}
///////////////////////////////////////////////////////////////////////////////

proc permrow (matrix A, int r1, int r2)
"USAGE:   permrow(A,r1,r2);  A matrix, r1,r2 positive integers
RETURN:  matrix,  A being modified by permuting row r1 and r2
EXAMPLE: example permrow; shows an example
"
{
   matrix B=A;
   B[r1,1..ncols(B)]=A[r2,1..ncols(A)];
   B[r2,1..ncols(B)]=A[r1,1..ncols(A)];
   return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(x,y,z),lp;
   matrix A[3][3]=1,2,3,x,y,z,7,8,9;
   print(A);
   print(permrow(A,2,1));
}
///////////////////////////////////////////////////////////////////////////////

proc rowred (matrix A,list #)
"USAGE:   rowred(A[,any]);  A matrix
RETURN:  matrix B, being the row reduced form of A,
         if any defined: a list of two matrices B,C, such that B = C * A
EXAMPLE: example rowred; shows an example
"
{
   int m,n=nrows(A),ncols(A);
   int i,j,k,l,rk;
   poly p;
   matrix d[m][n];
   for (i=1;i<=m;i=i+1)
   {  for (j=1;j<=n;j=j+1)
      {  p = A[i,j];
         if (ord(p)==0)
         {  if (deg(p)==0) { d[i,j]=p; }
         }
      }
   }
   matrix b = A;
   if (size(#)) { b = concat(b,unitmat(m)); }
   for (l=1;l<=n;l=l+1)
   {
      k  = findfirst(ideal(d[l]),rk+1);
      if (k)
      {  rk = rk+1;
         b  = permrow(b,rk,k); d = permrow(d,rk,k);
         p  = b[rk,l];         p = 1/p;
         b  = multrow(b,rk,p); d = multrow(d,rk,p);
         for (i=1;i<=m;i=i+1)
         {
            if (rk-i) { b = addrow(b,rk,-b[i,l],i);
                        d = addrow(d,rk,-d[i,l],i);
                      }
         }
      }
   }
   d = submat(b,1..m,1..n);
   if (size(#)) {return(d,submat(b,1..m,n+1..n+m));}
   return(d);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(A,B,C),dp;
   matrix m[12][14]=
   AC,BC,-3BC,0, B2-A2,-3AC,B2,B2,0, 0, -C2,0, 0, 0,
   2C,0, 0,   B, -A,   -4C, 2A,0, 0, 0, 0,  0, 0, 0,
   0, 2C,-4C, -A,B,    0,   B, 3B,AB,B2,0,  0, 0, 0,
   0, 0, A,   0, 0,    B,   0, 0, 0, 0, 0,  0, 0, -C2,
   0, 0, 0,   0, 0,    0,   2, 0, 2A,0, 0,  0, 0, 0,
   0, 0, 0,   0, 0,    0,   0, 2, A, 3B,0,  B2,0, 0,
   0, 0, 0,   0, 0,    0,   0, 0, 2, 0, 0,  0, 0, 0,
   0, 0, 0,   0, 0,    0,   0, 0, 0, 2, 0,  3B,B2,0,
   0, 0, 0,   0, 0,    0,   0, 0, 0, 0, 1,  0, 0, 0,
   0, 0, 0,   0, 0,    0,   0, 0, 0, 0, 0,  2, 3B,0,
   0, 0, 0,   0, 0,    0,   0, 0, 0, 0, 0,  0, 2, 0,
   0, 0, 0,   0, 0,    0,   0, 0, 0, 0, 0,  0, 0, 1;
   print(rowred(m));
   list L=rowred(m,1);
   print(L[1]);"---------";print(L[2]);
}
///////////////////////////////////////////////////////////////////////////////

proc colred (matrix A,list #)
"USAGE:   colred(A[,any]);  A matrix
RETURN:  matrix B, being the column reduced form of A,
         if any defined: a list of two matrices B,C, such that B =  A * C
EXAMPLE: example colred; shows an example
"
{
   A = transpose(A);
   if (size(#))
   { list L = rowred(A,1); return(transpose(L[1]),transpose(L[2]));}
   else
   { return(transpose(rowred(A)));}
}
example
{ "EXAMPLE:"; echo = 2;
 ring r=0,(A,B,C),dp;
   matrix m[14][12]=
   AC,    2C, 0,  0,  0, 0, 0,0, 0,0, 0,0,
   BC,    0,  2C, 0,  0, 0, 0,0, 0,0, 0,0,
   -3BC,  0,  -4C,A,  0, 0, 0,0, 0,0, 0,0,
   0,     B,  -A, 0,  0, 0, 0,0, 0,0, 0,0,
   -A2+B2,-A, B,  0,  0, 0, 0,0, 0,0, 0,0,
   -3AC,  -4C,0,  B,  0, 0, 0,0, 0,0, 0,0,
   B2,    2A, B,  0,  2, 0, 0,0, 0,0, 0,0,
   B2,    0,  3B, 0,  0, 2, 0,0, 0,0, 0,0,
   0,     0,  AB, 0,  2A,A, 2,0, 0,0, 0,0,
   0,     0,  B2, 0,  0, 3B,0,2, 0,0, 0,0,
   -C2,   0,  0,  0,  0, 0, 0,0, 1,0, 0,0,
   0,     0,  0,  0,  0, B2,0,3B,0,2, 0,0,
   0,     0,  0,  0,  0, 0, 0,B2,0,3B,2,0,
   0,     0,  0,  -C2,0, 0, 0,0, 0,0, 0,1;
   print(colred(m));
   list L=colred(m,1);
   print(L[1]);"---------";print(L[2]);
}
//////////////////////////////////////////////////////////////////////////////

static proc findfirst (ideal i,int t)
{
   int n,k;
   for (n=t;n<=ncols(i);n=n+1)
   {
      if (i[n]!=0) { k=n;break;}
   }
   return(k);
}
//////////////////////////////////////////////////////////////////////////////

proc rm_unitcol(matrix A)
"USAGE:   rm_unitcol(A); A matrix (being row reduced)
RETURN:  matrix, obtained from A by deleting unit columns (having just one 1
         and else 0 as entries) and associated rows
EXAMPLE: example rm_unitcol; shows an example
"
{
 int l,j;
 intvec v;
 for (j=1;j<=ncols(A);j=j+1)
 {
    if (gen(l+1)==module(A)[j]) {l=l+1;}
    else { v=v,j;}
 }
 if (size(v)>1)
    {  v = v[2..size(v)];
       return(submat(A,l+1..nrows(A),v));
    }
 else
    { return(0);}
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(A,B,C),dp;
   matrix m[14][12]=
   1/2B2,   A,  1/2B,    0,  1,0,0,0,0,0,0,0,
   1/2B2,   0,  3/2B,    0,  0,1,0,0,0,0,0,0,
   -3/4AB2, -A2,-3/4AB,  0,  0,0,1,0,0,0,0,0,
   -3/4B3,  0,  -7/4B2,  0,  0,0,0,1,0,0,0,0,
   -C2,     0,  0,       0,  0,0,0,0,1,0,0,0,
   7/8B4,   0,  15/8B3,  0,  0,0,0,0,0,1,0,0,
   -15/16B5,0,  -31/16B4,0,  0,0,0,0,0,0,1,0,
   0,       0,  0,       -C2,0,0,0,0,0,0,0,1,
   -3BC,    0,  -4C,     A,  0,0,0,0,0,0,0,0,
   0,       B,  -A,      0,  0,0,0,0,0,0,0,0,
   -A2+B2,  -A, B,       0,  0,0,0,0,0,0,0,0,
   -3AC,    -4C,0,       B,  0,0,0,0,0,0,0,0,
   AC,      2C, 0,       0,  0,0,0,0,0,0,0,0,
   BC,      0,  2C,      0,  0,0,0,0,0,0,0,0;
   print(rm_unitcol(m));
}
//////////////////////////////////////////////////////////////////////////////

proc rm_unitrow (matrix A)
"USAGE:   rm_unitrow(A); A matrix (being colreduced)
RETURN:  matrix, obtained from A by deleting unit rows (having just one 1
         and else 0 as entries) and associated columns
EXAMPLE: example rm_unitrow; shows an example
"
{
 int l,j;
 intvec v;
 module M = transpose(A);
 for (j=1; j <= nrows(A); j=j+1)
 {
    if (gen(l+1) == M[j]) { l=l+1; }
    else { v=v,j; }
 }
 if (size(v) > 1)
    {  v = v[2..size(v)];
       return(submat(A,v,l+1..ncols(A)));
    }
 else
    { return(0);}
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(A,B,C),dp;
   matrix m[12][14]=
   1/2B2,1/2B2,-3/4AB2,-3/4B3,-C2,7/8B4, -15/16B5,0,  -3BC,0, B2-A2,-3AC,AC,BC,
   A,    0,    -A2,    0,     0,  0,     0,       0,  0,   B, -A,   -4C, 2C,0,
   1/2B, 3/2B, -3/4AB, -7/4B2,0,  15/8B3,-31/16B4,0,  -4C, -A,B,    0,   0, 2C,
   0,    0,    0,      0,     0,  0,     0,       -C2,A,   0, 0,    B,   0, 0,
   1,    0,    0,      0,     0,  0,     0,       0,  0,   0, 0,    0,   0, 0,
   0,    1,    0,      0,     0,  0,     0,       0,  0,   0, 0,    0,   0, 0,
   0,    0,    1,      0,     0,  0,     0,       0,  0,   0, 0,    0,   0, 0,
   0,    0,    0,      1,     0,  0,     0,       0,  0,   0, 0,    0,   0, 0,
   0,    0,    0,      0,     1,  0,     0,       0,  0,   0, 0,    0,   0, 0,
   0,    0,    0,      0,     0,  1,     0,       0,  0,   0, 0,    0,   0, 0,
   0,    0,    0,      0,     0,  0,     1,       0,  0,   0, 0,    0,   0, 0,
   0,    0,    0,      0,     0,  0,     0,       1,  0,   0, 0,    0,   0, 0;
   print(rm_unitrow(m));
}
//////////////////////////////////////////////////////////////////////////////