Initial population.

1 files changed, 413 insertions, 0 deletions

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}
+<<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}
+<<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}
+<<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}
+<<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>>
+
+<<package INEP InnerNumericEigenPackage>>
+<<package NREP NumericRealEigenPackage>>
+<<package NCEP NumericComplexEigenPackage>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}