\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra smith.spad} \author{Patrizia Gianni} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package SMITH SmithNormalForm} <>= )abbrev package SMITH SmithNormalForm ++ Author: Patrizia Gianni ++ Date Created: October 1992 ++ Date Last Updated: ++ Basic Operations: ++ Related Domains: Matrix(R) ++ Also See: ++ AMS Classifications: ++ Keywords: matrix, canonical forms, linear algebra ++ Examples: ++ References: ++ Description: ++ \spadtype{SmithNormalForm} is a package ++ which provides some standard canonical forms for matrices. SmithNormalForm(R,Row,Col,M) : Exports == Implementation where R : EuclideanDomain Row : FiniteLinearAggregate R Col : FiniteLinearAggregate R M : MatrixCategory(R,Row,Col) I ==> Integer NNI ==> NonNegativeInteger HermiteForm ==> Record(Hermite:M,eqMat:M) SmithForm ==> Record(Smith : M, leftEqMat : M, rightEqMat : M) PartialV ==> Union(Col, "failed") Both ==> Record(particular: PartialV, basis: List Col) Exports == with hermite : M -> M ++ \spad{hermite(m)} returns the Hermite normal form of the ++ matrix m. completeHermite : M -> HermiteForm ++ \spad{completeHermite} returns a record that contains ++ the Hermite normal form H of the matrix and the equivalence matrix ++ U such that U*m = H smith : M -> M ++ \spad{smith(m)} returns the Smith Normal form of the matrix m. completeSmith : M -> SmithForm ++ \spad{completeSmith} returns a record that contains ++ the Smith normal form H of the matrix and the left and right ++ equivalence matrices U and V such that U*m*v = H diophantineSystem : (M,Col) -> Both ++ \spad{diophantineSystem(A,B)} returns a particular integer solution and ++ an integer basis of the equation \spad{AX = B}. Implementation == add MATCAT1 ==> MatrixCategoryFunctions2(R,Row,Col,M,QF,Row2,Col2,M2) MATCAT2 ==> MatrixCategoryFunctions2(QF,Row2,Col2,M2,R,Row,Col,M) QF ==> Fraction R Row2 ==> Vector QF Col2 ==> Vector QF M2 ==> Matrix QF ------ Local Functions ----- elRow1 : (M,I,I) -> M elRow2 : (M,R,I,I) -> M elColumn2 : (M,R,I,I) -> M isDiagonal? : M -> Boolean ijDivide : (SmithForm ,I,I) -> SmithForm lastStep : SmithForm -> SmithForm test1 : (M,Col,NNI) -> Union(NNI, "failed") test2 : (M, Col,NNI,NNI) -> Union( Col, "failed") -- inconsistent system : case 0 = c -- test1(sm:M,b:Col,m1 : NNI) : Union(NNI , "failed") == km:=m1 while zero? sm(km,km) repeat if not zero?(b(km)) then return "failed" km:= (km - 1) :: NNI km if Col has ShallowlyMutableAggregate R then test2(sm : M ,b : Col, n1:NNI,dk:NNI) : Union( Col, "failed") == -- test divisibility -- sol:Col := new(n1,0) for k in 1..dk repeat if (c:=(b(k) exquo sm(k,k))) case "failed" then return "failed" sol(k):= c::R sol -- test if the matrix is diagonal or pseudo-diagonal -- isDiagonal?(m : M) : Boolean == m1:= nrows m n1:= ncols m for i in 1..m1 repeat for j in 1..n1 | (j ~= i) repeat if not zero?(m(i,j)) then return false true -- 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 -- modify SmithForm in such a way that the term m(i,i) -- -- divides the term m(j,j). m is diagonal -- ijDivide(sf : SmithForm , i : I,j : I) : SmithForm == m:=sf.Smith mii:=m(i,i) mjj:=m(j,j) extGcd:=extendedEuclidean(mii,mjj) d := extGcd.generator mii:=(mii exquo d)::R mjj := (mjj exquo d) :: R -- add to row j extGcd.coef1*row i -- lMat:=elRow2(sf.leftEqMat,extGcd.coef1,j,i) -- switch rows i and j -- lMat:=elRow1(lMat,i,j) -- add to row j -mii*row i -- lMat := elRow2(lMat,-mii,j,i) -- lMat := ijModify(mii,mjj,extGcd.coef1,extGcd.coef2,sf.leftEqMat,i,j) m(j,j):= m(i,i) * mjj m(i,i):= d -- add to column i extGcd.coef2 * column j -- rMat := elColumn2(sf.rightEqMat,extGcd.coef2,i,j) -- add to column j -mjj*column i -- rMat:=elColumn2(rMat,-mjj,j,i) -- multiply by -1 column j -- setColumn!(rMat,j,map(-1 * #1,column(rMat,j))) [m,lMat,rMat] -- given a diagonal matrix compute its Smith form -- lastStep(sf : SmithForm) : SmithForm == m:=sf.Smith m1:=min(nrows m,ncols m) for i in 1..m1 while not zero?(mii:=m(i,i)) repeat for j in i+1..m1 repeat if (m(j,j) exquo mii) case "failed" then return lastStep(ijDivide(sf,i,j)) sf -- given m and t row-equivalent matrices, with t in upper triangular -- -- form compute the matrix u such that u*m=t -- findEqMat(m : M,t : M) : Record(Hermite : M, eqMat : M) == m1:=nrows m n1:=ncols m "and"/[zero? t(m1,j) for j in 1..n1] => -- there are 0 rows if "and"/[zero? t(1,j) for j in 1..n1] then return [m,scalarMatrix(m1,1)] -- m is the zero matrix mm:=horizConcat(m,scalarMatrix(m1,1)) mmh:=rowEchelon mm [subMatrix(mmh,1,m1,1,n1), subMatrix(mmh,1,m1,n1+1,n1+m1)] u:M:=zero(m1,m1) j: Integer :=1 while t(1,j)=0 repeat j:=j+1 -- there are 0 columns t1:=copy t mm:=copy m if j>1 then t1:=subMatrix(t,1,m1,j,n1) mm:=subMatrix(m,1,m1,j,n1) t11:=t1(1,1) for i in 1..m1 repeat u(i,1) := (mm(i,1) exquo t11) :: R for j: local in 2..m1 repeat j0:=j tjj : R while zero?(tjj:=t1(j,j0)) repeat j0:=j0+1 u(i,j) :=((mm(i,j0) - ("+"/[u(i,k) * t1(k,j0) for k in 1..(j-1)])) exquo tjj) :: R u1:M2:= map(#1 :: QF,u)$MATCAT1 [t,map(retract$QF,(inverse u1)::M2)$MATCAT2] --- Hermite normal form of m --- hermite(m:M) : M == rowEchelon m -- Hermite normal form and equivalence matrix -- completeHermite(m : M) : Record(Hermite : M, eqMat : M) == findEqMat(m,rowEchelon m) smith(m : M) : M == completeSmith(m).Smith completeSmith(m : M) : Record(Smith : M, leftEqMat : M, rightEqMat : M) == cm1:=completeHermite m leftm:=cm1.eqMat m1:=cm1.Hermite isDiagonal? m1 => lastStep([m1,leftm,scalarMatrix(ncols m,1)]) nr:=nrows m cm1:=completeHermite transpose m1 rightm:= transpose cm1.eqMat m1:=cm1.Hermite isDiagonal? m1 => cm2:=lastStep([m1,leftm,rightm]) nrows(m:=cm2.Smith) = nr => cm2 [transpose m,cm2.leftEqMat, cm2.rightEqMat] cm2:=completeSmith m1 cm2:=lastStep([cm2.Smith,transpose(cm2.rightEqMat)*leftm, rightm*transpose(cm2.leftEqMat)]) nrows(m:=cm2.Smith) = nr => cm2 [transpose m, cm2.leftEqMat, cm2.rightEqMat] -- Find the solution in R of the linear system mX = b -- diophantineSystem(m : M, b : Col) : Both == sf:=completeSmith m sm:=sf.Smith m1:=nrows sm lm:=sf.leftEqMat b1:Col:= lm* b (t1:=test1(sm,b1,m1)) case "failed" => ["failed",empty()] dk:=t1 :: NNI n1:=ncols sm (t2:=test2(sm,b1,n1,dk)) case "failed" => ["failed",empty()] rm := sf.rightEqMat sol:=rm*(t2 :: Col) -- particular solution dk = n1 => [sol,list new(n1,0)] lsol:List Col := [column(rm,i) for i in (dk+1)..n1] [sol,lsol] @ \section{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. @ <<*>>= <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}