Initial population.

1 files changed, 904 insertions, 0 deletions

diff --git a/src/algebra/matcat.spad.pamphlet b/src/algebra/matcat.spad.pamphlet
new file mode 100644
index 00000000..f7bc1ef8
--- /dev/null
+++ b/src/algebra/matcat.spad.pamphlet
@@ -0,0 +1,904 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra matcat.spad}
+\author{Johannes Grabmeier, Oswald Gschnitzer, Clifton J. Williamson}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{category MATCAT MatrixCategory}
+<<category MATCAT MatrixCategory>>=
+)abbrev category MATCAT MatrixCategory
+++ Authors: Grabmeier, Gschnitzer, Williamson
+++ Date Created: 1987
+++ Date Last Updated: July 1990
+++ Basic Operations:
+++ Related Domains: Matrix(R)
+++ Also See:
+++ AMS Classifications:
+++ Keywords: matrix, linear algebra
+++ Examples:
+++ References:
+++ Description:
+++   \spadtype{MatrixCategory} is a general matrix category which allows
+++   different representations and indexing schemes.  Rows and
+++   columns may be extracted with rows returned as objects of
+++   type Row and colums returned as objects of type Col.
+++   A domain belonging to this category will be shallowly mutable.
+++   The index of the 'first' row may be obtained by calling the
+++   function \spadfun{minRowIndex}.  The index of the 'first' column may
+++   be obtained by calling the function \spadfun{minColIndex}.  The index of
+++   the first element of a Row is the same as the index of the
+++   first column in a matrix and vice versa.
+MatrixCategory(R,Row,Col): Category == Definition where
+  R   : Ring
+  Row : FiniteLinearAggregate R
+  Col : FiniteLinearAggregate R
+
+  Definition ==> TwoDimensionalArrayCategory(R,Row,Col) with
+     shallowlyMutable
+       ++ One may destructively alter matrices
+
+     finiteAggregate
+       ++ matrices are finite
+
+--% Predicates
+
+     square?  : % -> Boolean
+       ++ \spad{square?(m)} returns true if m is a square matrix
+       ++ (i.e. if m has the same number of rows as columns) and false otherwise.
+     diagonal?: % -> Boolean
+       ++ \spad{diagonal?(m)} returns true if the matrix m is square and
+       ++ diagonal (i.e. all entries of m not on the diagonal are zero) and
+       ++ false otherwise.
+     symmetric?: % -> Boolean
+       ++ \spad{symmetric?(m)} returns true if the matrix m is square and
+       ++ symmetric (i.e. \spad{m[i,j] = m[j,i]} for all i and j) and false
+       ++ otherwise.
+     antisymmetric?: % -> Boolean
+       ++ \spad{antisymmetric?(m)} returns true if the matrix m is square and
+       ++ antisymmetric (i.e. \spad{m[i,j] = -m[j,i]} for all i and j) and false
+       ++ otherwise.
+
+--% Creation
+
+     zero: (NonNegativeInteger,NonNegativeInteger) -> %
+       ++ \spad{zero(m,n)} returns an m-by-n zero matrix.
+     matrix: List List R -> %
+       ++ \spad{matrix(l)} converts the list of lists l to a matrix, where the
+       ++ list of lists is viewed as a list of the rows of the matrix.
+     scalarMatrix: (NonNegativeInteger,R) -> %
+       ++ \spad{scalarMatrix(n,r)} returns an n-by-n matrix with r's on the
+       ++ diagonal and zeroes elsewhere.
+     diagonalMatrix: List R -> %
+       ++ \spad{diagonalMatrix(l)} returns a diagonal matrix with the elements
+       ++ of l on the diagonal.
+     diagonalMatrix: List % -> %
+       ++ \spad{diagonalMatrix([m1,...,mk])} creates a block diagonal matrix
+       ++ M with block matrices {\em m1},...,{\em mk} down the diagonal,
+       ++ with 0 block matrices elsewhere.
+       ++ More precisly: if \spad{ri := nrows mi}, \spad{ci := ncols mi},
+       ++ then m is an (r1+..+rk) by (c1+..+ck) - matrix  with entries
+       ++ \spad{m.i.j = ml.(i-r1-..-r(l-1)).(j-n1-..-n(l-1))}, if
+       ++ \spad{(r1+..+r(l-1)) < i <= r1+..+rl} and
+       ++ \spad{(c1+..+c(l-1)) < i <= c1+..+cl},
+       ++ \spad{m.i.j} = 0  otherwise.
+     coerce: Col -> %
+       ++ \spad{coerce(col)} converts the column col to a column matrix.
+     transpose: Row -> %
+       ++ \spad{transpose(r)} converts the row r to a row matrix.
+
+--% Creation of new matrices from old
+
+     transpose: % -> %
+       ++ \spad{transpose(m)} returns the transpose of the matrix m.
+     squareTop: % -> %
+       ++ \spad{squareTop(m)} returns an n-by-n matrix consisting of the first
+       ++ n rows of the m-by-n matrix m. Error: if
+       ++ \spad{m < n}.
+     horizConcat: (%,%) -> %
+       ++ \spad{horizConcat(x,y)} horizontally concatenates two matrices with
+       ++ an equal number of rows. The entries of y appear to the right
+       ++ of the entries of x.  Error: if the matrices
+       ++ do not have the same number of rows.
+     vertConcat: (%,%) -> %
+       ++ \spad{vertConcat(x,y)} vertically concatenates two matrices with an
+       ++ equal number of columns. The entries of y appear below
+       ++ of the entries of x.  Error: if the matrices
+       ++ do not have the same number of columns.
+
+--% Part extractions/assignments
+
+     listOfLists: % -> List List R
+       ++ \spad{listOfLists(m)} returns the rows of the matrix m as a list
+       ++ of lists.
+     elt: (%,List Integer,List Integer) -> %
+       ++ \spad{elt(x,rowList,colList)} returns an m-by-n matrix consisting
+       ++ of elements of x, where \spad{m = # rowList} and \spad{n = # colList}.
+       ++ If \spad{rowList = [i<1>,i<2>,...,i<m>]} and \spad{colList =
+       ++ [j<1>,j<2>,...,j<n>]}, then the \spad{(k,l)}th entry of
+       ++ \spad{elt(x,rowList,colList)} is \spad{x(i<k>,j<l>)}.
+     setelt: (%,List Integer,List Integer, %) -> %
+       ++ \spad{setelt(x,rowList,colList,y)} destructively alters the matrix x.
+       ++ If y is \spad{m}-by-\spad{n}, \spad{rowList = [i<1>,i<2>,...,i<m>]}
+       ++ and \spad{colList = [j<1>,j<2>,...,j<n>]}, then \spad{x(i<k>,j<l>)}
+       ++ is set to \spad{y(k,l)} for \spad{k = 1,...,m} and \spad{l = 1,...,n}.
+     swapRows_!: (%,Integer,Integer) -> %
+       ++ \spad{swapRows!(m,i,j)} interchanges the \spad{i}th and \spad{j}th
+       ++ rows of m. This destructively alters the matrix.
+     swapColumns_!: (%,Integer,Integer) -> %
+       ++ \spad{swapColumns!(m,i,j)} interchanges the \spad{i}th and \spad{j}th
+       ++ columns of m. This destructively alters the matrix.
+     subMatrix: (%,Integer,Integer,Integer,Integer) -> %
+       ++ \spad{subMatrix(x,i1,i2,j1,j2)} extracts the submatrix
+       ++ \spad{[x(i,j)]} where the index i ranges from \spad{i1} to \spad{i2}
+       ++ and the index j ranges from \spad{j1} to \spad{j2}.
+     setsubMatrix_!: (%,Integer,Integer,%) -> %
+       ++ \spad{setsubMatrix(x,i1,j1,y)} destructively alters the
+       ++ matrix x. Here \spad{x(i,j)} is set to \spad{y(i-i1+1,j-j1+1)} for
+       ++ \spad{i = i1,...,i1-1+nrows y} and \spad{j = j1,...,j1-1+ncols y}.
+
+--% Arithmetic
+
+     "+": (%,%) -> %
+       ++ \spad{x + y} is the sum of the matrices x and y.
+       ++ Error: if the dimensions are incompatible.
+     "-": (%,%) -> %
+       ++ \spad{x - y} is the difference of the matrices x and y.
+       ++ Error: if the dimensions are incompatible.
+     "-":  %    -> %
+       ++ \spad{-x} returns the negative of the matrix x.
+     "*": (%,%) -> %
+       ++ \spad{x * y} is the product of the matrices x and y.
+       ++ Error: if the dimensions are incompatible.
+     "*": (R,%) -> %
+       ++ \spad{r*x} is the left scalar multiple of the scalar r and the
+       ++ matrix x.
+     "*": (%,R) -> %
+       ++ \spad{x * r} is the right scalar multiple of the scalar r and the
+       ++ matrix x.
+     "*": (Integer,%) -> %
+       ++ \spad{n * x} is an integer multiple.
+     "*": (%,Col) -> Col
+       ++ \spad{x * c} is the product of the matrix x and the column vector c.
+       ++ Error: if the dimensions are incompatible.
+     "*": (Row,%) -> Row
+       ++ \spad{r * x} is the product of the row vector r and the matrix x.
+       ++ Error: if the dimensions are incompatible.
+     "**": (%,NonNegativeInteger) -> %
+       ++ \spad{x ** n} computes a non-negative integral power of the matrix x.
+       ++ Error: if the matrix is not square.
+     if R has IntegralDomain then
+       "exquo": (%,R) -> Union(%,"failed")
+         ++ \spad{exquo(m,r)} computes the exact quotient of the elements
+         ++ of m by r, returning \axiom{"failed"} if this is not possible.
+     if R has Field then
+       "/": (%,R) -> %
+         ++ \spad{m/r} divides the elements of m by r. Error: if \spad{r = 0}.
+
+--% Linear algebra
+
+     if R has EuclideanDomain then
+       rowEchelon: % -> %
+         ++ \spad{rowEchelon(m)} returns the row echelon form of the matrix m.
+     if R has IntegralDomain then
+       rank: % -> NonNegativeInteger
+         ++ \spad{rank(m)} returns the rank of the matrix m.
+       nullity: % -> NonNegativeInteger
+         ++ \spad{nullity(m)} returns the nullity of the matrix m. This is
+         ++ the dimension of the null space of the matrix m.
+       nullSpace: % -> List Col
+         ++ \spad{nullSpace(m)} returns a basis for the null space of
+         ++ the matrix m.
+     if R has commutative("*") then
+       determinant: % -> R
+         ++ \spad{determinant(m)} returns the determinant of the matrix m.
+         ++ Error: if the matrix is not square.
+       minordet: % -> R
+         ++ \spad{minordet(m)} computes the determinant of the matrix m using
+         ++ minors. Error: if the matrix is not square.
+     if R has Field then
+       inverse: % -> Union(%,"failed")
+         ++ \spad{inverse(m)} returns the inverse of the matrix m.
+         ++ If the matrix is not invertible, "failed" is returned.
+         ++ Error: if the matrix is not square.
+       "**": (%,Integer) -> %
+         ++ \spad{m**n} computes an integral power of the matrix m.
+         ++ Error: if matrix is not square or if the matrix
+         ++ is square but not invertible.
+
+   add
+     minr ==> minRowIndex
+     maxr ==> maxRowIndex
+     minc ==> minColIndex
+     maxc ==> maxColIndex
+     mini ==> minIndex
+     maxi ==> maxIndex
+
+--% Predicates
+
+     square? x == nrows x = ncols x
+
+     diagonal? x ==
+       not square? x => false
+       for i in minr x .. maxr x repeat
+         for j in minc x .. maxc x | (j - minc x) ^= (i - minr x) repeat
+           not zero? qelt(x, i, j) => return false
+       true
+
+     symmetric? x ==
+       (nRows := nrows x) ^= ncols x => false
+       mr := minRowIndex x; mc := minColIndex x
+       for i in 0..(nRows - 1) repeat
+         for j in (i + 1)..(nRows - 1) repeat
+           qelt(x,mr + i,mc + j) ^= qelt(x,mr + j,mc + i) => return false
+       true
+
+     antisymmetric? x ==
+       (nRows := nrows x) ^= ncols x => false
+       mr := minRowIndex x; mc := minColIndex x
+       for i in 0..(nRows - 1) repeat
+         for j in i..(nRows - 1) repeat
+           qelt(x,mr + i,mc + j) ^= -qelt(x,mr + j,mc + i) =>
+             return false
+       true
+
+--% Creation of matrices
+
+     zero(rows,cols) == new(rows,cols,0)
+
+     matrix(l: List List R) ==
+       null l => new(0,0,0)
+       -- error check: this is a top level function
+       rows : NonNegativeInteger := 1; cols := # first l
+       cols = 0 => error "matrices with zero columns are not supported"
+       for ll in rest l repeat
+         cols ^= # ll => error "matrix: rows of different lengths"
+         rows := rows + 1
+       ans := new(rows,cols,0)
+       for i in minr(ans)..maxr(ans) for ll in l repeat
+         for j in minc(ans)..maxc(ans) for r in ll repeat
+           qsetelt_!(ans,i,j,r)
+       ans
+
+     scalarMatrix(n,r) ==
+       ans := zero(n,n)
+       for i in minr(ans)..maxr(ans) for j in minc(ans)..maxc(ans) repeat
+         qsetelt_!(ans,i,j,r)
+       ans
+
+     diagonalMatrix(l: List R) ==
+       n := #l; ans := zero(n,n)
+       for i in minr(ans)..maxr(ans) for j in minc(ans)..maxc(ans) _
+           for r in l repeat qsetelt_!(ans,i,j,r)
+       ans
+
+     diagonalMatrix(list: List %) ==
+       rows : NonNegativeInteger := 0
+       cols : NonNegativeInteger := 0
+       for mat in list repeat
+         rows := rows + nrows mat
+         cols := cols + ncols mat
+       ans := zero(rows,cols)
+       loR := minr ans; loC := minc ans
+       for mat in list repeat
+         hiR := loR + nrows(mat) - 1; hiC := loC + nrows(mat) - 1
+         for i in loR..hiR for k in minr(mat)..maxr(mat) repeat
+           for j in loC..hiC for l in minc(mat)..maxc(mat) repeat
+             qsetelt_!(ans,i,j,qelt(mat,k,l))
+         loR := hiR + 1; loC := hiC + 1
+       ans
+
+     coerce(v:Col) ==
+       x := new(#v,1,0)
+       one := minc(x)
+       for i in minr(x)..maxr(x) for k in mini(v)..maxi(v) repeat
+         qsetelt_!(x,i,one,qelt(v,k))
+       x
+
+     transpose(v:Row) ==
+       x := new(1,#v,0)
+       one := minr(x)
+       for j in minc(x)..maxc(x) for k in mini(v)..maxi(v) repeat
+         qsetelt_!(x,one,j,qelt(v,k))
+       x
+
+     transpose(x:%) ==
+       ans := new(ncols x,nrows x,0)
+       for i in minr(ans)..maxr(ans) repeat
+         for j in minc(ans)..maxc(ans) repeat
+           qsetelt_!(ans,i,j,qelt(x,j,i))
+       ans
+
+     squareTop x ==
+       nrows x < (cols := ncols x) =>
+         error "squareTop: number of columns exceeds number of rows"
+       ans := new(cols,cols,0)
+       for i in minr(x)..(minr(x) + cols - 1) repeat
+         for j in minc(x)..maxc(x) repeat
+           qsetelt_!(ans,i,j,qelt(x,i,j))
+       ans
+
+     horizConcat(x,y) ==
+       (rows := nrows x) ^= nrows y =>
+         error "HConcat: matrices must have same number of rows"
+       ans := new(rows,(cols := ncols x) + ncols y,0)
+       for i in minr(x)..maxr(x) repeat
+         for j in minc(x)..maxc(x) repeat
+           qsetelt_!(ans,i,j,qelt(x,i,j))
+       for i in minr(y)..maxr(y) repeat
+         for j in minc(y)..maxc(y) repeat
+           qsetelt_!(ans,i,j + cols,qelt(y,i,j))
+       ans
+
+     vertConcat(x,y) ==
+       (cols := ncols x) ^= ncols y =>
+         error "HConcat: matrices must have same number of columns"
+       ans := new((rows := nrows x) + nrows y,cols,0)
+       for i in minr(x)..maxr(x) repeat
+         for j in minc(x)..maxc(x) repeat
+           qsetelt_!(ans,i,j,qelt(x,i,j))
+       for i in minr(y)..maxr(y) repeat
+         for j in minc(y)..maxc(y) repeat
+           qsetelt_!(ans,i + rows,j,qelt(y,i,j))
+       ans
+
+--% Part extraction/assignment
+
+     listOfLists x ==
+       ll : List List R := nil()
+       for i in maxr(x)..minr(x) by -1 repeat
+         l : List R := nil()
+         for j in maxc(x)..minc(x) by -1 repeat
+           l := cons(qelt(x,i,j),l)
+         ll := cons(l,ll)
+       ll
+
+     swapRows_!(x,i1,i2) ==
+       (i1 < minr(x)) or (i1 > maxr(x)) or (i2 < minr(x)) or _
+           (i2 > maxr(x)) => error "swapRows!: index out of range"
+       i1 = i2 => x
+       for j in minc(x)..maxc(x) repeat
+         r := qelt(x,i1,j)
+         qsetelt_!(x,i1,j,qelt(x,i2,j))
+         qsetelt_!(x,i2,j,r)
+       x
+
+     swapColumns_!(x,j1,j2) ==
+       (j1 < minc(x)) or (j1 > maxc(x)) or (j2 < minc(x)) or _
+           (j2 > maxc(x)) => error "swapColumns!: index out of range"
+       j1 = j2 => x
+       for i in minr(x)..maxr(x) repeat
+         r := qelt(x,i,j1)
+         qsetelt_!(x,i,j1,qelt(x,i,j2))
+         qsetelt_!(x,i,j2,r)
+       x
+
+     elt(x:%,rowList:List Integer,colList:List Integer) ==
+       for ei in rowList repeat
+         (ei < minr(x)) or (ei > maxr(x)) =>
+           error "elt: index out of range"
+       for ej in colList repeat
+         (ej < minc(x)) or (ej > maxc(x)) =>
+           error "elt: index out of range"
+       y := new(# rowList,# colList,0)
+       for ei in rowList for i in minr(y)..maxr(y) repeat
+         for ej in colList for j in minc(y)..maxc(y) repeat
+           qsetelt_!(y,i,j,qelt(x,ei,ej))
+       y
+
+     setelt(x:%,rowList:List Integer,colList:List Integer,y:%) ==
+       for ei in rowList repeat
+         (ei < minr(x)) or (ei > maxr(x)) =>
+           error "setelt: index out of range"
+       for ej in colList repeat
+         (ej < minc(x)) or (ej > maxc(x)) =>
+           error "setelt: index out of range"
+       ((# rowList) ^= (nrows y)) or ((# colList) ^= (ncols y)) =>
+         error "setelt: matrix has bad dimensions"
+       for ei in rowList for i in minr(y)..maxr(y) repeat
+         for ej in colList for j in minc(y)..maxc(y) repeat
+           qsetelt_!(x,ei,ej,qelt(y,i,j))
+       y
+
+     subMatrix(x,i1,i2,j1,j2) ==
+       (i2 < i1) => error "subMatrix: bad row indices"
+       (j2 < j1) => error "subMatrix: bad column indices"
+       (i1 < minr(x)) or (i2 > maxr(x)) =>
+         error "subMatrix: index out of range"
+       (j1 < minc(x)) or (j2 > maxc(x)) =>
+         error "subMatrix: index out of range"
+       rows := (i2 - i1 + 1) pretend NonNegativeInteger
+       cols := (j2 - j1 + 1) pretend NonNegativeInteger
+       y := new(rows,cols,0)
+       for i in minr(y)..maxr(y) for k in i1..i2 repeat
+         for j in minc(y)..maxc(y) for l in j1..j2 repeat
+           qsetelt_!(y,i,j,qelt(x,k,l))
+       y
+
+     setsubMatrix_!(x,i1,j1,y) ==
+       i2 := i1 + nrows(y) -1
+       j2 := j1 + ncols(y) -1
+       (i1 < minr(x)) or (i2 > maxr(x)) =>
+         error "setsubMatrix!: inserted matrix too big, use subMatrix to restrict it"
+       (j1 < minc(x)) or (j2 > maxc(x)) =>
+         error "setsubMatrix!: inserted matrix too big, use subMatrix to restrict it"
+       for i in minr(y)..maxr(y) for k in i1..i2 repeat
+         for j in minc(y)..maxc(y) for l in j1..j2 repeat
+           qsetelt_!(x,k,l,qelt(y,i,j))
+       x
+
+--% Arithmetic
+
+     x + y ==
+       ((r := nrows x) ^= nrows y) or ((c := ncols x) ^= ncols y) =>
+         error "can't add matrices of different dimensions"
+       ans := new(r,c,0)
+       for i in minr(x)..maxr(x) repeat
+         for j in minc(x)..maxc(x) repeat
+           qsetelt_!(ans,i,j,qelt(x,i,j) + qelt(y,i,j))
+       ans
+
+     x - y ==
+       ((r := nrows x) ^= nrows y) or ((c := ncols x) ^= ncols y) =>
+         error "can't subtract matrices of different dimensions"
+       ans := new(r,c,0)
+       for i in minr(x)..maxr(x) repeat
+         for j in minc(x)..maxc(x) repeat
+           qsetelt_!(ans,i,j,qelt(x,i,j) - qelt(y,i,j))
+       ans
+
+     - x == map(- #1,x)
+
+     a:R * x:% == map(a * #1,x)
+     x:% * a:R == map(#1 * a,x)
+     m:Integer * x:% == map(m * #1,x)
+
+     x:% * y:% ==
+       (ncols x ^= nrows y) =>
+         error "can't multiply matrices of incompatible dimensions"
+       ans := new(nrows x,ncols y,0)
+       for i in minr(x)..maxr(x) repeat
+         for j in minc(y)..maxc(y) repeat
+           entry :=
+             sum : R := 0
+             for k in minr(y)..maxr(y) for l in minc(x)..maxc(x) repeat
+               sum := sum + qelt(x,i,l) * qelt(y,k,j)
+             sum
+           qsetelt_!(ans,i,j,entry)
+       ans
+
+     positivePower:(%,Integer) -> %
+     positivePower(x,n) ==
+--       one? n => x
+       (n = 1) => x
+       odd? n => x * positivePower(x,n - 1)
+       y := positivePower(x,n quo 2)
+       y * y
+
+     x:% ** n:NonNegativeInteger ==
+       not((nn:= nrows x) = ncols x) => error "**: matrix must be square"
+       zero? n => scalarMatrix(nn,1)
+       positivePower(x,n)
+
+     --if R has ConvertibleTo InputForm then
+       --convert(x:%):InputForm ==
+         --convert [convert("matrix"::Symbol)@InputForm,
+                  --convert listOfLists x]$List(InputForm)
+
+     if Col has shallowlyMutable then
+
+       x:% * v:Col ==
+         ncols(x) ^= #v =>
+           error "can't multiply matrix A and vector v if #cols A ^= #v"
+         w : Col := new(nrows x,0)
+         for i in minr(x)..maxr(x) for k in mini(w)..maxi(w) repeat
+           w.k :=
+             sum : R := 0
+             for j in minc(x)..maxc(x) for l in mini(v)..maxi(v) repeat
+               sum := sum + qelt(x,i,j) * v(l)
+             sum
+         w
+
+     if Row has shallowlyMutable then
+
+       v:Row * x:% ==
+         nrows(x) ^= #v =>
+           error "can't multiply vector v and matrix A if #rows A ^= #v"
+         w : Row := new(ncols x,0)
+         for j in minc(x)..maxc(x) for k in mini(w)..maxi(w) repeat
+           w.k :=
+             sum : R := 0
+             for i in minr(x)..maxr(x) for l in mini(v)..maxi(v) repeat
+               sum := sum + qelt(x,i,j) * v(l)
+             sum
+         w
+
+     if R has IntegralDomain then
+       x exquo a ==
+         ans := new(nrows x,ncols x,0)
+         for i in minr(x)..maxr(x) repeat
+           for j in minc(x)..maxc(x) repeat
+             entry :=
+               (r := (qelt(x,i,j) exquo a)) case "failed" =>
+                 return "failed"
+               r :: R
+             qsetelt_!(ans,i,j,entry)
+         ans
+
+     if R has Field then
+       x / r == map(#1 / r,x)
+
+       x:% ** n:Integer ==
+         not((nn:= nrows x) = ncols x) => error "**: matrix must be square"
+         zero? n => scalarMatrix(nn,1)
+         positive? n => positivePower(x,n)
+         (xInv := inverse x) case "failed" =>
+           error "**: matrix must be invertible"
+         positivePower(xInv :: %,-n)
+
+@
+\section{category RMATCAT RectangularMatrixCategory}
+<<category RMATCAT RectangularMatrixCategory>>=
+)abbrev category RMATCAT RectangularMatrixCategory
+++ Authors: Grabmeier, Gschnitzer, Williamson
+++ Date Created: 1987
+++ Date Last Updated: July 1990
+++ Basic Operations:
+++ Related Domains: RectangularMatrix(m,n,R)
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ Examples:
+++ References:
+++ Description:
+++   \spadtype{RectangularMatrixCategory} is a category of matrices of fixed
+++   dimensions.  The dimensions of the matrix will be parameters of the
+++   domain.  Domains in this category will be R-modules and will be
+++   non-mutable.
+RectangularMatrixCategory(m,n,R,Row,Col): Category == Definition where
+  m,n : NonNegativeInteger
+  R   : Ring
+  Row : DirectProductCategory(n,R)
+  Col : DirectProductCategory(m,R)
+
+  Definition ==> Join(BiModule(R,R),HomogeneousAggregate(R)) with
+
+    finiteAggregate
+      ++ matrices are finite
+
+    if R has CommutativeRing then Module(R)
+
+--% Matrix creation
+
+    matrix: List List R -> %
+      ++ \spad{matrix(l)} converts the list of lists l to a matrix, where the
+      ++ list of lists is viewed as a list of the rows of the matrix.
+
+--% Predicates
+
+    square?  : % -> Boolean
+      ++ \spad{square?(m)} returns true if m is a square matrix (i.e. if m
+      ++ has the same number of rows as columns) and false otherwise.
+    diagonal?: % -> Boolean
+      ++ \spad{diagonal?(m)} returns true if the matrix m is square and diagonal
+      ++ (i.e. all entries of m not on the diagonal are zero) and false
+      ++ otherwise.
+    symmetric?: % -> Boolean
+      ++ \spad{symmetric?(m)} returns true if the matrix m is square and
+      ++ symmetric (i.e. \spad{m[i,j] = m[j,i]} for all \spad{i} and j) and
+      ++ false otherwise.
+    antisymmetric?: % -> Boolean
+      ++ \spad{antisymmetric?(m)} returns true if the matrix m is square and
+      ++ antisymmetric (i.e. \spad{m[i,j] = -m[j,i]} for all \spad{i} and j)
+      ++ and false otherwise.
+
+--% Size inquiries
+
+    minRowIndex : % -> Integer
+      ++ \spad{minRowIndex(m)} returns the index of the 'first' row of the
+      ++ matrix m.
+    maxRowIndex : % -> Integer
+      ++ \spad{maxRowIndex(m)} returns the index of the 'last' row of the
+      ++ matrix m.
+    minColIndex : % -> Integer
+      ++ \spad{minColIndex(m)} returns the index of the 'first' column of the
+      ++ matrix m.
+    maxColIndex : % -> Integer
+      ++ \spad{maxColIndex(m)} returns the index of the 'last' column of the
+      ++ matrix m.
+    nrows : % -> NonNegativeInteger
+      ++ \spad{nrows(m)} returns the number of rows in the matrix m.
+    ncols : % -> NonNegativeInteger
+      ++ \spad{ncols(m)} returns the number of columns in the matrix m.
+
+--% Part extractions
+
+    listOfLists: % -> List List R
+      ++ \spad{listOfLists(m)} returns the rows of the matrix m as a list
+      ++ of lists.
+    elt: (%,Integer,Integer) -> R
+      ++ \spad{elt(m,i,j)} returns the element in the \spad{i}th row and
+      ++ \spad{j}th column of the matrix m.
+      ++ Error: if indices are outside the proper
+      ++ ranges.
+    qelt: (%,Integer,Integer) -> R
+      ++ \spad{qelt(m,i,j)} returns the element in the \spad{i}th row and
+      ++ \spad{j}th column of
+      ++ the matrix m. Note: there is NO error check to determine if indices are
+      ++ in the proper ranges.
+    elt: (%,Integer,Integer,R) -> R
+      ++ \spad{elt(m,i,j,r)} returns the element in the \spad{i}th row and
+      ++ \spad{j}th column of the matrix m, if m has an \spad{i}th row and a
+      ++ \spad{j}th column, and returns r otherwise.
+    row: (%,Integer) -> Row
+      ++ \spad{row(m,i)} returns the \spad{i}th row of the matrix m.
+      ++ Error: if the index is outside the proper range.
+    column: (%,Integer) -> Col
+      ++ \spad{column(m,j)} returns the \spad{j}th column of the matrix m.
+      ++ Error: if the index outside the proper range.
+
+--% Map and Zip
+
+    map: (R -> R,%) -> %
+      ++ \spad{map(f,a)} returns b, where \spad{b(i,j) = a(i,j)} for all i, j.
+    map:((R,R) -> R,%,%) -> %
+      ++ \spad{map(f,a,b)} returns c, where c is such that
+      ++ \spad{c(i,j) = f(a(i,j),b(i,j))} for all \spad{i}, j.
+
+--% Arithmetic
+
+    if R has IntegralDomain then
+      "exquo": (%,R) -> Union(%,"failed")
+        ++ \spad{exquo(m,r)} computes the exact quotient of the elements
+        ++ of m by r, returning \axiom{"failed"} if this is not possible.
+    if R has Field then
+      "/": (%,R) -> %
+        ++ \spad{m/r} divides the elements of m by r. Error: if \spad{r = 0}.
+
+--% Linear algebra
+
+    if R has EuclideanDomain then
+      rowEchelon: % -> %
+        ++ \spad{rowEchelon(m)} returns the row echelon form of the matrix m.
+    if R has IntegralDomain then
+      rank: % -> NonNegativeInteger
+        ++ \spad{rank(m)} returns the rank of the matrix m.
+      nullity: % -> NonNegativeInteger
+        ++ \spad{nullity(m)} returns the nullity of the matrix m. This is
+        ++ the dimension of the null space of the matrix m.
+      nullSpace: % -> List Col
+        ++ \spad{nullSpace(m)}+ returns a basis for the null space of
+        ++ the matrix m.
+   add
+     nrows x == m
+     ncols x == n
+     square? x == m = n
+
+     diagonal? x ==
+       not square? x => false
+       for i in minRowIndex x .. maxRowIndex x repeat
+         for j in minColIndex x .. maxColIndex x
+           | (j - minColIndex x) ^= (i - minRowIndex x) repeat
+             not zero? qelt(x, i, j) => return false
+       true
+
+     symmetric? x ==
+       m ^= n => false
+       mr := minRowIndex x; mc := minColIndex x
+       for i in 0..(n - 1) repeat
+         for j in (i + 1)..(n - 1) repeat
+           qelt(x,mr + i,mc + j) ^= qelt(x,mr + j,mc + i) => return false
+       true
+
+     antisymmetric? x ==
+       m ^= n => false
+       mr := minRowIndex x; mc := minColIndex x
+       for i in 0..(n - 1) repeat
+         for j in i..(n - 1) repeat
+           qelt(x,mr + i,mc + j) ^= -qelt(x,mr + j,mc + i) =>
+             return false
+       true
+
+@
+\section{category SMATCAT SquareMatrixCategory}
+<<category SMATCAT SquareMatrixCategory>>=
+)abbrev category SMATCAT SquareMatrixCategory
+++ Authors: Grabmeier, Gschnitzer, Williamson
+++ Date Created: 1987
+++ Date Last Updated: July 1990
+++ Basic Operations:
+++ Related Domains: SquareMatrix(ndim,R)
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ Examples:
+++ References:
+++ Description:
+++   \spadtype{SquareMatrixCategory} is a general square matrix category which
+++   allows different representations and indexing schemes.  Rows and
+++   columns may be extracted with rows returned as objects of
+++   type Row and colums returned as objects of type Col.
+SquareMatrixCategory(ndim,R,Row,Col): Category == Definition where
+  ndim : NonNegativeInteger
+  R    : Ring
+  Row  : DirectProductCategory(ndim,R)
+  Col  : DirectProductCategory(ndim,R)
+  I ==> Integer
+
+  Definition ==> Join(DifferentialExtension R, BiModule(R, R),_
+                      RectangularMatrixCategory(ndim,ndim,R,Row,Col),_
+                      FullyRetractableTo R,_
+                      FullyLinearlyExplicitRingOver R) with
+    if R has CommutativeRing then Module(R)
+    scalarMatrix: R -> %
+      ++ \spad{scalarMatrix(r)} returns an n-by-n matrix with r's on the
+      ++ diagonal and zeroes elsewhere.
+    diagonalMatrix: List R -> %
+      ++ \spad{diagonalMatrix(l)} returns a diagonal matrix with the elements
+      ++ of l on the diagonal.
+    diagonal: % -> Row
+      ++ \spad{diagonal(m)} returns a row consisting of the elements on the
+      ++ diagonal of the matrix m.
+    trace: % -> R
+      ++ \spad{trace(m)} returns the trace of the matrix m. this is the sum
+      ++ of the elements on the diagonal of the matrix m.
+    diagonalProduct: % -> R
+      ++ \spad{diagonalProduct(m)} returns the product of the elements on the
+      ++ diagonal of the matrix m.
+    "*": (%,Col) -> Col
+      ++ \spad{x * c} is the product of the matrix x and the column vector c.
+      ++ Error: if the dimensions are incompatible.
+    "*": (Row,%) -> Row
+      ++ \spad{r * x} is the product of the row vector r and the matrix x.
+      ++ Error: if the dimensions are incompatible.
+
+--% Linear algebra
+
+    if R has commutative("*") then
+      Algebra R
+      determinant: % -> R
+        ++ \spad{determinant(m)} returns the determinant of the matrix m.
+      minordet: % -> R
+        ++ \spad{minordet(m)} computes the determinant of the matrix m
+        ++ using minors.
+    if R has Field then
+      inverse: % -> Union(%,"failed")
+        ++ \spad{inverse(m)} returns the inverse of the matrix m, if that
+        ++ matrix is invertible and returns "failed" otherwise.
+      "**": (%,Integer) -> %
+        ++ \spad{m**n} computes an integral power of the matrix m.
+        ++ Error: if the matrix is not invertible.
+
+   add
+    minr ==> minRowIndex
+    maxr ==> maxRowIndex
+    minc ==> minColIndex
+    maxc ==> maxColIndex
+    mini ==> minIndex
+    maxi ==> maxIndex
+
+    positivePower:(%,Integer) -> %
+    positivePower(x,n) ==
+--      one? n => x
+      (n = 1) => x
+      odd? n => x * positivePower(x,n - 1)
+      y := positivePower(x,n quo 2)
+      y * y
+
+    x:% ** n:NonNegativeInteger ==
+      zero? n => scalarMatrix 1
+      positivePower(x,n)
+
+    coerce(r:R) == scalarMatrix r
+
+    equation2R: Vector % -> Matrix R
+
+    differentiate(x:%,d:R -> R) == map(d,x)
+
+    diagonal x ==
+      v:Vector(R) := new(ndim,0)
+      for i in minr x .. maxr x
+        for j in minc x .. maxc x
+          for k in minIndex v .. maxIndex v repeat
+            qsetelt_!(v, k, qelt(x, i, j))
+      directProduct v
+
+    retract(x:%):R ==
+      diagonal? x => retract diagonal x
+      error "Not retractable"
+
+    retractIfCan(x:%):Union(R, "failed") ==
+      diagonal? x => retractIfCan diagonal x
+      "failed"
+
+    equation2R v ==
+      ans:Matrix(Col) := new(ndim,#v,0)
+      for i in minr ans .. maxr ans repeat
+        for j in minc ans .. maxc ans repeat
+          qsetelt_!(ans, i, j, column(qelt(v, j), i))
+      reducedSystem ans
+
+    reducedSystem(x:Matrix %):Matrix(R) ==
+      empty? x => new(0,0,0)
+      reduce(vertConcat, [equation2R row(x, i)
+                               for i in minr x .. maxr x])$List(Matrix R)
+
+    reducedSystem(m:Matrix %, v:Vector %):
+     Record(mat:Matrix R, vec:Vector R) ==
+      vh:Vector(R) :=
+        empty? v => new(0,0)
+        rh := reducedSystem(v::Matrix %)@Matrix(R)
+        column(rh, minColIndex rh)
+      [reducedSystem(m)@Matrix(R), vh]
+
+    trace x ==
+      tr : R := 0
+      for i in minr(x)..maxr(x) for j in minc(x)..maxc(x) repeat
+        tr := tr + x(i,j)
+      tr
+
+    diagonalProduct x ==
+      pr : R := 1
+      for i in minr(x)..maxr(x) for j in minc(x)..maxc(x) repeat
+        pr := pr * x(i,j)
+      pr
+
+    if R has Field then
+
+      x:% ** n:Integer ==
+        zero? n => scalarMatrix 1
+        positive? n => positivePower(x,n)
+        (xInv := inverse x) case "failed" =>
+          error "**: matrix must be invertible"
+        positivePower(xInv :: %,-n)
+
+@
+\section{License}
+<<license>>=
+--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
+--All rights reserved.
+--
+--Redistribution and use in source and binary forms, with or without
+--modification, are permitted provided that the following conditions are
+--met:
+--
+--    - Redistributions of source code must retain the above copyright
+--      notice, this list of conditions and the following disclaimer.
+--
+--    - Redistributions in binary form must reproduce the above copyright
+--      notice, this list of conditions and the following disclaimer in
+--      the documentation and/or other materials provided with the
+--      distribution.
+--
+--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
+--      names of its contributors may be used to endorse or promote products
+--      derived from this software without specific prior written permission.
+--
+--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+@
+<<*>>=
+<<license>>
+
+<<category MATCAT MatrixCategory>>
+<<category RMATCAT RectangularMatrixCategory>>
+<<category SMATCAT SquareMatrixCategory>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}