From ab8cc85adde879fb963c94d15675783f2cf4b183 Mon Sep 17 00:00:00 2001 From: dos-reis Date: Tue, 14 Aug 2007 05:14:52 +0000 Subject: Initial population. --- src/algebra/numeigen.spad.pamphlet | 413 +++++++++++++++++++++++++++++++++++++ 1 file changed, 413 insertions(+) create mode 100644 src/algebra/numeigen.spad.pamphlet (limited to 'src/algebra/numeigen.spad.pamphlet') diff --git a/src/algebra/numeigen.spad.pamphlet b/src/algebra/numeigen.spad.pamphlet new file mode 100644 index 00000000..310fe945 --- /dev/null +++ b/src/algebra/numeigen.spad.pamphlet @@ -0,0 +1,413 @@ +\documentclass{article} +\usepackage{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} -- cgit v1.2.3