/////////////////////////////////////////////////////////////////////////////// version="$Id: matrix.lib,v 1.26.2.2 2002/10/07 14:05:02 lossen Exp $"; 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 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); if ( size(m) == 0) { intmat B; } else { 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 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=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][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=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; 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=i+1 ) { 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=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: 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" ) { "// 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[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=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 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=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 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))+")") ) { changeord("@R","C,dp",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 L[1] = L[2] * A 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 " { 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: multrow(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[,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). 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=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(#) != 0) { 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); p = b[rk,l]; p = 1/p; b = multrow(b,rk,p); for (i=1;i<=m;i=i+1) { if (rk-i) { b = addrow(b,rk,-b[i,l],i);} } d = 0; for (i=rk+1;i<=m;i=i+1) { for (j=l+1;j<=n;j=j+1) { 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). 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]); } ////////////////////////////////////////////////////////////////////////////// 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[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=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[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)); } //////////////////////////////////////////////////////////////////////////////