///////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Linear Algebra";
info="
LIBRARY:  matrix.lib    Elementary Matrix Operations

PROCEDURES:
 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 a matrix into col-reduced Gauss normal form
 gauss_row(A);          transform a matrix into row-reduced Gauss normal form
 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
 linear_relations(E);   find linear relations between homogeneous vectors
 rm_unitrow(A);         remove unit rows and associated columns of A
 rm_unitcol(A);         remove unit columns and associated rows of A
 headStand(A);          A[n-i+1,m-j+1]:=A[i,j]
 symmetricBasis(n,k[,s]); basis of k-th symmetric power of n-dim v.space
 exteriorBasis(n,k[,s]); basis of k-th exterior power of n-dim v.space
 symmetricPower(A,k);   k-th symmetric power of a module/matrix A
 exteriorPower(A,k);    k-th exterior power of a module/matrix A
          (parameters in square brackets [] are optional)
";

LIB "inout.lib";
LIB "ring.lib";
LIB "random.lib";
LIB "general.lib"; // for sort
LIB "nctools.lib"; // for superCommutative

///////////////////////////////////////////////////////////////////////////////

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);
      if ( size(m) == 0)
      { intmat B; }
      else
      { intmat B[nrows(A)][size(m)]; }
      int i,j;
      for( i=1; i<=ncols(A); i++ )
      {
         if( m[i]!=[0] )
         {
            j++;
            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;
   for( i=size(#);i>0; i-- ) { #[i]=module(#[i]); }
   module B=#[1..size(#)];
   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 unit matrix 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++ ) { 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][2] = 1,2,3,4,5,6;
   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++ ) { N=N+nrows(#[i]); }
   for( i=1; i<=size(#); i++ )
   {
      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;
   print(A);
   print(B);
   print(dsum(A,B));
}
///////////////////////////////////////////////////////////////////////////////

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[2][3] = 1,2,x,y,z,7;
   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.
NOTE:    if id has less than nxm elements, the matrix is filled with 0's,
         (default: id=maxideal(1)).
         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(3,3));      // the generic 3x3 matrix
   ring R1 = 0,(a,b,c,d),dp;
   matrix A = genericmat(3,4,maxideal(1)^3);
   print(A);
   int n,m = 3,2;
   ideal i = ideal(randommat(1,n*m,maxideal(1),9));
   print(genericmat(n,m,i));    // matrix of generic linear forms
}
///////////////////////////////////////////////////////////////////////////////

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.
         If printlevel > 0, the position where c is not a complex is shown.
EXAMPLE: example is_complex; shows an example
"
{
   int i;
   module @test;
   for( i=1; i<=size(c)-1; i++ )
   {
      c[i]=matrix(c[i]); c[i+1]=matrix(c[i+1]);
      @test=c[i]*c[i+1];
      if (size(@test)!=0)
      {
        dbprint(printlevel-voice+2,"// not a complex at position " +string(i));
         return(0);
      }
   }
   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 (tensor) 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=ncols(A);i>0; i-- )
     {
       for( j=1; j<=nrows(A); j++ )
       {
          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>=0
RETURN:  intmat resp. matrix, the n-th power of A
NOTE:    for A=intmat 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" )
   {
      ERROR("no matrix or intmat!");
   }
   if( ncols(A) != nrows(A) )
   {
      ERROR("not a square matrix!");
   }
//---------------------------- 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[3][3]=0,x,y,z,0,0,y,z,0;
   print(power(B,3));"";
}
///////////////////////////////////////////////////////////////////////////////

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))
         skewmat(n); creates the generic skew-symmetric matrix
NOTE:    if id has less than n*(n-1)/2 elements, the matrix is
         filled with 0's,
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++ )
   {
      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
   ring R1 = 0,(a,b,c),dp;
   matrix A=skewmat(4,maxideal(1)^2);
   print(A);
   int n=3;
   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++ )
   {
      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++ ) {  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
   ring R1 = 0,(a,b,c),dp;
   matrix A=symmat(4,maxideal(1)^3);
   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
"
{
   if (ncols(A)==0)
   {
     int q=nrows(A)*nrows(B);
     matrix D[q][0];
     return(D);
   }

   int i,j;
   matrix C,D;
   for( i=1; i<=nrows(A); i++ )
   {
     C = A[i,1]*B;
     for( j=2; j<=ncols(A); j++ )
     {
       C = concat(C,A[i,j]*B);
     }
     D = concat(D,transpose(C));
   }
   D = transpose(D);
   return(submat(D,2..nrows(D),1..ncols(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, list #)
"USAGE:   gauss_col(A[,e]); A a matrix, e any type
RETURN:  - a matrix B, if called with one argument; B is the complete column-
           reduced upper-triangular normal form of A if A is constant,
           (resp. as far as this is possible if A is a polynomial matrix;
           no division by polynomials).
@*       - a list L of two matrices, if called with two arguments;
           L satisfies L[1] = A * L[2] with L[1] the column-reduced form of A
           and L[2] the transformation matrix.
NOTE:    * The procedure just applies interred to A with ordering (C,dp).
           The transformation matrix is obtained by applying 'lift'.
           This should be faster than the procedure colred.
@*       * It should only be used with exact coefficient field (there is no
           pivoting and rounding error treatment).
@*       * Parameters are allowed. Hence, if the entries of A are parameters,
           B is the column-reduced form of A over the rational function field.
SEE ALSO:  colred
EXAMPLE: example gauss_col; shows an example
"
{
   def R=basering; int u;
   string mp = string(minpoly);
   int n = nrows(A);
   int m = ncols(A);
   module M = A;
   intvec v = option(get);
//------------------------ change ordering if necessary ----------------------
   if( ordstr(R) != ("C,dp("+string(nvars(R))+")") )
   {
     def @R=changeord("C,dp",R); setring @R; u=1;
     execute("minpoly="+mp+";");
     matrix A = imap(R,A);
     module M = A;
   }
//------------------------------ start computation ---------------------------
   option(redSB);
   M = simplify(interred(M),1);
   if(size(#) != 0)
   {
      module N = lift(A,M);
   }
//--------------- reset ring and options and return --------------------------
   if ( u==1 )
   {
      setring R;
      M=imap(@R,M);
      if (size(#) != 0)
      {
         module N = imap(@R,N);
      }
      kill @R;
   }
   option(set,v);
   // M = sort(M,size(M)..1)[1];
   A = matrix(M,n,m);
   if (size(#) != 0)
   {
     list L= A,matrix(N,m,m);
     return(L);
   }
   return(matrix(M,n,m));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=(0,a,b),(A,B,C),dp;
   matrix m[8][6]=
   0,    2*C, 0,    0,  0,   0,
   0,    -4*C,a*A,  0,  0,   0,
   b*B,  -A,  0,    0,  0,   0,
   -A,   B,   0,    0,  0,   0,
   -4*C, 0,   B,    2,  0,   0,
   2*A,  B,   0,    0,  0,   0,
   0,    3*B, 0,    0,  2b,  0,
   0,    AB,  0,    2*A,A,   2a;"";
   list L=gauss_col(m,1);
   print(L[1]);
   print(L[2]);

   ring S=0,x,(c,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, list #)
"USAGE:  gauss_row(A [,e]); A matrix, e any type
RETURN: - a matrix B, if called with one argument; B is the complete row-
          reduced lower-triangular normal form of A if A is constant,
          (resp. as far as this is possible if A is a polynomial matrix;
          no division by polynomials).
@*      - a list L of two matrices, if called with two arguments;
          L satisfies transpose(L[2])*A=transpose(L[1])
          with L[1] the row-reduced form of A
          and L[2] the transformation matrix.
NOTE:   * This procedure just applies gauss_col to the transposed matrix.
          The transformation matrix is obtained by applying lift.
          This should be faster than the procedure rowred.
@*      * It should only be used with exact coefficient field (there is no
          pivoting and rounding error treatment).
@*      * Parameters are allowed. Hence, if the entries of A are parameters,
          B is the row-reduced form of A over the rational function field.
SEE ALSO: rowred
EXAMPLE: example gauss_row; shows an example
"
{
   if(size(#) > 0)
   {
     list L = gauss_col(transpose(A),1);
     return(L);
   }
   A = gauss_col(transpose(A));
   return(transpose(A));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=(0,a,b),(A,B,C),dp;
   matrix m[6][8]=
   0, 0,  b*B, -A,-4C,2A,0, 0,
   2C,-4C,-A,B, 0,  B, 3B,AB,
   0,a*A,  0, 0, B,  0, 0, 0,
   0, 0,  0, 0, 2,  0, 0, 2A,
   0, 0,  0, 0, 0,  0, 2b, A,
   0, 0,  0, 0, 0,  0, 0, 2a;"";
   print(gauss_row(m));"";
   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;
   list L=gauss_row(A,1);
   print(L[1]);
   print(L[2]);
}
///////////////////////////////////////////////////////////////////////////////

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
"
{
   int k=nrows(A);
   A[1..k,c2]=A[1..k,c2]+p*A[1..k,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:   addrow(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
"
{
   int k=ncols(A);
   A[r2,1..k]=A[r2,1..k]+p*A[r1,1..k];
   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:   multcol(A,c,p);  A matrix, p poly, c positive integer
RETURN:  matrix,  A being modified by multiplying column c by p
EXAMPLE: example multcol; shows an example
"
{
   int k=nrows(A);
   A[1..k,c]=p*A[1..k,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:   multrow(A,r,p);  A matrix, p poly, r positive integer
RETURN:  matrix,  A being modified by multiplying row r by p
EXAMPLE: example multrow; shows an example
"
{
   int k=ncols(A);
   A[r,1..k]=p*A[r,1..k];
   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 columns c1 and c2
EXAMPLE: example permcol; shows an example
"
{
   matrix B=A;
   int k=nrows(B);
   B[1..k,c1]=A[1..k,c2];
   B[1..k,c2]=A[1..k,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 rows r1 and r2
EXAMPLE: example permrow; shows an example
"
{
   matrix B=A;
   int k=ncols(B);
   B[r1,1..k]=A[r2,1..k];
   B[r2,1..k]=A[r1,1..k];
   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[,e]);  A matrix, e any type
RETURN:  - a matrix B, being the row reduced form of A, if rowred is called
           with one argument.
           (as far as this is possible over the polynomial ring; no division
           by polynomials)
@*       - a list L of two matrices, such that L[1] = L[2] * A with L[1]
           the row-reduced form of A and L[2] the transformation matrix
           (if rowred is called with two arguments).
ASSUME:  The entries of A are in the base field. It is not verified whether
         this assumption holds.
NOTE:    * This procedure is designed for teaching purposes mainly.
@*       * The straight forward Gaussian algorithm is implemented in the
           library (no standard basis computation).
           The transformation matrix is obtained by concatenating a unit
           matrix to A. proc gauss_row should be faster.
@*       * It should only be used with exact coefficient field (there is no
           pivoting) over the polynomial ring (ordering lp or dp).
@*       * Parameters are allowed. Hence, if the entries of A are parameters
           the computation takes place over the field of rational functions.
SEE ALSO:  gauss_row
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++)
   {  for (j=1;j<=n;j++)
      {  p = A[i,j];
         if (ord(p)==0)
         {  if (deg(p)==0) { d[i,j]=p; }
         }
      }
   }
   matrix b = A;
   if (size(#) != 0) { b = concat(b,unitmat(m)); }
      for (l=1;l<=n;l=l+1)
   {  pmat(d);
      k  = findfirst(ideal(d[l]),rk+1);
      if (k)
      {  rk = rk+1;
         b  = permrow(b,rk,k);
         p  = b[rk,l];         p = 1/p;
         b  = multrow(b,rk,p);
         for (i=1;i<=m;i++)
         {
            if (rk-i) { b = addrow(b,rk,-b[i,l],i);}
         }
         d = 0;
         for (i=rk+1;i<=m;i++)
         {  for (j=l+1;j<=n;j++)
            {  p = b[i,j];
               if (ord(p)==0)
               {  if (deg(p)==0) { d[i,j]=p; }
               }
            }
         }

      }
   }
   d = submat(b,1..m,1..n);
   if (size(#))
   {
      list L=d,submat(b,1..m,n+1..n+m);
      return(L);
   }
   return(d);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=(0,a,b),(A,B,C),dp;
   matrix m[6][8]=
   0, 0,  b*B, -A,-4C,2A,0, 0,
   2C,-4C,-A,B, 0,  B, 3B,AB,
   0,a*A,  0, 0, B,  0, 0, 0,
   0, 0,  0, 0, 2,  0, 0, 2A,
   0, 0,  0, 0, 0,  0, 2b, A,
   0, 0,  0, 0, 0,  0, 0, 2a;"";
   print(rowred(m));"";
   list L=rowred(m,1);
   print(L[1]);
   print(L[2]);
}
///////////////////////////////////////////////////////////////////////////////

proc colred (matrix A,list #)
"USAGE:   colred(A[,e]);  A matrix, e any type
RETURN:  - a matrix B, being the column reduced form of A, if colred is
           called with one argument.
           (as far as this is possible over the polynomial ring;
           no division by polynomials)
@*       - a list L of two matrices, such that L[1] = A * L[2] with L[1]
           the column-reduced form of A and L[2] the transformation matrix
           (if colred is called with two arguments).
ASSUME:  The entries of A are in the base field. It is not verified whether
         this assumption holds.
NOTE:    * This procedure is designed for teaching purposes mainly.
@*       * It applies rowred to the transposed matrix.
           proc gauss_col should be faster.
@*       * It should only be used with exact coefficient field (there is no
           pivoting) over the polynomial ring (ordering lp or dp).
@*       * Parameters are allowed. Hence, if the entries of A are parameters
           the computation takes place over the field of rational functions.
SEE ALSO:  gauss_col
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),(A,B,C),dp;
   matrix m[8][6]=
   0,    2*C, 0,    0,  0,   0,
   0,    -4*C,a*A,  0,  0,   0,
   b*B,  -A,  0,    0,  0,   0,
   -A,   B,   0,    0,  0,   0,
   -4*C, 0,   B,    2,  0,   0,
   2*A,  B,   0,    0,  0,   0,
   0,    3*B, 0,    0,  2b,  0,
   0,    AB,  0,    2*A,A,   2a;"";
   print(colred(m));"";
   list L=colred(m,1);
   print(L[1]);
   print(L[2]);
}
//////////////////////////////////////////////////////////////////////////////

proc linear_relations(module M)
"USAGE:   linear_relations(M);
         M: a module
ASSUME:  All non-zero entries of M are homogeneous polynomials of the same
         positive degree. The base field must be an exact field (not real
         or complex).
         It is not checked whether these assumptions hold.
RETURN:  a maximal module R such that M*R is formed by zero vectors.
EXAMPLE: example linear_relations; shows an example.
"
{ int n = ncols(M);
  def BaseR = basering;
  def br = changeord("dp",basering);
  setring br;
  module M = imap(BaseR,M);
  ideal vars = maxideal(1);
  ring tmpR = 0, ('y(1..n)), dp;
  def newR = br + tmpR;
  setring newR;
  module M = imap(br,M);
  ideal vars = imap(br,vars);
  attrib(vars,"isSB",1);
  for (int i = 1; i<=n; i++) {
    M[i] = M[i] + 'y(i)*gen(1);
  }
  M = interred(M);
  module redM = NF(M,vars);
  module REL;
  int sizeREL;
  int j;
  for (i=1; i<=n; i++) {
    if (M[i][1]==redM[i][1]) { //-- relation found!
      sizeREL++;
      REL[sizeREL]=0;
      for (j=1; j<=n; j++) {
        REL[sizeREL] = REL[sizeREL] + (M[i][1]/'y(j))*gen(j);
      }
    }
  }
  setring BaseR;
  return(minbase(imap(newR,REL)));
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (3,w), (a,b,c,d),dp;
  minpoly = w2-w-1;
  module M = [a2,b2],[wab,w2c2+2b2],[(w-2)*a2+wab,wb2+w2c2];
  module REL = linear_relations(M);
  pmat(REL);
  pmat(M*REL);
}

//////////////////////////////////////////////////////////////////////////////

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++)
 {
    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[6][8]=
   0,  0,    A,   0, 1,0,  0,0,
   0,  0,  -C2,   0, 0,0,  1,0,
   0,  0,    0,1/2B, 0,0,  0,1,
   0,  0,    B,  -A, 0,2A, 0,0,
   2C,-4C,  -A,   B, 0,B,  0,0,
   0,  A,    0,   0, 0,0,  0,0;
   print(rm_unitcol(m));
}
//////////////////////////////////////////////////////////////////////////////

proc rm_unitrow (matrix A)
"USAGE:   rm_unitrow(A); A matrix (being col-reduced)
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++)
 {
    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[8][6]=
   0,0,  0,   0, 2C, 0,
   0,0,  0,   0, -4C,A,
   A,-C2,0,   B, -A, 0,
   0,0,  1/2B,-A,B,  0,
   1,0,  0,   0, 0,  0,
   0,0,  0,   2A,B,  0,
   0,1,  0,   0, 0,  0,
   0,0,  1,   0, 0,  0;
   print(rm_unitrow(m));
}
//////////////////////////////////////////////////////////////////////////////
proc headStand(matrix M)
"USAGE:   headStand(M);  M matrix
RETURN:  matrix B such that B[i][j]=M[n-i+1,m-j+1], n=nrows(M), m=ncols(M)
EXAMPLE: example headStand; shows an example
"
{
  int i,j;
  int n=nrows(M);
  int m=ncols(M);
  matrix B[n][m];
  for(i=1;i<=n;i++)
  {
     for(j=1;j<=m;j++)
     {
        B[n-i+1,m-j+1]=M[i,j];
     }
  }
  return(B);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(A,B,C),dp;
   matrix M[2][3]=
   0,A,  B,
   A2, B2, C;
   print(M);
   print(headStand(M));
}
//////////////////////////////////////////////////////////////////////////////

// Symmetric/Exterior powers thanks to Oleksandr Iena for his persistence ;-)

proc symmetricBasis(int n, int k, list #)
"USAGE:    symmetricBasis(n, k[,s]); n int, k int, s string
RETURN:   ring, poynomial ring containing the ideal \"symBasis\",
          being a basis of the k-th symmetric power of an n-dim vector space.
NOTE:     The output polynomial ring has characteristics 0 and n variables
          named \"S(i)\", where the base variable name S is either given by the
          optional string argument(which must not contain brackets) or equal to
          "e" by default.
SEE ALSO: exteriorBasis
KEYWORDS: symmetric basis
EXAMPLE:  example symmetricBasis; shows an example"
{
//------------------------ handle optional base variable name---------------
  string S = "e";
  if( size(#) > 0 )
  {
    if( typeof(#[1]) != "string" )
    {
      ERROR("Wrong optional argument: must be a string");
    }
    S = #[1];
    if( (find(S, "(") + find(S, ")")) > 0 )
    {
      ERROR("Wrong optional argument: must be a string without brackets");
    }
  }

//------------------------- create ring container for symmetric power basis-
  execute("ring @@@SYM_POWER_RING_NAME=(0),("+S+"(1.."+string(n)+")),dp;");

//------------------------- choose symmetric basis -------------------------
  ideal symBasis = maxideal(k);

//------------------------- export and return      -------------------------
  export symBasis;
  return(basering);
}
example
{ "EXAMPLE:"; echo = 2;

// basis of the 3-rd symmetricPower of a 4-dim vector space:
def R = symmetricBasis(4, 3, "@e"); setring R;
R;  // container ring:
symBasis; // symmetric basis:
}

//////////////////////////////////////////////////////////////////////////////

proc exteriorBasis(int n, int k, list #)
"USAGE:    exteriorBasis(n, k[,s]); n int, k int, s string
RETURN:   qring, an exterior algebra containing the ideal \"extBasis\",
          being a basis of the k-th exterior power of an n-dim vector space.
NOTE:     The output polynomial ring has characteristics 0 and n variables
          named \"S(i)\", where the base variable name S is either given by the
          optional string argument(which must not contain brackets) or equal to
          "e" by default.
SEE ALSO: symmetricBasis
KEYWORDS: exterior basis
EXAMPLE:  example exteriorBasis; shows an example"
{
//------------------------ handle optional base variable name---------------
  string S = "e";
  if( size(#) > 0 )
  {
    if( typeof(#[1]) != "string" )
    {
      ERROR("Wrong optional argument: must be a string");
    }
    S = #[1];
    if( (find(S, "(") + find(S, ")")) > 0 )
    {
      ERROR("Wrong optional argument: must be a string without brackets");
    }
  }

//------------------------- create ring container for symmetric power basis-
  execute("ring @@@EXT_POWER_RING_NAME=(0),("+S+"(1.."+string(n)+")),dp;");

//------------------------- choose exterior basis -------------------------
  def T = superCommutative(); setring T;
  ideal extBasis = simplify( NF(maxideal(k), std(0)), 1 + 2 + 8 );

//------------------------- export and return      -------------------------
  export extBasis;
  return(basering);
}
example
{ "EXAMPLE:"; echo = 2;
// basis of the 3-rd symmetricPower of a 4-dim vector space:
def r = exteriorBasis(4, 3, "@e"); setring r;
r; // container ring:
extBasis; // exterior basis:
}

//////////////////////////////////////////////////////////////////////////////

static proc chooseSafeVarName(string prefix, string suffix)
"USAGE: give appropreate prefix for variable names
RETURN: safe variable name (repeated prefix + suffix)
"
{
  string V = varstr(basering);
  string S = suffix;
  while( find(V, S) > 0 )
  {
    S = prefix + S;
  }
  return(S);
}

//////////////////////////////////////////////////////////////////////////////

static proc mapPower(int p, module A, int k, def Tn, def Tm)
"USAGE: by both symmetric- and exterior-Power"
NOTE: everything over the basering!
      module A (matrix of the map), int k (power)
      rings Tn is source- and Tm is image-ring with bases
          resp. Ink and Imk.
      M = max dim of Image, N - dim. of source
SEE ALSO: symmetricPower, exteriorPower"
{
  def save = basering;

  int n = nvars(save);
  int M = nrows(A);
  int N = ncols(A);

  int i, j;

//------------------------- compute matrix of single images ------------------
  def Rm = save + Tm;  setring Rm;
  dbprint(p-2, "Temporary Working Ring", Rm);

  module A = imap(save, A);

  ideal B; poly t;

  for( i = N; i > 0; i-- )
  {
    t = 0;
    for( j = M; j > 0; j-- )
    {
      t = t + A[i][j] * var(n + j);
    }

    B[i] = t;
  }

  dbprint(p-1, "Matrix of single images", B);

//------------------------- compute image ---------------------
  // apply S^k(A): Tn -> Rm  to Source basis vectors Ink:
  map TMap = Tn, B;

  ideal C = NF(TMap(Ink), std(0));
  dbprint(p-1, "Image Matrix: ", C);


//------------------------- write it in Image basis ---------------------
  ideal Imk = imap(Tm, Imk);

  module D; poly lm; vector tt;

  for( i = ncols(C); i > 0; i-- )
  {
    t = C[i];
    tt = 0;

    while( t != 0 )
    {
      lm = leadmonom(t);
      //    lm;
      for( j = ncols(Imk); j > 0; j-- )
      {
        if( lm / Imk[j] != 0 )
        {
          tt = tt + (lead(t) / Imk[j]) * gen(j);
          break;
        }
      }
      t = t - lead(t);
    }

    D[i] = tt;
  }

//------------------------- map it back and return  ---------------------
  setring save;
  return( imap(Rm, D) );
}


//////////////////////////////////////////////////////////////////////////////

proc symmetricPower(module A, int k)
"USAGE:    symmetricPower(A, k); A module, k int
RETURN:   module: the k-th symmetric power of A
NOTE:     the chosen bases and most of intermediate data will be shown if
          printlevel is big enough
SEE ALSO: exteriorPower
KEYWORDS: symmetric power
EXAMPLE:  example symmetricPower; shows an example"
{
  int p = printlevel - voice + 2;

  def save = basering;

  int M = nrows(A);
  int N = ncols(A);

  string S = chooseSafeVarName("@", "@_e");

//------------------------- choose source basis -------------------------
  def Tn = symmetricBasis(N, k, S); setring Tn;
  ideal Ink = symBasis;
  export Ink;
  dbprint(p-3, "Temporary Source Ring", basering);
  dbprint(p, "S^k(Source Basis)", Ink);

//------------------------- choose image basis -------------------------
  def Tm = symmetricBasis(M, k, S); setring Tm;
  ideal Imk = symBasis;
  export Imk;
  dbprint(p-3, "Temporary Image Ring", basering);
  dbprint(p, "S^k(Image Basis)", Imk);

//------------------------- compute and return S^k(A) in chosen bases --
  setring save;

  return(mapPower(p, A, k, Tn, Tm));
}
example
{ "EXAMPLE:"; echo = 2;

ring r = (0),(a, b, c, d), dp; r;
module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2); print(B);

// symmetric power over a commutative K-algebra:
print(symmetricPower(B, 2));
print(symmetricPower(B, 3));

// symmetric power over an exterior algebra:
def g = superCommutative(); setring g; g;

module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2); print(B);

print(symmetricPower(B, 2)); // much smaller!
print(symmetricPower(B, 3)); // zero! (over an exterior algebra!)

}

//////////////////////////////////////////////////////////////////////////////

proc exteriorPower(module A, int k)
"USAGE:    exteriorPower(A, k); A module, k int
RETURN:   module: the k-th exterior power of A
NOTE:     the chosen bases and most of intermediate data will be shown if
          printlevel is big enough. Last rows will be invisible if zero.
SEE ALSO: symmetricPower
KEYWORDS: exterior power
EXAMPLE:  example exteriorPower; shows an example"
{
  int p = printlevel - voice + 2;
  def save = basering;

  int M = nrows(A);
  int N = ncols(A);

  string S = chooseSafeVarName("@", "@_e");

//------------------------- choose source basis -------------------------
  def Tn = exteriorBasis(N, k, S); setring Tn;
  ideal Ink = extBasis;
  export Ink;
  dbprint(p-3, "Temporary Source Ring", basering);
  dbprint(p, "E^k(Source Basis)", Ink);

//------------------------- choose image basis -------------------------
  def Tm = exteriorBasis(M, k, S); setring Tm;
  ideal Imk = extBasis;
  export Imk;
  dbprint(p-3, "Temporary Image Ring", basering);
  dbprint(p, "E^k(Image Basis)", Imk);

//------------------------- compute and return E^k(A) in chosen bases --
  setring save;
  return(mapPower(p, A, k, Tn, Tm));
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (0),(a, b, c, d, e, f), dp;
  r; "base ring:";

  module B = a*gen(1) + c*gen(2) + e*gen(3),
             b*gen(1) + d*gen(2) + f*gen(3),
                        e*gen(1) + f*gen(3);

  print(B);
  print(exteriorPower(B, 2));
  print(exteriorPower(B, 3));

  def g = superCommutative(); setring g; g;

  module A = a*gen(1), b * gen(1), c*gen(2), d * gen(2);
  print(A);

  print(exteriorPower(A, 2));

  module B = a*gen(1) + c*gen(2) + e*gen(3),
             b*gen(1) + d*gen(2) + f*gen(3),
                        e*gen(1) + f*gen(3);
  print(B);

  print(exteriorPower(B, 2));
  print(exteriorPower(B, 3));

}

//////////////////////////////////////////////////////////////////////////////