Initial population.

1 files changed, 530 insertions, 0 deletions

diff --git a/src/algebra/matrix.spad.pamphlet b/src/algebra/matrix.spad.pamphlet
new file mode 100644
index 00000000..3704e00e
--- /dev/null
+++ b/src/algebra/matrix.spad.pamphlet
@@ -0,0 +1,530 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra matrix.spad}
+\author{Johannes Grabmeier, Oswald Gschnitzer, Clifton J. Williamson}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{domain IMATRIX IndexedMatrix}
+<<domain IMATRIX IndexedMatrix>>=
+)abbrev domain IMATRIX IndexedMatrix
+++ Author: Grabmeier, Gschnitzer, Williamson
+++ Date Created: 1987
+++ Date Last Updated: July 1990
+++ Basic Operations:
+++ Related Domains: Matrix, RectangularMatrix, SquareMatrix,
+++   StorageEfficientMatrixOperations
+++ Also See:
+++ AMS Classifications:
+++ Keywords: matrix, linear algebra
+++ Examples:
+++ References:
+++ Description:
+++   An \spad{IndexedMatrix} is a matrix where the minimal row and column
+++   indices are parameters of the type.  The domains Row and Col
+++   are both IndexedVectors.
+++   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.
+IndexedMatrix(R,mnRow,mnCol): Exports == Implementation where
+  R : Ring
+  mnRow, mnCol : Integer
+  Row ==> IndexedVector(R,mnCol)
+  Col ==> IndexedVector(R,mnRow)
+  MATLIN ==> MatrixLinearAlgebraFunctions(R,Row,Col,$)
+ 
+  Exports ==> MatrixCategory(R,Row,Col)
+ 
+  Implementation ==>
+    InnerIndexedTwoDimensionalArray(R,mnRow,mnCol,Row,Col) add
+ 
+      swapRows_!(x,i1,i2) ==
+        (i1 < minRowIndex(x)) or (i1 > maxRowIndex(x)) or _
+           (i2 < minRowIndex(x)) or (i2 > maxRowIndex(x)) =>
+             error "swapRows!: index out of range"
+        i1 = i2 => x
+        minRow := minRowIndex x
+        xx := x pretend PrimitiveArray(PrimitiveArray(R))
+        n1 := i1 - minRow; n2 := i2 - minRow
+        row1 := qelt(xx,n1)
+        qsetelt_!(xx,n1,qelt(xx,n2))
+        qsetelt_!(xx,n2,row1)
+        xx pretend $
+ 
+      if R has commutative("*") then
+ 
+        determinant x == determinant(x)$MATLIN
+        minordet    x == minordet(x)$MATLIN
+ 
+      if R has EuclideanDomain then
+ 
+        rowEchelon  x == rowEchelon(x)$MATLIN
+ 
+      if R has IntegralDomain then
+ 
+        rank        x == rank(x)$MATLIN
+        nullity     x == nullity(x)$MATLIN
+        nullSpace   x == nullSpace(x)$MATLIN
+ 
+      if R has Field then
+ 
+        inverse     x == inverse(x)$MATLIN
+
+@
+\section{domain MATRIX Matrix}
+<<domain MATRIX Matrix>>=
+)abbrev domain MATRIX Matrix
+++ Author: Grabmeier, Gschnitzer, Williamson
+++ Date Created: 1987
+++ Date Last Updated: July 1990
+++ Basic Operations:
+++ Related Domains: IndexedMatrix, RectangularMatrix, SquareMatrix
+++ Also See:
+++ AMS Classifications:
+++ Keywords: matrix, linear algebra
+++ Examples:
+++ References:
+++ Description:
+++   \spadtype{Matrix} is a matrix domain where 1-based indexing is used
+++   for both rows and columns.
+Matrix(R): Exports == Implementation where
+  R : Ring
+  Row ==> Vector R
+  Col ==> Vector R
+  mnRow ==> 1
+  mnCol ==> 1
+  MATLIN ==> MatrixLinearAlgebraFunctions(R,Row,Col,$)
+  MATSTOR ==> StorageEfficientMatrixOperations(R)
+ 
+  Exports ==> MatrixCategory(R,Row,Col) with
+    diagonalMatrix: Vector R -> $
+      ++ \spad{diagonalMatrix(v)} returns a diagonal matrix where the elements
+      ++ of v appear on the diagonal.
+
+    if R has ConvertibleTo InputForm then ConvertibleTo InputForm
+
+    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.
+--     matrix: Vector Vector R -> $
+--       ++ \spad{matrix(v)} converts the vector of vectors v to a matrix, where
+--       ++ the vector of vectors is viewed as a vector of the rows of the
+--       ++ matrix
+--     diagonalMatrix: Vector $ -> $
+--       ++ \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.
+--     vectorOfVectors: $ -> Vector Vector R
+--       ++ \spad{vectorOfVectors(m)} returns the rows of the matrix m as a
+--       ++ vector of vectors
+ 
+  Implementation ==>
+   InnerIndexedTwoDimensionalArray(R,mnRow,mnCol,Row,Col) add
+    minr ==> minRowIndex
+    maxr ==> maxRowIndex
+    minc ==> minColIndex
+    maxc ==> maxColIndex
+    mini ==> minIndex
+    maxi ==> maxIndex
+ 
+    minRowIndex x == mnRow
+    minColIndex x == mnCol
+ 
+    swapRows_!(x,i1,i2) ==
+        (i1 < minRowIndex(x)) or (i1 > maxRowIndex(x)) or _
+           (i2 < minRowIndex(x)) or (i2 > maxRowIndex(x)) =>
+             error "swapRows!: index out of range"
+        i1 = i2 => x
+        minRow := minRowIndex x
+        xx := x pretend PrimitiveArray(PrimitiveArray(R))
+        n1 := i1 - minRow; n2 := i2 - minRow
+        row1 := qelt(xx,n1)
+        qsetelt_!(xx,n1,qelt(xx,n2))
+        qsetelt_!(xx,n2,row1)
+        xx pretend $
+ 
+    positivePower:($,Integer,NonNegativeInteger) -> $
+    positivePower(x,n,nn) ==
+--      one? n => x
+      (n = 1) => x
+      -- no need to allocate space for 3 additional matrices
+      n = 2 => x * x
+      n = 3 => x * x * x
+      n = 4 => (y := x * x; y * y)
+      a := new(nn,nn,0) pretend Matrix(R)
+      b := new(nn,nn,0) pretend Matrix(R)
+      c := new(nn,nn,0) pretend Matrix(R)
+      xx := x pretend Matrix(R)
+      power_!(a,b,c,xx,n :: NonNegativeInteger)$MATSTOR pretend $
+ 
+    x:$ ** n:NonNegativeInteger ==
+      not((nn := nrows x) = ncols x) =>
+        error "**: matrix must be square"
+      zero? n => scalarMatrix(nn,1)
+      positivePower(x,n,nn)
+ 
+    if R has commutative("*") then
+ 
+        determinant x == determinant(x)$MATLIN
+        minordet    x == minordet(x)$MATLIN
+ 
+    if R has EuclideanDomain then
+ 
+        rowEchelon  x == rowEchelon(x)$MATLIN
+ 
+    if R has IntegralDomain then
+ 
+        rank        x == rank(x)$MATLIN
+        nullity     x == nullity(x)$MATLIN
+        nullSpace   x == nullSpace(x)$MATLIN
+ 
+    if R has Field then
+ 
+        inverse     x == inverse(x)$MATLIN
+ 
+        x:$ ** n:Integer ==
+          nn := nrows x
+          not(nn = ncols x) =>
+            error "**: matrix must be square"
+          zero? n => scalarMatrix(nn,1)
+          positive? n => positivePower(x,n,nn)
+          (xInv := inverse x) case "failed" =>
+            error "**: matrix must be invertible"
+          positivePower(xInv :: $,-n,nn)
+ 
+--     matrix(v: Vector Vector R) ==
+--       (rows := # v) = 0 => new(0,0,0)
+--       -- error check: this is a top level function
+--       cols := # v.mini(v)
+--       for k in (mini(v) + 1)..maxi(v) repeat
+--         cols ^= # v.k => error "matrix: rows of different lengths"
+--       ans := new(rows,cols,0)
+--       for i in minr(ans)..maxr(ans) for k in mini(v)..maxi(v) repeat
+--         vv := v.k
+--         for j in minc(ans)..maxc(ans) for l in mini(vv)..maxi(vv) repeat
+--           ans(i,j) := vv.l
+--       ans
+ 
+    diagonalMatrix(v: Vector R) ==
+      n := #v; ans := zero(n,n)
+      for i in minr(ans)..maxr(ans) for j in minc(ans)..maxc(ans) _
+          for k in mini(v)..maxi(v) repeat qsetelt_!(ans,i,j,qelt(v,k))
+      ans
+ 
+--     diagonalMatrix(vec: Vector $) ==
+--       rows : NonNegativeInteger := 0
+--       cols : NonNegativeInteger := 0
+--       for r in mini(vec)..maxi(vec) repeat
+--         mat := vec.r
+--         rows := rows + nrows mat; cols := cols + ncols mat
+--       ans := zero(rows,cols)
+--       loR := minr ans; loC := minc ans
+--       for r in mini(vec)..maxi(vec) repeat
+--         mat := vec.r
+--         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
+--             ans(i,j) := mat(k,l)
+--         loR := hiR + 1; loC := hiC + 1
+--       ans
+ 
+--     vectorOfVectors x ==
+--       vv : Vector Vector R := new(nrows x,0)
+--       cols := ncols x
+--       for k in mini(vv)..maxi(vv) repeat
+--         vv.k := new(cols,0)
+--       for i in minr(x)..maxr(x) for k in mini(vv)..maxi(vv) repeat
+--         v := vv.k
+--         for j in minc(x)..maxc(x) for l in mini(v)..maxi(v) repeat
+--           v.l := x(i,j)
+--       vv
+ 
+    if R has ConvertibleTo InputForm then
+      convert(x:$):InputForm ==
+         convert [convert("matrix"::Symbol)@InputForm,
+                  convert listOfLists x]$List(InputForm)
+
+@
+\section{domain RMATRIX RectangularMatrix}
+<<domain RMATRIX RectangularMatrix>>=
+)abbrev domain RMATRIX RectangularMatrix
+++ Author: Grabmeier, Gschnitzer, Williamson
+++ Date Created: 1987
+++ Date Last Updated: July 1990
+++ Basic Operations:
+++ Related Domains: IndexedMatrix, Matrix, SquareMatrix
+++ Also See:
+++ AMS Classifications:
+++ Keywords: matrix, linear algebra
+++ Examples:
+++ References:
+++ Description:
+++   \spadtype{RectangularMatrix} is a matrix domain where the number of rows
+++   and the number of columns are parameters of the domain.
+RectangularMatrix(m,n,R): Exports == Implementation where
+  m,n : NonNegativeInteger
+  R   : Ring
+  Row ==> DirectProduct(n,R)
+  Col ==> DirectProduct(m,R)
+  Exports ==> Join(RectangularMatrixCategory(m,n,R,Row,Col),_
+                   CoercibleTo Matrix R) with
+ 
+    if R has Field then VectorSpace R
+ 
+    if R has ConvertibleTo InputForm then ConvertibleTo InputForm
+
+    rectangularMatrix: Matrix R -> $
+      ++ \spad{rectangularMatrix(m)} converts a matrix of type \spadtype{Matrix}
+      ++ to a matrix of type \spad{RectangularMatrix}.
+    coerce: $ -> Matrix R
+      ++ \spad{coerce(m)} converts a matrix of type \spadtype{RectangularMatrix}
+      ++ to a matrix of type \spad{Matrix}.
+ 
+  Implementation ==> Matrix R add
+    minr ==> minRowIndex
+    maxr ==> maxRowIndex
+    minc ==> minColIndex
+    maxc ==> maxColIndex
+    mini ==> minIndex
+    maxi ==> maxIndex
+ 
+    ZERO := new(m,n,0)$Matrix(R) pretend $
+    0    == ZERO
+ 
+    coerce(x:$):OutputForm == coerce(x pretend Matrix R)$Matrix(R)
+
+    matrix(l: List List R) ==
+      -- error check: this is a top level function
+      #l ^= m => error "matrix: wrong number of rows"
+      for ll in l repeat
+        #ll ^= n => error "matrix: wrong number of columns"
+      ans : Matrix R := new(m,n,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 pretend $
+ 
+    row(x,i)    == directProduct row(x pretend Matrix(R),i)
+    column(x,j) == directProduct column(x pretend Matrix(R),j)
+ 
+    coerce(x:$):Matrix(R) == copy(x pretend Matrix(R))
+ 
+    rectangularMatrix x ==
+      (nrows(x) ^= m) or (ncols(x) ^= n) =>
+        error "rectangularMatrix: matrix of bad dimensions"
+      copy(x) pretend $
+ 
+    if R has EuclideanDomain then
+ 
+      rowEchelon x == rowEchelon(x pretend Matrix(R)) pretend $
+ 
+    if R has IntegralDomain then
+ 
+      rank x    == rank(x pretend Matrix(R))
+      nullity x == nullity(x pretend Matrix(R))
+      nullSpace x ==
+        [directProduct c for c in nullSpace(x pretend Matrix(R))]
+ 
+    if R has Field then
+ 
+      dimension() == (m * n) :: CardinalNumber
+ 
+    if R has ConvertibleTo InputForm then
+      convert(x:$):InputForm ==
+         convert [convert("rectangularMatrix"::Symbol)@InputForm,
+                  convert(x::Matrix(R))]$List(InputForm)
+
+@
+\section{domain SQMATRIX SquareMatrix}
+<<domain SQMATRIX SquareMatrix>>=
+)abbrev domain SQMATRIX SquareMatrix
+++ Author: Grabmeier, Gschnitzer, Williamson
+++ Date Created: 1987
+++ Date Last Updated: July 1990
+++ Basic Operations:
+++ Related Domains: IndexedMatrix, Matrix, RectangularMatrix
+++ Also See:
+++ AMS Classifications:
+++ Keywords: matrix, linear algebra
+++ Examples:
+++ References:
+++ Description:
+++   \spadtype{SquareMatrix} is a matrix domain of square matrices, where the
+++   number of rows (= number of columns) is a parameter of the type.
+SquareMatrix(ndim,R): Exports == Implementation where
+  ndim : NonNegativeInteger
+  R    : Ring
+  Row ==> DirectProduct(ndim,R)
+  Col ==> DirectProduct(ndim,R)
+  MATLIN ==> MatrixLinearAlgebraFunctions(R,Row,Col,$)
+ 
+  Exports ==> Join(SquareMatrixCategory(ndim,R,Row,Col),_
+                   CoercibleTo Matrix R) with
+ 
+    transpose: $ -> $
+      ++ \spad{transpose(m)} returns the transpose of the matrix m.
+    squareMatrix: Matrix R -> $
+      ++ \spad{squareMatrix(m)} converts a matrix of type \spadtype{Matrix}
+      ++ to a matrix of type \spadtype{SquareMatrix}.
+    coerce: $ -> Matrix R
+      ++ \spad{coerce(m)} converts a matrix of type \spadtype{SquareMatrix}
+      ++ to a matrix of type \spadtype{Matrix}.
+--  symdecomp : $ -> Record(sym:$,antisym:$)
+--    ++ \spad{symdecomp(m)} decomposes the matrix m as a sum of a symmetric
+--    ++ matrix \spad{m1} and an antisymmetric matrix \spad{m2}. The object
+--    ++ returned is the Record \spad{[m1,m2]}
+--  if R has commutative("*") then
+--    minorsVect: -> Vector(Union(R,"uncomputed")) --range: 1..2**n-1
+--      ++ \spad{minorsVect(m)} returns a vector of the minors of the matrix m
+    if R has commutative("*") then central
+      ++ the elements of the Ring R, viewed as diagonal matrices, commute
+      ++ with all matrices and, indeed, are the only matrices which commute
+      ++ with all matrices.
+    if R has commutative("*") and R has unitsKnown then unitsKnown
+      ++ the invertible matrices are simply the matrices whose determinants
+      ++ are units in the Ring R.
+    if R has ConvertibleTo InputForm then ConvertibleTo InputForm
+ 
+  Implementation ==> Matrix R add
+    minr ==> minRowIndex
+    maxr ==> maxRowIndex
+    minc ==> minColIndex
+    maxc ==> maxColIndex
+    mini ==> minIndex
+    maxi ==> maxIndex
+ 
+    ZERO := scalarMatrix 0
+    0    == ZERO
+    ONE  := scalarMatrix 1
+    1    == ONE
+
+    characteristic() == characteristic()$R
+ 
+    matrix(l: List List R) ==
+      -- error check: this is a top level function
+      #l ^= ndim => error "matrix: wrong number of rows"
+      for ll in l repeat
+        #ll ^= ndim => error "matrix: wrong number of columns"
+      ans : Matrix R := new(ndim,ndim,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 pretend $
+ 
+    row(x,i)    == directProduct row(x pretend Matrix(R),i)
+    column(x,j) == directProduct column(x pretend Matrix(R),j)
+    coerce(x:$):OutputForm == coerce(x pretend Matrix R)$Matrix(R)
+ 
+    scalarMatrix r == scalarMatrix(ndim,r)$Matrix(R) pretend $
+ 
+    diagonalMatrix l ==
+      #l ^= ndim =>
+        error "diagonalMatrix: wrong number of entries in list"
+      diagonalMatrix(l)$Matrix(R) pretend $
+ 
+    coerce(x:$):Matrix(R) == copy(x pretend Matrix(R))
+ 
+    squareMatrix x ==
+      (nrows(x) ^= ndim) or (ncols(x) ^= ndim) =>
+        error "squareMatrix: matrix of bad dimensions"
+      copy(x) pretend $
+ 
+    x:$ * v:Col ==
+      directProduct((x pretend Matrix(R)) * (v :: Vector(R)))
+ 
+    v:Row * x:$ ==
+      directProduct((v :: Vector(R)) * (x pretend Matrix(R)))
+ 
+    x:$ ** n:NonNegativeInteger ==
+      ((x pretend Matrix(R)) ** n) pretend $
+ 
+    if R has commutative("*") then
+ 
+      determinant x == determinant(x pretend Matrix(R))
+      minordet x    == minordet(x pretend Matrix(R))
+ 
+    if R has EuclideanDomain then
+ 
+      rowEchelon x == rowEchelon(x pretend Matrix(R)) pretend $
+ 
+    if R has IntegralDomain then
+ 
+      rank x    == rank(x pretend Matrix(R))
+      nullity x == nullity(x pretend Matrix(R))
+      nullSpace x ==
+        [directProduct c for c in nullSpace(x pretend Matrix(R))]
+ 
+    if R has Field then
+ 
+      dimension() == (m * n) :: CardinalNumber
+ 
+      inverse x ==
+        (u := inverse(x pretend Matrix(R))) case "failed" => "failed"
+        (u :: Matrix(R)) pretend $
+ 
+      x:$ ** n:Integer ==
+        ((x pretend Matrix(R)) ** n) pretend $
+ 
+      recip x == inverse x
+ 
+    if R has ConvertibleTo InputForm then
+      convert(x:$):InputForm ==
+         convert [convert("squareMatrix"::Symbol)@InputForm,
+                  convert(x::Matrix(R))]$List(InputForm)
+
+
+@
+\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>>
+
+<<domain IMATRIX IndexedMatrix>>
+<<domain MATRIX Matrix>>
+<<domain RMATRIX RectangularMatrix>>
+<<domain SQMATRIX SquareMatrix>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}