\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra numeigen.spad} \author{Patrizia Gianni} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package INEP InnerNumericEigenPackage} <>= )abbrev package INEP InnerNumericEigenPackage ++ Author:P. Gianni ++ Date Created: Summer 1990 ++ Date Last Updated:Spring 1991 ++ Basic Functions: ++ Related Constructors: ModularField ++ Also See: ++ AMS Classifications: ++ Keywords: ++ References: ++ Description: ++ This package is the inner package to be used by NumericRealEigenPackage ++ and NumericComplexEigenPackage for the computation of numeric ++ eigenvalues and eigenvectors. InnerNumericEigenPackage(K,F,Par) : C == T where F : Field -- this is the field where the answer will be -- for dealing with the complex case K : Field -- type of the input Par : Join(Field,OrderedRing) -- it will be NF or RN SE ==> Symbol() RN ==> Fraction Integer I ==> Integer NF ==> Float CF ==> Complex Float GRN ==> Complex RN GI ==> Complex Integer PI ==> PositiveInteger NNI ==> NonNegativeInteger MRN ==> Matrix RN MK ==> Matrix K PK ==> Polynomial K MF ==> Matrix F SUK ==> SparseUnivariatePolynomial K SUF ==> SparseUnivariatePolynomial F SUP ==> SparseUnivariatePolynomial MSUK ==> Matrix SUK PEigenForm ==> Record(algpol:SUK,almult:Integer,poleigen:List(MSUK)) outForm ==> Record(outval:F,outmult:Integer,outvect:List MF) IntForm ==> Union(outForm,PEigenForm) UFactor ==> (SUK -> Factored SUK) C == with charpol : MK -> SUK ++ charpol(m) computes the characteristic polynomial of a matrix ++ m with entries in K. ++ This function returns a polynomial ++ over K, while the general one (that is in EiegenPackage) returns ++ Fraction P K solve1 : (SUK, Par) -> List F ++ solve1(pol, eps) finds the roots of the univariate polynomial ++ polynomial pol to precision eps. If K is \spad{Fraction Integer} ++ then only the real roots are returned, if K is ++ \spad{Complex Fraction Integer} then all roots are found. innerEigenvectors : (MK,Par,UFactor) -> List(outForm) ++ innerEigenvectors(m,eps,factor) computes explicitly ++ the eigenvalues and the correspondent eigenvectors ++ of the matrix m. The parameter eps determines the type of ++ the output, factor is the univariate factorizer to br used ++ to reduce the characteristic polynomial into irreducible factors. T == add numeric(r:K):F == K is RN => F is NF => convert(r)$RN F is RN => r F is CF => r :: RN :: CF F is GRN => r::RN::GRN K is GRN => F is GRN => r F is CF => convert(convert r) error "unsupported coefficient type" ---- next functions neeeded for defining ModularField ---- monicize(f:SUK) : SUK == (a:=leadingCoefficient f) =1 => f inv(a)*f reduction(u:SUK,p:SUK):SUK == u rem p merge(p:SUK,q:SUK):Union(SUK,"failed") == p = q => p p = 0 => q q = 0 => p "failed" exactquo(u:SUK,v:SUK,p:SUK):Union(SUK,"failed") == val:=extendedEuclidean(v,p,u) val case "failed" => "failed" val.coef1 ---- eval a vector of F in a radical expression ---- evalvect(vect:MSUK,alg:F) : MF == n:=nrows vect w:MF:=zero(n,1)$MF for i in 1..n repeat polf:=map(numeric, vect(i,1))$UnivariatePolynomialCategoryFunctions2(K,SUK,F,SUF) v:F:=elt(polf,alg) setelt(w,i,1,v) w ---- internal function for the computation of eigenvectors ---- inteigen(A:MK,p:SUK,fact:UFactor) : List(IntForm) == dimA:NNI:= nrows A MM:=ModularField(SUK,SUK,reduction,merge,exactquo) AM:=Matrix(MM) lff:=factors fact(p) res: List IntForm :=[] lr : List MF:=[] for ff in lff repeat pol:SUK:= ff.factor if (degree pol)=1 then alpha:K:=-coefficient(pol,0)/leadingCoefficient pol -- compute the eigenvectors, rational case B1:MK := zero(dimA,dimA)$MK for i in 1..dimA repeat for j in 1..dimA repeat B1(i,j):=A(i,j) B1(i,i):= B1(i,i) - alpha lr:=[] for vecr in nullSpace B1 repeat wf:MF:=zero(dimA,1) for i in 1..dimA repeat wf(i,1):=numeric vecr.i lr:=cons(wf,lr) res:=cons([numeric alpha,ff.exponent,lr]$outForm,res) else ppol:=monicize pol alg:MM:= reduce(monomial(1,1),ppol) B:AM:= zero(dimA,dimA)$AM for i in 1..dimA repeat for j in 1..dimA repeat B(i,j):=reduce(A(i,j) ::SUK,ppol) B(i,i):=B(i,i) - alg sln2:=nullSpace B soln:List MSUK :=[] for vec in sln2 repeat wk:MSUK:=zero(dimA,1) for i in 1..dimA repeat wk(i,1):=(vec.i)::SUK soln:=cons(wk,soln) res:=cons([ff.factor,ff.exponent,soln]$PEigenForm, res) res if K is RN then solve1(up:SUK, eps:Par) : List(F) == denom := "lcm"/[denom(c::RN) for c in coefficients up] up:=denom*up upi := map(numer,up)$UnivariatePolynomialCategoryFunctions2(RN,SUP RN,I,SUP I) innerSolve1(upi, eps)$InnerNumericFloatSolvePackage(I,F,Par) else if K is GRN then solve1(up:SUK, eps:Par) : List(F) == denom := "lcm"/[lcm(denom real(c::GRN), denom imag(c::GRN)) for c in coefficients up] up:=denom*up upgi := map(complex(numer(real #1), numer(imag #1)), up)$UnivariatePolynomialCategoryFunctions2(GRN,SUP GRN,GI,SUP GI) innerSolve1(upgi, eps)$InnerNumericFloatSolvePackage(GI,F,Par) else error "unsupported matrix type" ---- the real eigenvectors expressed as floats ---- innerEigenvectors(A:MK,eps:Par,fact:UFactor) : List outForm == pol:= charpol A sln1:List(IntForm):=inteigen(A,pol,fact) n:=nrows A sln:List(outForm):=[] for lev in sln1 repeat lev case outForm => sln:=cons(lev,sln) leva:=lev::PEigenForm lval:List(F):= solve1(leva.algpol,eps) lvect:=leva.poleigen lmult:=leva.almult for alg in lval repeat nsl:=[alg,lmult,[evalvect(ep,alg) for ep in lvect]]$outForm sln:=cons(nsl,sln) sln charpol(A:MK) : SUK == dimA :PI := (nrows A):PI dimA ~= ncols A => error " The matrix is not square" B:Matrix SUK :=zero(dimA,dimA) for i in 1..dimA repeat for j in 1..dimA repeat B(i,j):=A(i,j)::SUK B(i,i) := B(i,i) - monomial(1,1)$SUK determinant B @ \section{package NREP NumericRealEigenPackage} <>= )abbrev package NREP NumericRealEigenPackage ++ Author:P. Gianni ++ Date Created:Summer 1990 ++ Date Last Updated:Spring 1991 ++ Basic Functions: ++ Related Constructors: FloatingRealPackage ++ Also See: ++ AMS Classifications: ++ Keywords: ++ References: ++ Description: ++ This package computes explicitly eigenvalues and eigenvectors of ++ matrices with entries over the Rational Numbers. ++ The results are expressed as floating numbers or as rational numbers ++ depending on the type of the parameter Par. NumericRealEigenPackage(Par) : C == T where Par : Join(Field,OrderedRing) -- Float or RationalNumber SE ==> Symbol() RN ==> Fraction Integer I ==> Integer NF ==> Float CF ==> Complex Float GRN ==> Complex RN GI ==> Complex Integer PI ==> PositiveInteger NNI ==> NonNegativeInteger MRN ==> Matrix RN MPar ==> Matrix Par outForm ==> Record(outval:Par,outmult:Integer,outvect:List MPar) C == with characteristicPolynomial : MRN -> Polynomial RN ++ characteristicPolynomial(m) returns the characteristic polynomial ++ of the matrix m expressed as polynomial ++ over RN with a new symbol as variable. -- while the function in EigenPackage returns Fraction P RN. characteristicPolynomial : (MRN,SE) -> Polynomial RN ++ characteristicPolynomial(m,x) returns the characteristic polynomial ++ of the matrix m expressed as polynomial ++ over RN with variable x. -- while the function in EigenPackage returns ++ Fraction P RN. realEigenvalues : (MRN,Par) -> List Par ++ realEigenvalues(m,eps) computes the eigenvalues of the matrix ++ m to precision eps. The eigenvalues are expressed as floats or ++ rational numbers depending on the type of eps (float or rational). realEigenvectors : (MRN,Par) -> List(outForm) ++ realEigenvectors(m,eps) returns a list of ++ records each one containing ++ a real eigenvalue, its algebraic multiplicity, and a list of ++ associated eigenvectors. All these results ++ are computed to precision eps as floats or rational ++ numbers depending on the type of eps . T == add import InnerNumericEigenPackage(RN, Par, Par) characteristicPolynomial(m:MRN) : Polynomial RN == x:SE:=new()$SE multivariate(charpol(m),x) ---- characteristic polynomial of a matrix A ---- characteristicPolynomial(A:MRN,x:SE):Polynomial RN == multivariate(charpol(A),x) realEigenvalues(m:MRN,eps:Par) : List Par == solve1(charpol m, eps) realEigenvectors(m:MRN,eps:Par) :List outForm == innerEigenvectors(m,eps,factor$GenUFactorize(RN)) @ \section{package NCEP NumericComplexEigenPackage} <>= )abbrev package NCEP NumericComplexEigenPackage ++ Author: P. Gianni ++ Date Created: Summer 1990 ++ Date Last Updated: Spring 1991 ++ Basic Functions: ++ Related Constructors: FloatingComplexPackage ++ Also See: ++ AMS Classifications: ++ Keywords: ++ References: ++ Description: ++ This package computes explicitly eigenvalues and eigenvectors of ++ matrices with entries over the complex rational numbers. ++ The results are expressed either as complex floating numbers or as ++ complex rational numbers ++ depending on the type of the precision parameter. NumericComplexEigenPackage(Par) : C == T where Par : Join(Field,OrderedRing) -- Float or RationalNumber SE ==> Symbol() RN ==> Fraction Integer I ==> Integer NF ==> Float CF ==> Complex Float GRN ==> Complex RN GI ==> Complex Integer PI ==> PositiveInteger NNI ==> NonNegativeInteger MRN ==> Matrix RN MCF ==> Matrix CF MGRN ==> Matrix GRN MCPar ==> Matrix Complex Par SUPGRN ==> SparseUnivariatePolynomial GRN outForm ==> Record(outval:Complex Par,outmult:Integer,outvect:List MCPar) C == with characteristicPolynomial : MGRN -> Polynomial GRN ++ characteristicPolynomial(m) returns the characteristic polynomial ++ of the matrix m expressed as polynomial ++ over complex rationals with a new symbol as variable. -- while the function in EigenPackage returns Fraction P GRN. characteristicPolynomial : (MGRN,SE) -> Polynomial GRN ++ characteristicPolynomial(m,x) returns the characteristic polynomial ++ of the matrix m expressed as polynomial ++ over Complex Rationals with variable x. -- while the function in EigenPackage returns Fraction P GRN. complexEigenvalues : (MGRN,Par) -> List Complex Par ++ complexEigenvalues(m,eps) computes the eigenvalues of the matrix ++ m to precision eps. The eigenvalues are expressed as complex floats or ++ complex rational numbers depending on the type of eps (float or rational). complexEigenvectors : (MGRN,Par) -> List(outForm) ++ complexEigenvectors(m,eps) returns a list of ++ records each one containing ++ a complex eigenvalue, its algebraic multiplicity, and a list of ++ associated eigenvectors. All these results ++ are computed to precision eps and are expressed as complex floats ++ or complex rational numbers depending on the type of eps (float or rational). T == add import InnerNumericEigenPackage(GRN,Complex Par,Par) characteristicPolynomial(m:MGRN) : Polynomial GRN == x:SE:=new()$SE multivariate(charpol m, x) ---- characteristic polynomial of a matrix A ---- characteristicPolynomial(A:MGRN,x:SE):Polynomial GRN == multivariate(charpol A, x) complexEigenvalues(m:MGRN,eps:Par) : List Complex Par == solve1(charpol m, eps) complexEigenvectors(m:MGRN,eps:Par) :List outForm == innerEigenvectors(m,eps,factor$ComplexFactorization(RN,SUPGRN)) @ \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}