Dear Martin, Waldek, and all.
I take an old matfuns.spad and I test this patch in order to insert theses
functions :
basis : List Vector Field -> List Vector Field
--- extracts a basis from a list of vectors
columnBasis : Matrix Fied -> List Vector Field
--- extracts a basis from the columns of a matrix
rowBasis : Matrix Fied -> List Vector Field
--- extracts a basis from the rows of a matrix
sumBasis : (List Vector Field, List Vector Field) -> List Vector Field
sumBasis : List List Vector Field -> List Vector Field
--- extracts a basis from the sum of subspaces
intBasis : (List Vector Field, List Vector Field) -> List Vector Field
intBasis : List List Vector Field -> List Vector Field
--- extracts a basis from the intersect of subspaces
This patch supposes that the basis of the null space is the empty list,
even if this old matfuns.spad file doesn't compute so.
I also send this short *.input file.
Problems remain :
1/ I don't know if this computation is possible over modules,
with Ring, not Field.
2/ I donif 't understand why nullSpace add a condition over these matrix :
shallowlyMutable. intBasis calls nullSpace, so I let this condition.
3/ nullSpace is possible over a Ring, I don't know what I can get for
theses basis functions over a Ring.
4/ The use of theses functions are heavy : see theses examples.input
Is it possible to have a better use ?
----------------------------------------------------------------------
--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.
-- This file and MATRIX SPAD must be compiled in bootstrap mode.
)abbrev package IMATLIN InnerMatrixLinearAlgebraFunctions
++ Author: Clifton J. Williamson, P.Gianni
++ Date Created: 13 November 1989
++ Date Last Updated: September 1993
++ Basic Operations:
++ Related Domains: IndexedMatrix(R,minRow,minCol), Matrix(R),
++ RectangularMatrix(n,m,R), SquareMatrix(n,R)
++ Also See:
++ AMS Classifications:
++ Keywords: matrix, canonical forms, linear algebra
++ Examples:
++ References:
++ Description:
++ \spadtype{InnerMatrixLinearAlgebraFunctions} is an internal package
++ which provides standard linear algebra functions on domains in
++ \spad{MatrixCategory}
InnerMatrixLinearAlgebraFunctions(R,Row,Col,M):_
Exports == Implementation where
R : Field
Row : FiniteLinearAggregate R
Col : FiniteLinearAggregate R
M : MatrixCategory(R,Row,Col)
I ==> Integer
Exports ==> with
rowEchelon: M -> M
++ \spad{rowEchelon(m)} returns the row echelon form of the matrix m.
rank: M -> NonNegativeInteger
++ \spad{rank(m)} returns the rank of the matrix m.
nullity: M -> NonNegativeInteger
++ \spad{nullity(m)} returns the mullity of the matrix m. This is the
++ dimension of the null space of the matrix m.
columnBasis: M -> List Col
++ \spad{basis m} returns a basis of the columns of the matrix m.
++ This basis is extracted from the columns of m.
rowBasis: M -> List Col
++ \spad{basis m} returns a basis of the rows of the matrix m.
++ This basis is extracted from the rows of m.
basis: List Col -> List Col
++ \spad{basis lv} extracts a basis of the subspace spanned by a
++ family of vectors.
sumBasis: (List Col, List Col) -> List Col
++ \spad{sumBasis (lv1,lv2)} extracts a basis of the sum of two
++ subspaces defined by two lists of vectors.
sumBasis: List List Col -> List Col
++ \spad{sumBasis llv} extracts a basis of the sum of a list
++ of subspaces defined by a list of lists of vectors.
if Col has shallowlyMutable then
nullSpace: M -> List Col
++ \spad{nullSpace(m)} returns a basis for the null space of the
++ matrix m.
intBasis: (List Col, List Col) -> List Col
++ \spad{intBasis (lv1,lv2)} extracts a basis of the intersect of two
++ subspaces defined by two lists of vectors.
intBasis: List List Col -> List Col
++ \spad{intBasis llv} extracts a basis of the intersect of a list
++ of subspaces defined by a list of lists of vectors.
determinant: M -> R
++ \spad{determinant(m)} returns the determinant of the matrix m.
++ an error message is returned if the matrix is not square.
generalizedInverse: M -> M
++ \spad{generalizedInverse(m)} returns the generalized (Moore--Penrose)
++ inverse of the matrix m, i.e. the matrix h such that
++ m*h*m=h, h*m*h=m, m*h and h*m are both symmetric matrices.
inverse: M -> Union(M,"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.
Implementation ==> add
rowAllZeroes?: (M,I) -> Boolean
rowAllZeroes?(x,i) ==
-- determines if the ith row of x consists only of zeroes
-- internal function: no check on index i
for j in minColIndex(x)..maxColIndex(x) repeat
qelt(x,i,j) ^= 0 => return false
true
colAllZeroes?: (M,I) -> Boolean
colAllZeroes?(x,j) ==
-- determines if the ith column of x consists only of zeroes
-- internal function: no check on index j
for i in minRowIndex(x)..maxRowIndex(x) repeat
qelt(x,i,j) ^= 0 => return false
true
rowEchelon y ==
-- row echelon form via Gaussian elimination
x := copy y
minR := minRowIndex x; maxR := maxRowIndex x
minC := minColIndex x; maxC := maxColIndex x
i := minR
n: I := minR - 1
for j in minC..maxC repeat
i > maxR => return x
n := minR - 1
-- n = smallest k such that k >= i and x(k,j) ^= 0
for k in i..maxR repeat
if qelt(x,k,j) ^= 0 then leave (n := k)
n = minR - 1 => "no non-zeroes"
-- put nth row in ith position
if i ^= n then swapRows_!(x,i,n)
-- divide ith row by its first non-zero entry
b := inv qelt(x,i,j)
qsetelt_!(x,i,j,1)
for k in (j+1)..maxC repeat qsetelt_!(x,i,k,b * qelt(x,i,k))
-- perform row operations so that jth column has only one 1
for k in minR..maxR repeat
if k ^= i and qelt(x,k,j) ^= 0 then
for k1 in (j+1)..maxC repeat
qsetelt_!(x,k,k1,qelt(x,k,k1) - qelt(x,k,j) * qelt(x,i,k1))
qsetelt_!(x,k,j,0)
-- increment i
i := i + 1
x
rank x ==
y :=
(rk := nrows x) > (rh := ncols x) =>
rk := rh
transpose x
copy x
y := rowEchelon y; i := maxRowIndex y
while rk > 0 and rowAllZeroes?(y,i) repeat
i := i - 1
rk := (rk - 1) :: NonNegativeInteger
rk :: NonNegativeInteger
sameSizeVectors? : List Col -> Boolean
sameSizeVectors2 : (Integer, List Col) -> Boolean
sameSizeVectors? Lv ==
null Lv => true
sameSizeVectors2 (#(first Lv), rest Lv)
sameSizeVectors2 (n,Lv) ==
null Lv => true
#(first Lv)=n => sameSizeVectors2 (n, rest Lv)
false
columnBasis mat ==
mat2 := rowEchelon mat
basis : List Col := []
indrow : Integer := 1
n : Integer := ncols mat
m : Integer := nrows mat
for k in 1..n repeat
if indrow <= m and mat2.(indrow,k) ^= 0
then
basis := cons (column (mat, k), basis)
indrow := indrow + 1
reverse basis
rowBasis mat == columnBasis transpose mat
basis Lv ==
null Lv => []
not (sameSizeVectors? Lv) => error "vectors have not the same size"
mm := new(#(first Lv),#Lv,0)
for k in 1..#Lv repeat setColumn! (mm, k, Lv.k)
columnBasis mm
sumBasis LLv == basis (reduce (concat, LLv))
sumBasis (Lv1, Lv2) == basis (concat (Lv1, Lv2))
nullity x == (ncols x - rank x) :: NonNegativeInteger
if Col has shallowlyMutable then
nullSpace y ==
x := rowEchelon y
minR := minRowIndex x; maxR := maxRowIndex x
minC := minColIndex x; maxC := maxColIndex x
nrow := nrows x; ncol := ncols x
basis : List Col := nil()
rk := nrow; row := maxR
-- compute rank = # rows - # rows of all zeroes
while rk > 0 and rowAllZeroes?(x,row) repeat
rk := (rk - 1) :: NonNegativeInteger
row := (row - 1) :: NonNegativeInteger
-- if maximal rank, return zero vector
ncol <= nrow and rk = ncol => [new(ncol,0)]
-- if rank = 0, return standard basis vectors
rk = 0 =>
for j in minC..maxC repeat
w : Col := new(ncol,0)
qsetelt_!(w,j,1)
basis := cons(w,basis)
basis
-- v contains information about initial 1's in the rows of x
-- if the ith row has an initial 1 in the jth column, then
-- v.j = i; v.j = minR - 1, otherwise
v : IndexedOneDimensionalArray(I,minC) := new(ncol,minR - 1)
for i in minR..(minR + rk - 1) repeat
for j in minC.. while qelt(x,i,j) = 0 repeat j
qsetelt_!(v,j,i)
j := maxC; l := minR + ncol - 1
while j >= minC repeat
w : Col := new(ncol,0)
-- if there is no row with an initial 1 in the jth column,
-- create a basis vector with a 1 in the jth row
if qelt(v,j) = minR - 1 then
colAllZeroes?(x,j) =>
qsetelt_!(w,l,1)
basis := cons(w,basis)
for k in minC..(j-1) for ll in minR..(l-1) repeat
if qelt(v,k) ^= minR - 1 then
qsetelt_!(w,ll,-qelt(x,qelt(v,k),j))
qsetelt_!(w,l,1)
basis := cons(w,basis)
j := j - 1; l := l - 1
basis
subVector (v:Col, a:Integer, b:Integer):Col ==
vv : Col := new ((b-a+1)::NonNegativeInteger, 0)
for k in 1..b-a+1 repeat vv.k := v.(k+a-1)
vv
linearSum (t:Col, Lv:List Col):Col ==
vv : Col := new (#(Lv.1),0)
for k in 1..#Lv repeat
v2 := Lv.k
t2 := t.k
for j in 1..#vv repeat
vv.j := vv.j + t2*v2.j
vv
-- reduce ("+", [t.i*Lv.i for i in 1..#t])
intBasis (Lv1, Lv2) ==
Lb1 := basis Lv1
Lb2 := basis Lv2
null Lb1 => []
null Lb2 => []
d1 := #Lb1
d2 := #Lb2
#(first Lb1) ^= #(first Lb2) => error "vectors have not the same size"
mm := new(#(first Lb1), d1+d2, 0)
for k in 1..d2 repeat setColumn! (mm, k, Lb2.k)
for k in 1..d1 repeat setColumn! (mm, d2+k, Lb1.k)
lkv := nullSpace mm
LcoeffV1 := [subVector (kv, d2+1, d1+d2) for kv in lkv]
[linearSum (cc, Lb1) for cc in LcoeffV1]
intBasis LLv ==
#LLv = 0 => error "no space to intersect"
#LLv = 1 => LLv.1 -- or reduce (intBasis, LLv)
intBasis (LLv.1, intBasis rest LLv)
determinant y ==
(ndim := nrows y) ^= (ncols y) =>
error "determinant: matrix must be square"
-- Gaussian Elimination
ndim = 1 => qelt(y,minRowIndex y,minColIndex y)
x := copy y
minR := minRowIndex x; maxR := maxRowIndex x
minC := minColIndex x; maxC := maxColIndex x
ans : R := 1
for i in minR..(maxR - 1) for j in minC..(maxC - 1) repeat
if qelt(x,i,j) = 0 then
rown := minR - 1
for k in (i+1)..maxR repeat
qelt(x,k,j) ^= 0 => leave (rown := k)
if rown = minR - 1 then return 0
swapRows_!(x,i,rown); ans := -ans
ans := qelt(x,i,j) * ans; b := -inv qelt(x,i,j)
for l in (j+1)..maxC repeat qsetelt_!(x,i,l,b * qelt(x,i,l))
for k in (i+1)..maxR repeat
if (b := qelt(x,k,j)) ^= 0 then
for l in (j+1)..maxC repeat
qsetelt_!(x,k,l,qelt(x,k,l) + b * qelt(x,i,l))
qelt(x,maxR,maxC) * ans
generalizedInverse(x) ==
SUP:=SparseUnivariatePolynomial R
FSUP := Fraction SUP
VFSUP := Vector FSUP
MATCAT2 := MatrixCategoryFunctions2(R, Row, Col, M,
FSUP, VFSUP, VFSUP, Matrix FSUP)
MATCAT22 := MatrixCategoryFunctions2(FSUP, VFSUP, VFSUP, Matrix FSUP,
R, Row, Col, M)
y:= map(coerce(coerce(#1)$SUP)$(Fraction SUP),x)$MATCAT2
ty:=transpose y
yy:=ty*y
nc:=ncols yy
var:=monomial(1,1)$SUP ::(Fraction SUP)
yy:=inverse(yy+scalarMatrix(ncols yy,var))::Matrix(FSUP)*ty
map(elt(#1,0),yy)$MATCAT22
inverse x ==
(ndim := nrows x) ^= (ncols x) =>
error "inverse: matrix must be square"
ndim = 2 =>
ans2 : M := zero(ndim, ndim)
zero?(det := x(1,1)*x(2,2)-x(1,2)*x(2,1)) => "failed"
detinv := inv det
ans2(1,1) := x(2,2)*detinv
ans2(1,2) := -x(1,2)*detinv
ans2(2,1) := -x(2,1)*detinv
ans2(2,2) := x(1,1)*detinv
ans2
AB : M := zero(ndim,ndim + ndim)
minR := minRowIndex x; maxR := maxRowIndex x
minC := minColIndex x; maxC := maxColIndex x
kmin := minRowIndex AB; kmax := kmin + ndim - 1
lmin := minColIndex AB; lmax := lmin + ndim - 1
for i in minR..maxR for k in kmin..kmax repeat
for j in minC..maxC for l in lmin..lmax repeat
qsetelt_!(AB,k,l,qelt(x,i,j))
qsetelt_!(AB,k,lmin + ndim + k - kmin,1)
AB := rowEchelon AB
elt(AB,kmax,lmax) = 0 => "failed"
subMatrix(AB,kmin,kmax,lmin + ndim,lmax + ndim)
)abbrev package MATCAT2 MatrixCategoryFunctions2
++ Author: Clifton J. Williamson
++ Date Created: 21 November 1989
++ Date Last Updated: 21 March 1994
++ Basic Operations:
++ Related Domains: IndexedMatrix(R,minRow,minCol), Matrix(R),
++ RectangularMatrix(n,m,R), SquareMatrix(n,R)
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Keywords: matrix, map, reduce
++ Examples:
++ References:
++ Description:
++ \spadtype{MatrixCategoryFunctions2} provides functions between two matrix
++ domains. The functions provided are \spadfun{map} and \spadfun{reduce}.
MatrixCategoryFunctions2(R1,Row1,Col1,M1,R2,Row2,Col2,M2):_
Exports == Implementation where
R1 : Ring
Row1 : FiniteLinearAggregate R1
Col1 : FiniteLinearAggregate R1
M1 : MatrixCategory(R1,Row1,Col1)
R2 : Ring
Row2 : FiniteLinearAggregate R2
Col2 : FiniteLinearAggregate R2
M2 : MatrixCategory(R2,Row2,Col2)
Exports ==> with
map: (R1 -> R2,M1) -> M2
++ \spad{map(f,m)} applies the function f to the elements of the matrix m.
map: (R1 -> Union(R2,"failed"),M1) -> Union(M2,"failed")
++ \spad{map(f,m)} applies the function f to the elements of the matrix m.
reduce: ((R1,R2) -> R2,M1,R2) -> R2
++ \spad{reduce(f,m,r)} returns a matrix n where
++ \spad{n[i,j] = f(m[i,j],r)} for all indices i and j.
Implementation ==> add
minr ==> minRowIndex
maxr ==> maxRowIndex
minc ==> minColIndex
maxc ==> maxColIndex
map(f:(R1->R2),m:M1):M2 ==
ans : M2 := new(nrows m,ncols m,0)
for i in minr(m)..maxr(m) for k in minr(ans)..maxr(ans) repeat
for j in minc(m)..maxc(m) for l in minc(ans)..maxc(ans) repeat
qsetelt_!(ans,k,l,f qelt(m,i,j))
ans
map(f:(R1 -> (Union(R2,"failed"))),m:M1):Union(M2,"failed") ==
ans : M2 := new(nrows m,ncols m,0)
for i in minr(m)..maxr(m) for k in minr(ans)..maxr(ans) repeat
for j in minc(m)..maxc(m) for l in minc(ans)..maxc(ans) repeat
(r := f qelt(m,i,j)) = "failed" => return "failed"
qsetelt_!(ans,k,l,r::R2)
ans
reduce(f,m,ident) ==
s := ident
for i in minr(m)..maxr(m) repeat
for j in minc(m)..maxc(m) repeat
s := f(qelt(m,i,j),s)
s
)abbrev package RMCAT2 RectangularMatrixCategoryFunctions2
++ Author: Clifton J. Williamson
++ Date Created: 21 November 1989
++ Date Last Updated: 12 June 1991
++ Basic Operations:
++ Related Domains: IndexedMatrix(R,minRow,minCol), Matrix(R),
++ RectangularMatrix(n,m,R), SquareMatrix(n,R)
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Keywords: matrix, map, reduce
++ Examples:
++ References:
++ Description:
++ \spadtype{RectangularMatrixCategoryFunctions2} provides functions between
++ two matrix domains. The functions provided are \spadfun{map} and
\spadfun{reduce}.
RectangularMatrixCategoryFunctions2(m,n,R1,Row1,Col1,M1,R2,Row2,Col2,M2):_
Exports == Implementation where
m,n : NonNegativeInteger
R1 : Ring
Row1 : DirectProductCategory(n, R1)
Col1 : DirectProductCategory(m, R1)
M1 : RectangularMatrixCategory(m,n,R1,Row1,Col1)
R2 : Ring
Row2 : DirectProductCategory(n, R2)
Col2 : DirectProductCategory(m, R2)
M2 : RectangularMatrixCategory(m,n,R2,Row2,Col2)
Exports ==> with
map: (R1 -> R2,M1) -> M2
++ \spad{map(f,m)} applies the function \spad{f} to the elements of the
++ matrix \spad{m}.
reduce: ((R1,R2) -> R2,M1,R2) -> R2
++ \spad{reduce(f,m,r)} returns a matrix \spad{n} where
++ \spad{n[i,j] = f(m[i,j],r)} for all indices spad{i} and \spad{j}.
Implementation ==> add
minr ==> minRowIndex
maxr ==> maxRowIndex
minc ==> minColIndex
maxc ==> maxColIndex
map(f,mat) ==
ans : M2 := new(m,n,0)$Matrix(R2) pretend M2
for i in minr(mat)..maxr(mat) for k in minr(ans)..maxr(ans) repeat
for j in minc(mat)..maxc(mat) for l in minc(ans)..maxc(ans) repeat
qsetelt_!(ans pretend Matrix R2,k,l,f qelt(mat,i,j))
ans
reduce(f,mat,ident) ==
s := ident
for i in minr(mat)..maxr(mat) repeat
for j in minc(mat)..maxc(mat) repeat
s := f(qelt(mat,i,j),s)
s
)abbrev package IMATQF InnerMatrixQuotientFieldFunctions
++ Author: Clifton J. Williamson
++ Date Created: 22 November 1989
++ Date Last Updated: 22 November 1989
++ Basic Operations:
++ Related Domains: IndexedMatrix(R,minRow,minCol), Matrix(R),
RectangularMatrix(n,m,R), SquareMatrix(n,R)
++ Also See:
++ AMS Classifications:
++ Keywords: matrix, inverse, integral domain
++ Examples:
++ References:
++ Description:
++ \spadtype{InnerMatrixQuotientFieldFunctions} provides functions on matrices
++ over an integral domain which involve the quotient field of that integral
++ domain. The functions rowEchelon and inverse return matrices with
++ entries in the quotient field.
InnerMatrixQuotientFieldFunctions(R,Row,Col,M,QF,Row2,Col2,M2):_
Exports == Implementation where
R : IntegralDomain
Row : FiniteLinearAggregate R
Col : FiniteLinearAggregate R
M : MatrixCategory(R,Row,Col)
QF : QuotientFieldCategory R
Row2 : FiniteLinearAggregate QF
Col2 : FiniteLinearAggregate QF
M2 : MatrixCategory(QF,Row2,Col2)
IMATLIN ==> InnerMatrixLinearAlgebraFunctions(QF,Row2,Col2,M2)
MATCAT2 ==> MatrixCategoryFunctions2(R,Row,Col,M,QF,Row2,Col2,M2)
CDEN ==> InnerCommonDenominator(R,QF,Col,Col2)
Exports ==> with
rowEchelon: M -> M2
++ \spad{rowEchelon(m)} returns the row echelon form of the matrix m.
++ the result will have entries in the quotient field.
inverse: M -> Union(M2,"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.
++ Note: the result will have entries in the quotient field.
if Col2 has shallowlyMutable then
nullSpace : M -> List Col
++ \spad{nullSpace(m)} returns a basis for the null space of the
++ matrix m.
Implementation ==> add
qfMat: M -> M2
qfMat m == map(#1 :: QF,m)$MATCAT2
rowEchelon m == rowEchelon(qfMat m)$IMATLIN
inverse m ==
(inv := inverse(qfMat m)$IMATLIN) case "failed" => "failed"
inv :: M2
if Col2 has shallowlyMutable then
nullSpace m ==
[clearDenominator(v)$CDEN for v in nullSpace(qfMat m)$IMATLIN]
)abbrev package MATLIN MatrixLinearAlgebraFunctions
++ Author: Clifton J. Williamson, P.Gianni
++ Date Created: 13 November 1989
++ Date Last Updated: December 1992
++ Basic Operations:
++ Related Domains: IndexedMatrix(R,minRow,minCol), Matrix(R),
++ RectangularMatrix(n,m,R), SquareMatrix(n,R)
++ Also See:
++ AMS Classifications:
++ Keywords: matrix, canonical forms, linear algebra
++ Examples:
++ References:
++ Description:
++ \spadtype{MatrixLinearAlgebraFunctions} provides functions to compute
++ inverses and canonical forms.
MatrixLinearAlgebraFunctions(R,Row,Col,M):Exports == Implementation where
R : CommutativeRing
Row : FiniteLinearAggregate R
Col : FiniteLinearAggregate R
M : MatrixCategory(R,Row,Col)
I ==> Integer
Exports ==> with
determinant: M -> R
++ \spad{determinant(m)} returns the determinant of the matrix m.
++ an error message is returned if the matrix is not square.
minordet: M -> R
++ \spad{minordet(m)} computes the determinant of the matrix m using
++ minors. Error: if the matrix is not square.
elRow1! : (M,I,I) -> M
++ elRow1!(m,i,j) swaps rows i and j of matrix m : elementary operation
++ of first kind
elRow2! : (M,R,I,I) -> M
++ elRow2!(m,a,i,j) adds to row i a*row(m,j) : elementary operation of
++ second kind. (i ^=j)
elColumn2! : (M,R,I,I) -> M
++ elColumn2!(m,a,i,j) adds to column i a*column(m,j) : elementary
++ operation of second kind. (i ^=j)
if R has IntegralDomain then
rank: M -> NonNegativeInteger
++ \spad{rank(m)} returns the rank of the matrix m.
nullity: M -> NonNegativeInteger
++ \spad{nullity(m)} returns the mullity of the matrix m. This is
++ the dimension of the null space of the matrix m.
nullSpace: M -> List Col
++ \spad{nullSpace(m)} returns a basis for the null space of the
++ matrix m.
fractionFreeGauss! : M -> M
++ \spad{fractionFreeGauss(m)} performs the fraction
++ free gaussian elimination on the matrix m.
invertIfCan : M -> Union(M,"failed")
++ \spad{invertIfCan(m)} returns the inverse of m over R
adjoint : M -> Record(adjMat:M, detMat:R)
++ \spad{adjoint(m)} returns the ajoint matrix of m (i.e. the matrix
++ n such that m*n = determinant(m)*id) and the detrminant of m.
if R has EuclideanDomain then
rowEchelon: M -> M
++ \spad{rowEchelon(m)} returns the row echelon form of the matrix m.
normalizedDivide: (R, R) -> Record(quotient:R, remainder:R)
++ normalizedDivide(n,d) returns a normalized quotient and
++ remainder such that consistently unique representatives
++ for the residue class are chosen, e.g. positive remainders
if R has Field then
inverse: M -> Union(M,"failed")
++ \spad{inverse(m)} returns the inverse of the matrix.
++ If the matrix is not invertible, "failed" is returned.
++ Error: if the matrix is not square.
columnBasis: M -> List Col
rowBasis: M -> List Col
basis: List Col -> List Col
sumBasis: (List Col, List Col) -> List Col
sumBasis: List List Col -> List Col
intBasis: (List Col, List Col) -> List Col
intBasis: List List Col -> List Col
Implementation ==> add
rowAllZeroes?: (M,I) -> Boolean
rowAllZeroes?(x,i) ==
-- determines if the ith row of x consists only of zeroes
-- internal function: no check on index i
for j in minColIndex(x)..maxColIndex(x) repeat
qelt(x,i,j) ^= 0 => return false
true
colAllZeroes?: (M,I) -> Boolean
colAllZeroes?(x,j) ==
-- determines if the ith column of x consists only of zeroes
-- internal function: no check on index j
for i in minRowIndex(x)..maxRowIndex(x) repeat
qelt(x,i,j) ^= 0 => return false
true
minorDet:(M,I,List I,I,PrimitiveArray(Union(R,"uncomputed")))-> R
minorDet(x,m,l,i,v) ==
z := v.m
z case R => z
ans : R := 0; rl : List I := nil()
j := first l; l := rest l; pos := true
minR := minRowIndex x; minC := minColIndex x;
repeat
if qelt(x,j + minR,i + minC) ^= 0 then
ans :=
md := minorDet(x,m - 2**(j :: NonNegativeInteger),_
concat_!(reverse rl,l),i + 1,v) *_
qelt(x,j + minR,i + minC)
pos => ans + md
ans - md
null l =>
v.m := ans
return ans
pos := not pos; rl := cons(j,rl); j := first l; l := rest l
minordet x ==
(ndim := nrows x) ^= (ncols x) =>
error "determinant: matrix must be square"
-- minor expansion with (s---loads of) memory
n1 : I := ndim - 1
v : PrimitiveArray(Union(R,"uncomputed")) :=
new((2**ndim - 1) :: NonNegativeInteger,"uncomputed")
minR := minRowIndex x; maxC := maxColIndex x
for i in 0..n1 repeat
qsetelt_!(v,(2**i - 1),qelt(x,i + minR,maxC))
minorDet(x, 2**ndim - 2, [i for i in 0..n1], 0, v)
-- elementary operation of first kind: exchange two rows --
elRow1!(m:M,i:I,j:I) : M ==
vec:=row(m,i)
setRow!(m,i,row(m,j))
setRow!(m,j,vec)
m
-- elementary operation of second kind: add to row i--
-- a*row j (i^=j) --
elRow2!(m : M,a:R,i:I,j:I) : M ==
vec:= map(a*#1,row(m,j))
vec:=map("+",row(m,i),vec)
setRow!(m,i,vec)
m
-- elementary operation of second kind: add to column i --
-- a*column j (i^=j) --
elColumn2!(m : M,a:R,i:I,j:I) : M ==
vec:= map(a*#1,column(m,j))
vec:=map("+",column(m,i),vec)
setColumn!(m,i,vec)
m
if R has IntegralDomain then
-- Fraction-Free Gaussian Elimination
fractionFreeGauss! x ==
(ndim := nrows x) = 1 => x
ans := b := 1$R
minR := minRowIndex x; maxR := maxRowIndex x
minC := minColIndex x; maxC := maxColIndex x
i := minR
for j in minC..maxC repeat
if qelt(x,i,j) = 0 then -- candidate for pivot = 0
rown := minR - 1
for k in (i+1)..maxR repeat
if qelt(x,k,j) ^= 0 then
rown := k -- found a pivot
leave
if rown > minR - 1 then
swapRows_!(x,i,rown)
ans := -ans
(c := qelt(x,i,j)) = 0 => "next j" -- try next column
for k in (i+1)..maxR repeat
if qelt(x,k,j) = 0 then
for l in (j+1)..maxC repeat
qsetelt_!(x,k,l,(c * qelt(x,k,l) exquo b) :: R)
else
pv := qelt(x,k,j)
qsetelt_!(x,k,j,0)
for l in (j+1)..maxC repeat
val := c * qelt(x,k,l) - pv * qelt(x,i,l)
qsetelt_!(x,k,l,(val exquo b) :: R)
b := c
(i := i+1)>maxR => leave
if ans=-1 then
lasti := i-1
for j in 1..maxC repeat x(lasti, j) := -x(lasti,j)
x
--
lastStep(x:M) : M ==
ndim := nrows x
minR := minRowIndex x; maxR := maxRowIndex x
minC := minColIndex x; maxC := minC+ndim -1
exCol:=maxColIndex x
det:=x(maxR,maxC)
maxR1:=maxR-1
maxC1:=maxC+1
minC1:=minC+1
iRow:=maxR
iCol:=maxC-1
for i in maxR1..1 by -1 repeat
for j in maxC1..exCol repeat
ss:=+/[x(i,iCol+k)*x(i+k,j) for k in 1..(maxR-i)]
x(i,j) := _exquo((det * x(i,j) - ss),x(i,iCol))::R
iCol:=iCol-1
subMatrix(x,minR,maxR,maxC1,exCol)
invertIfCan(y) ==
(nr:=nrows y) ^= (ncols y) =>
error "invertIfCan: matrix must be square"
adjRec := adjoint y
(den:=recip(adjRec.detMat)) case "failed" => "failed"
den::R * adjRec.adjMat
adjoint(y) ==
(nr:=nrows y) ^= (ncols y) => error "adjoint: matrix must be square"
maxR := maxRowIndex y
maxC := maxColIndex y
x := horizConcat(copy y,scalarMatrix(nr,1$R))
ffr:= fractionFreeGauss!(x)
det:=ffr(maxR,maxC)
[lastStep(ffr),det]
if R has Field then
VR ==> Vector R
IMATLIN ==> InnerMatrixLinearAlgebraFunctions(R,Row,Col,M)
MMATLIN ==> InnerMatrixLinearAlgebraFunctions(R,VR,VR,Matrix R)
FLA2 ==> FiniteLinearAggregateFunctions2(R, VR, R, Col)
MAT2 ==> MatrixCategoryFunctions2(R,Row,Col,M,R,VR,VR,Matrix R)
rowEchelon y == rowEchelon(y)$IMATLIN
rank y == rank(y)$IMATLIN
nullity y == nullity(y)$IMATLIN
determinant y == determinant(y)$IMATLIN
inverse y == inverse(y)$IMATLIN
basis lv == basis(lv)$IMATLIN
columnBasis m == columnBasis(m)$IMATLIN
rowBasis m == columnBasis(m)$IMATLIN
sumBasis (lv1,lv2) == sumBasis(lv1,lv2)$IMATLIN
sumBasis llv == sumBasis(llv)$IMATLIN
if Col has shallowlyMutable then
nullSpace y == nullSpace(y)$IMATLIN
intBasis (lv1,lv2) == intBasis(lv1,lv2)$IMATLIN
intBasis llv == intBasis(llv)$IMATLIN
else
nullSpace y ==
[map(#1, v)$FLA2 for v in nullSpace(map(#1, y)$MAT2)$MMATLIN]
else if R has IntegralDomain then
QF ==> Fraction R
Row2 ==> Vector QF
Col2 ==> Vector QF
M2 ==> Matrix QF
IMATQF ==> InnerMatrixQuotientFieldFunctions(R,Row,Col,M,QF,Row2,Col2,M2)
nullSpace m == nullSpace(m)$IMATQF
determinant y ==
(nrows y) ^= (ncols y) => error "determinant: matrix must be square"
fm:=fractionFreeGauss!(copy y)
fm(maxRowIndex fm,maxColIndex fm)
rank x ==
y :=
(rk := nrows x) > (rh := ncols x) =>
rk := rh
transpose x
copy x
y := fractionFreeGauss! y
i := maxRowIndex y
while rk > 0 and rowAllZeroes?(y,i) repeat
i := i - 1
rk := (rk - 1) :: NonNegativeInteger
rk :: NonNegativeInteger
nullity x == (ncols x - rank x) :: NonNegativeInteger
if R has EuclideanDomain then
if R has IntegerNumberSystem then
normalizedDivide(n:R, d:R):Record(quotient:R, remainder:R) ==
qr := divide(n, d)
qr.remainder >= 0 => qr
d > 0 =>
qr.remainder := qr.remainder + d
qr.quotient := qr.quotient - 1
qr
qr.remainder := qr.remainder - d
qr.quotient := qr.quotient + 1
qr
else
normalizedDivide(n:R, d:R):Record(quotient:R, remainder:R) ==
divide(n, d)
rowEchelon y ==
x := copy y
minR := minRowIndex x; maxR := maxRowIndex x
minC := minColIndex x; maxC := maxColIndex x
n := minR - 1
i := minR
for j in minC..maxC repeat
if i > maxR then leave x
n := minR - 1
xnj: R
for k in i..maxR repeat
if not zero?(xkj:=qelt(x,k,j)) and ((n = minR - 1) _
or sizeLess?(xkj,xnj)) then
n := k
xnj := xkj
n = minR - 1 => "next j"
swapRows_!(x,i,n)
for k in (i+1)..maxR repeat
qelt(x,k,j) = 0 => "next k"
aa := extendedEuclidean(qelt(x,i,j),qelt(x,k,j))
(a,b,d) := (aa.coef1,aa.coef2,aa.generator)
b1 := (qelt(x,i,j) exquo d) :: R
a1 := (qelt(x,k,j) exquo d) :: R
-- a*b1+a1*b = 1
for k1 in (j+1)..maxC repeat
val1 := a * qelt(x,i,k1) + b * qelt(x,k,k1)
val2 := -a1 * qelt(x,i,k1) + b1 * qelt(x,k,k1)
qsetelt_!(x,i,k1,val1); qsetelt_!(x,k,k1,val2)
qsetelt_!(x,i,j,d); qsetelt_!(x,k,j,0)
un := unitNormal qelt(x,i,j)
qsetelt_!(x,i,j,un.canonical)
if un.associate ^= 1 then for jj in (j+1)..maxC repeat
qsetelt_!(x,i,jj,un.associate * qelt(x,i,jj))
xij := qelt(x,i,j)
for k in minR..(i-1) repeat
qelt(x,k,j) = 0 => "next k"
qr := normalizedDivide(qelt(x,k,j), xij)
qsetelt_!(x,k,j,qr.remainder)
for k1 in (j+1)..maxC repeat
qsetelt_!(x,k,k1,qelt(x,k,k1) - qr.quotient * qelt(x,i,k1))
i := i + 1
x
else determinant x == minordet x
---------------------------------------------------------------------------
---
--- )sys rm -r MAT* IMAT* RMCAT2*
--- )co matfuns
--- )cd /home/fmy/Axiom
)lib MATLIN
)lib IMATLIN
)expose MATLIN
M := matrix [[1,2,3],[4,5,6],[7,8,9]]::Matrix Fraction Integer
basis M
QQ ==> Fraction Integer
VQ ==> Vector QQ
MQQ ==> Matrix QQ
(basis$MATLIN(QQ,VQ,VQ,MQQ)) M
Sp1 := nullSpace (matrix [row (M, 1)])
Sp2 := nullSpace (matrix [row (M, 2)])
Sp3 := nullSpace (matrix [row (M, 3)])
(sumBasis$MATLIN(QQ,VQ,VQ,MQQ))(Sp1, Sp2)
(sumBasis$MATLIN(QQ,VQ,VQ,MQQ))[Sp1, Sp2, Sp3]
(intBasis$MATLIN(QQ,VQ,VQ,MQQ))(Sp1, Sp2)
(intBasis$MATLIN(QQ,VQ,VQ,MQQ))[Sp1, Sp2, Sp3]
nullSpace M
-------------------------------------------------------------------------
This SF.net email is sponsored by: Microsoft
Defy all challenges. Microsoft(R) Visual Studio 2008.
http://clk.atdmt.com/MRT/go/vse0120000070mrt/direct/01/
_______________________________________________
open-axiom-devel mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/open-axiom-devel
